ODE45 Returning Wrong Signed Answer

Hey, I have been struggling with this for the past few days. Basically, I have an equation that models the angular acceleration of a rigid body (pendulum). I want to integrate it twice (to get the plot of angle vs time). The angle should be decreasing, however, in my plot it is increasing. I can't seem to figure out why. The equation is correct, I suspect there is something I have implemented wrong.

Here is a sample of my code, first I have the function in its own file, then I use the main code with the input main. If you have any idea why it is producing the wrong answer please let me know

function dydt = phiddot
g=9.81; L=1; m=1;
dydt = [y(2); ((-6*g*sin(y(1)))/(L*(1+3*(sin(y(1))^2))))];
end

function maincode = main
t0=0;
T=10;
 y0 = [60;0];  
[ts,ys] = ode45(@phi_ddot,[t0,T],y0); 
plot(ts,ys(:,1),'r')

end

Answer 1

[David Goodmanson]

Hi h7,

I'm guessing that you want to start at 60 degrees and not 60 radians, so try

y0 = [60*pi/180;0];  
[ts,ys] = ode45(@phi_ddot,[t0,T],y0); 
plot(ts,ys(:,1)*180/pi,'r')   % plot in degrees

which does start out going negative.

In your example your function needed to be called phi_ddot in order to match the ode45 call.

[Steven Lord]

Instead of converting between degrees and radians, maybe the ODE function phiddot or phi_ddot should use the degree-based sind instead of the radian-based sin.

[hussein.j]

David! You are a lifesavor! I overlooked such a small detail! Thank you so much!!!

Stay safe!

[hussein.j]

Steven, I will also give that a try, didn't know about sind! Thanks

[David Goodmanson]

The problem with that idea is that now the size of the derivatives will be different.

 dsin(theta)/dt = cos(theta)              theta in radians, but 
 dsind(theta)/dt = (pi/180) cosd(theta)   theta in degrees  

so you have to make a correction anyway. Without a correction, you can see that happen if you go to 60 degrees as the starting point and use sind, in which case there are far fewer oscillations.

Since the differential eqn is set up for radians, it's safer to convert to radians at the beginning and to degrees after the calculation is over (if desired). I usually I define a constant C = pi/180 at the beginning and use that, since it's easier.

Nothing against sind and cosd in general, they are very useful most of the time such as with rotation matrices.