Solving system of second order differential equations with ode45

I have two second equations

theta1''= [-29.4sin(theta1)- 9.8sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(theta1'^(2)*cos(theta1-theta2)]/(3-cos(2*theta1-2*theta2))

theta2''=[2*sin(theta1-theta2)[(2*theta1')+19.6cos(theta1)+theta2'*cos(theta1-theta2)]]/(3-cos(2*theta1-2*theta1))

I think these should be written as a system of 4 first order equations, recast as a matrix and put into ode45 but I cannot figure out hwo to write these equatuons as 4 first first order due to the trig functions. I think I know to use ode45 once I have them in this form. Any help on how to rewrite these would be appreaciated.

[Thomas Holmes]

These theta1 and theta2 values represent the angles between the vertical and two point masses on a double pendulum, how can I use theta1 and theta2 to plot the path the bottom pendulum takes?

[Star Strider]

The angles will be the ‘y’ output of your ode45 call:

[t,y] = ode45(odesfcn, tspan, theta0);

With those and the geometry of your system, you can plot the trajectories of the point masses. You would have to model the masses as functions of the angles.

That is the best I can do in terms of a description.

Answer 1

[Star Strider]

I prefer to let the Symbolic MAth Toolbox do the ‘heavy lifting’:

syms theta1(t) theta2(t) t Y 
% theta1'' = -29.4sin(theta1)- 9.8sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(theta1'^(2)*cos(theta1-theta2)]/(3-cos(2*theta1-2*theta2))
% theta2'' = 2*sin(theta1-theta2)[(2*theta1')+19.6cos(theta1)+theta2'*cos(theta1-theta2)]]/(3-cos(2*theta1-2*theta1))
Dt11 = diff(theta1);
Dt12 = diff(Dt11);
Dt21 = diff(theta2);
Dt22 = diff(Dt21);
Eqn1 = Dt12 == (-29.4*sin(theta1)- 9.8*sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(Dt11^(2)*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta2));
Eqn2 = Dt22 == (2*sin(theta1-theta2)*((2*theta1')+19.6*cos(theta1)+Dt21*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta1));
Eqns = [Eqn1; Eqn2];
[VF,Sbs] = odeToVectorField(Eqns);
Sbs
odesfcn = matlabFunction(VF, 'Vars',{t, Y})

producing:

Sbs =
 
  theta2
 Dtheta2
  theta1
 Dtheta1
 

odesfcn =

  function_handle with value:

    @(t,Y)[Y(2);-sin(Y(1)-Y(3)).*(conj(Y(3)).*2.0+cos(Y(3)).*(9.8e+1./5.0)+cos(Y(1)-Y(3)).*Y(2));Y(4);-(sin(Y(1).*2.0-Y(3)).*(4.9e+1./5.0)-sin(Y(3)).*(1.47e+2./5.0)+cos(Y(1)-Y(3)).*sin(Y(1)-Y(3)).*Y(4).^2.*2.0)./(cos(Y(1).*2.0-Y(3).*2.0)-3.0)]

or:

    odesfcn = @(t,Y)[Y(2);-sin(Y(1)-Y(3)).*(conj(Y(3)).*2.0+cos(Y(3)).*(9.8e+1./5.0)+cos(Y(1)-Y(3)).*Y(2));Y(4);-(sin(Y(1).*2.0-Y(3)).*(4.9e+1./5.0)-sin(Y(3)).*(1.47e+2./5.0)+cos(Y(1)-Y(3)).*sin(Y(1)-Y(3)).*Y(4).^2.*2.0)./(cos(Y(1).*2.0-Y(3).*2.0)-3.0)];

Use the ‘Sbs’ vector in your legend call.

NOTE — I did not test this with an actual ode45 call. You may need to experiment with it to get it to work correctly.

[Thomas Holmes]

I dont understand what this code produces for the ode45 command. I tried saving the code as a function and then calling it in the ode45 command.

function dydt=Sbs(t,y)
syms theta1(t) theta2(t) t Y 
% theta1'' = -29.4sin(theta1)- 9.8sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(theta1'^(2)*cos(theta1-theta2)]/(3-cos(2*theta1-2*theta2))
% theta2'' = 2*sin(theta1-theta2)[(2*theta1')+19.6cos(theta1)+theta2'*cos(theta1-theta2)]]/(3-cos(2*theta1-2*theta1))
Dt11 = diff(theta1);
Dt12 = diff(Dt11);
Dt21 = diff(theta2);
Dt22 = diff(Dt21);
Eqn1 = Dt12 == (-29.4*sin(theta1)- 9.8*sin(theta1 -2*theta2) - 2*sin(theta1-theta2)*(Dt11^(2)*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta2));
Eqn2 = Dt22 == (2*sin(theta1-theta2)*((2*theta1')+19.6*cos(theta1)+Dt21*cos(theta1-theta2)))/(3-cos(2*theta1-2*theta1));
Eqns = [Eqn1; Eqn2];
[VF,Sbs] = odeToVectorField(Eqns);
Sbs
odesfcn = matlabFunction(VF, 'Vars',{t, Y}); 

Then i called this function in the ode45 command

clear all; close all; clc
opts=odeset('OutputFcn','odephas2');
[t,y] = ode45(@Sbs,[0 50],[pi/8,0,pi/8,0],opts)

I just get an error saying Sbs is a variable not a function. I feel like I am not understanding the code you provided

[Star Strider]

Use this with ode45:

odesfcn = @(t,Y)[Y(2);-sin(Y(1)-Y(3)).*(conj(Y(3)).*2.0+cos(Y(3)).*(9.8e+1./5.0)+cos(Y(1)-Y(3)).*Y(2));Y(4);-(sin(Y(1).*2.0-Y(3)).*(4.9e+1./5.0)-sin(Y(3)).*(1.47e+2./5.0)+cos(Y(1)-Y(3)).*sin(Y(1)-Y(3)).*Y(4).^2.*2.0)./(cos(Y(1).*2.0-Y(3).*2.0)-3.0)];

so:

[t,y] = ode45(odesfcn, tspan, theta0);

where ‘theta0’ is your vector of initial conditions.

The ‘Sbs’ output simply tells you that ‘Y(1) is theta2’ and so for the others.

[Thomas Holmes]

Thank you!

[Star Strider]

As always, my pleasure!