How to generate double pendulum using ode 45

function Math462hw2Q3parta 
m1=1;
m2=1;
l1=1;
l2=1;
g=9.8;
tspan=0:1:100;
%theta1initial=pi/8;
%theta2initial=pi/8;
theta1=pi/8;
theta2=pi/8;
theta1dot=0;
theta2dot=0;
thetainitial=[theta1 theta2 theta1dot theta2dot];

%Write out the ode45 function
[t,thetasoln]=ode45(@(t,thetainitial)derivativefuncs(t,thetainitial,m1,m2,l1,l2,g),tspan,thetainitial);

%Plot function
plot(t,thetasoln,'bl');
xlabel('time');
ylabel('value');
title('double pendulum');

    function dxdt=derivativefuncs(~,thetavalues,m1,m2,l1,l2,g)
        theta1=thetavalues(1);
        theta2=thetavalues(2);
        theta1dot=thetavalues(3);
        theta2dot=thetavalues(4);
        
        theta1doubledot=(m2*l2*theta2doubledot*cos(theta1-theta2)+m2*l2*theta2dot^2*sin(theta1-theta2)+(m1+m2)*g*sin(theta1))/((m1+m2)*l);
        theta2doubledot=l1*theta1doubledot*cos(theta1-theta2)-l1*theta1dot^2*sin(theta1-theta2)+g*sin(theta2);
        dxdt=[theta1doubledot;theta2doubledot];
        
    end 
end

Just got another question on how to use ode 45 to generate the double pendulum model. I think the problem is Matlab keeps saying the line with ode45 function is still wrong but there is no marks in the code. I think I did pass the parameters but it still does not work. Can some one maybe help?

Answer 1

[James Tursa]

You have a 4-element state vector, so your derivative needs to be a 4-element state vector. E.g.,

dxdt=[theta1dot;theta2dot;theta1doubledot;theta2doubledot];

Then you will need to pick off the appropriate columns of thetasoln to plot.

[Zhukun Wang]

Yeah, but it is still saying the line with ode45 is not right. The following is my lastest version.

function Math462hw2Q3parta

m1=1;

m2=1;

l1=1;

l2=1;

g=9.8;

tspan=0:1:100;

theta1=pi/8;

theta2=pi/8;

theta1dot=0;

theta2dot=0;

thetainitial=[theta1 theta2 theta1dot theta2dot];

%Write out the ode45 function

[t,thetasoln]=ode45(@(t,thetainitial)derivativefuncs(t,thetainitial,m1,m2,l1,l2,g),tspan,thetainitial);

%Plot function

plot(t,thetasoln(:,1),'bl');

xlabel('time');

ylabel('value');

title('double pendulum');

function dxdt=derivativefuncs(~,thetavalues,m1,m2,l1,l2,g)

theta1=thetavalues(1);

theta2=thetavalues(2);

theta1dot=thetavalues(3);

theta2dot=thetavalues(4);

theta1doubledot=(m2*l2*theta2doubledot*cos(theta1-theta2)+m2*l2*theta2dot^2*sin(theta1-theta2)+(m1+m2)*g*sin(theta1))/((m1+m2)*l);

theta2doubledot=l1*theta1doubledot*cos(theta1-theta2)-l1*theta1dot^2*sin(theta1-theta2)+g*sin(theta2);

dxdt=[theta1dot;theta2dot;theta1doubledot;theta2doubledot];

end

end

[James Tursa]

"not right" doesn't tell us much. Can you copy & post the entire error message including the offending line?

[Zhukun Wang]

出错 Math462hw2Q3parta (第 19 行)

[t,thetasoln]=ode45(@(t,thetainitial)derivativefuncs(t,thetainitial,m1,m2,l1,l2,g),tspan,thetainitial);

This is in Chinese by the way. The character you see only means error. Matlab does not tell me much about the error from the message.