Using ODE45 with 4 states on a double particle pendulum

Hello! I am trying to simulate a double particle pendulum with matlab code. I have already found the EOM's and put them in state-space form but I am having trouble solving my equations of motions which include two second order ODE's. I have tried using ode45 but I keep on running into problems trying to use it. After solving the ODE with ode45 I plan to plot the results then simulate them. I will paste my code below and any help will be much appriciated! For referance I would like to solve d2theta1 and d2theta2 in my code with ode45.

The error message is:

Error using odearguments (line 113)

Inputs must be floats, namely single or double.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in AE149_Project1 (line 30)

[t,y] = ode45(@(t,y) odeFunc, t_span, theta_0)

And the code is:

clear all; close all; clc;
%% Double particle pendulm 
%% Define variables and constants 
m1 = 30;
L1 = 10;
m2 = 30;
L2 = 10;
g = 9.81; % Gravity
% Intitial Conditions
theta_0 = [0];

%Time span
t_span = linspace(0,20,600);

% EOMs

syms theta1(t) theta2(t);
% % x = zeros(1,4);
% 
x1 = theta1;
x2 = diff(theta1);
x3 = theta2;
x4 = diff(theta2);

d2theta1 = (m2*L2*x4^2*sin(x3-x1)-(m1+m2)*g*sin(x1)+m2*cos(x3-x1)*(L1*x2^2*sin(x3-x1)+g*sin(x3)))/(L1*(m1+m2)-m2*L1*(cos(x3-x1)^2));
d2theta2 = (-cos(x3-x1)*(m2*L2*x4^2*sin(x3-x1)-(m1+m2)*g*sin(x1))+(m1+m2)*(-L1*x2^2*sin(x3-x1)-g*sin(x3)))/(L2*(m1+m2)-m2*L2*cos(x3-x1)^2);

odeFunc = matlabFunction(d2theta2)
% Solving ODE
[t,y] = ode45(@(t,y) odeFunc, t_span, theta_0)

Answer 1

[Paul]

Hi Raaed,

Why use symbolic when you alredy have a closed form expressions for the differential equations? Usually, the Symbolic Math Toolbox is used to derive these expressions. But you don't need to derive them, so no need to use the symbolic stuff. Besides that, there were a couple of issues in the code related to the initial condition, the definition of odeFunc, and the call to ode45.

odeFunn needs to take two inputs, t and x, where x in this case is 4x1 vector, and return a 4x1 output of derivatives. The call to ode45 should not include the @(t,y), which is used to define the function, not to use at as function handle input argument to ode45.

Alternatively, it might be clearer to define odeFunc as an m-function, as opposed to an anonymous function as done here.

Anyway, check out the modified code below, which at least runs w/o error ann I think implments the intended equations.

%% Double particle pendulm 
%% Define variables and constants 
m1 = 30;
L1 = 10;
m2 = 30;
L2 = 10;
g = 9.81; % Gravity
% Intitial Conditions
theta_0 = [0]; % need four initial conditions
%Time span
t_span = linspace(0,20,600);
% EOMs
%syms theta1 theta2;
x0 = [5 0 5 0]*pi/180; % need four initial conditions
% 
% x(1) = theta1;
% x(2) = diff(theta1);
% x(3) = theta2;
% x(4) = diff(theta2);
% define second derivatives as functions of x
d2theta1 = @(x) (m2*L2*x(4)^2*sin(x(3)-x(1))-(m1+m2)*g*sin(x(1))+m2*cos(x(3)-x(1))*(L1*x(2)^2*sin(x(3)-x(1))+g*sin(x(3))))/(L1*(m1+m2)-m2*L1*(cos(x(3)-x(1))^2));
d2theta2 = @(x) (-cos(x(3)-x(1))*(m2*L2*x(4)^2*sin(x(3)-x(1))-(m1+m2)*g*sin(x(1)))+(m1+m2)*(-L1*x(2)^2*sin(x(3)-x(1))-g*sin(x(3))))/(L2*(m1+m2)-m2*L2*cos(x(3)-x(1))^2);

% odeFunc needs to return a 4x1 vector xdot as a function of t and 4x1 vector x
odeFunc = @(t,x) ([x(2) ; d2theta1(x) ; x(4) ; d2theta2(x)]);

% just call ode45 with odeFunc
[t,y] = ode45(odeFunc, t_span, x0);
plot(t,y),grid