Returning original mathematical function values from ode45 in Matlab?

Question

My intention is to plot the original mathematical function values from the differential equation of the second order below:

I(thetadbldot)+md(g-o^2asin(ot))sin(theta)=0

where thetadbldot is the second derivative of theta with respect to t and m,d,I,g,a,o are given constants. Initial conditions are theta(0)=pi/2 and thetadot(0)=0.

My issue is that my knowledge and tutoring is limited to storing the values of the derivatives and returning them, not values from the original mathematical function in the equation. Below you can see a code that calculates the differential in Cauchy-form and gives me the derivatives. Does anyone have suggestions what to do? Thanks!

function xdot = penduluma(t,x)
% The function penduluma(t,x) calculates the differential 
% I(thetadbldot)+md(g-o^2asin(ot))sin(theta)=0 where thetadbldot is the second 
% derivative of theta with respect to t and m,d,I,g,a,o are given constants.
% For the state-variable form, x1=theta and x2=thetadot. x is a 2x1 vector on the form
% [theta,thetadot].
m=1;d=0.2;I=0.1;g=9.81;a=0.1;o=4;
xdot = [x(2);m*d*(o^2*a*sin(o*t)-g)*sin(x(1))/I];
end

options=odeset('RelTol', 1e-6);
[t,xa]=ode45(@penduluma,[0,20],[pi/2,0],options);
% Then the desired vector from xa is plotted to t. As it looks now the desired 
% values are not found in xa however.

There might be a misunderstanding from my side about what is returned from ode45 as xa. I am starting to lean towards it being x1 and x2, not x1dot and x2dot. — Christian

am304 am304 · Accepted Answer · 2013-12-16T12:25:43

Once you have the angle, you can calculate the angular velocity and acceleration using diff:

options=odeset('RelTol', 1e-6);
[t,xa]=ode45(@penduluma,[0,20],[pi/2,0],options);
x_ddot = zeros(size(t));
x_ddot(2:end) = diff(xa(:,2))./diff(t);
plot(t,xa,t,x_ddot)
legend('angle','angular velocity','angular acceleration')

which gives the following plot in Octave (should be the same in MATLAB):

enter image description here

Alternatively, you can work it out using your original differential equation:

x_ddot = -m*d*(o^2*a*sin(o*t)-g).*sin(xa(:,1))/I;

which gives a similar result:

enter image description here

Returning original mathematical function values from ode45 in Matlab?

1 Answers