ODE45 function in Matlab diverges from obvious solution, is it normal?
I wish to simulate the behaviour of a double-pendulum-like system. The system is a 2-degrees-of-freedom robot manipulator that is not actuated and will, therefore, behave mostly like a double-pendulum affected by gravity. 935
The only main difference with a double-pendulum is that it is composed of two rigid bodies with mass and inertia properties at their centers of mass.
Basically, I programmed ode45 under Matlab to solve a system of ODEs of the following type:
$$ \left[ \begin{array}{cccc} 1 & 0 & 0 & 0\ 0 & M_{11} & 0 & M_{12}\ 0 & 0 & 1 & 0\ 0 & M_{12} & 0 & M_{22} \end{array} \right] \left[ \begin{array}{c} \dot{x}_1\ \dot{x}_2\ \dot{x}_3\ \dot{x}_4 \end{array} \right]= \left[ \begin{array}{c} x_2\ -V_1-G_1\ x_4\ -V_2-G_2 \end{array} \right] $$
where $x_1$ is the angle of the first body with respect to the horizontal, $x_2$ is the angular velocity of the first body; $x_3$ is the angle of the second body with respect to the first body, and $x_4$ is the angular velocity of the second body. All of the coefficients are specified in the following code, in the rhs and fMass functions I created.
clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);
function F=rhs(t,x)
m=[1 1];
l=0.5;
a=[0.25 0.25];
g=9.81;
c1=cos(x(1));
s2=sin(x(3));
c12=cos(x(1)+x(3));
n1=m(2)*a(2)*l;
V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
V2=n1*s2*x(2)^2;
G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
G2=m(2)*g*a(2)*c12;
F(1)=x(2);
F(2)=-V1-G1;
F(3)=x(4);
F(4)=-V2-G2;
F=F';
end
function M=fMass(t,x)
m=[1 1];
l=0.5;
Izz=[0.11 0.11];
a=[0.25 0.25];
c2=cos(x(3));
n1=m(2)*a(2)*l;
M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
M12=m(2)*a(2)^2+n1*c2+Izz(2);
M22=m(2)*a(2)^2+Izz(2);
M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end
Notice how I set the initial condition of $x_1$ (angle of the first body with respect to the horizontal) so that the system starts in a completely vertical position. This way, since only gravity is acting, the obvious outcome is that the system should not move at all from that position.
NOTE: in all of the graphics below, I plotted the solutions $x_1$ and $x_3$ with respect to time. 804