Why does the numerical solution of an ODE move away from an unstable equilibrium?

Question

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. 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.

ODE45

When I run the simulation for 6 seconds with ode45, I get the expected solution with no problems at all, the system stays where it is and does not move:

However, when I run the simulation for 10 seconds, the system starts moving unreasonably:

ODE23

I then ran the simulation with ode23 to see if the problem persisted. I end up with the same behavior, only this time the divergence begins 1 second later:

ODE15s

I then ran the simulation with ode15s to see if the issue persisted and no, the system appears to be stable even during 100 seconds:

Then again, ode15s is only first order and note that there are only a few integrating steps. So I ran yet another simulation with ode15s during 10 seconds but a MaxStep size of $0.01$ to increase precision, and unfortunately, this leads to the same outcome as with both ode45 and ode23.

Normally, the obvious outcome of these simulations would be that the system stays at its initial position since nothing is perturbing it. Why is this divergence occurring? Does it have something to do with the fact that these type of systems are chaotic in nature? Is this a normal behavior for ode functions in Matlab?

Answer 1

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in.

To give a bit more detail how this plays out, consider your problem in the form of a general initial value problem

\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{M}^{-1} \mathbf{f}(t, \mathbf{y}(t)) \end{equation} where $\mathbf{y}(t) = (x_1(t), x_2(t), x_3(t), x_4(t))$ and

\begin{equation} \mathbf{f}(t, \mathbf{y}(t)) = \begin{bmatrix} x_2 \\ -V_1 - G_1 \\ x_4 \\ -V_2 - G_2 \end{bmatrix} \end{equation}

In exact arithmetic, you would have $\mathbf{f}(0, \mathbf{y}_0) = 0$ and thus $\dot{\mathbf{y}}(0) = 0 $ and nothing ever changes: your system stays in equilibrium. However, the arithmetic in a computer are not exact, meaning that for a variety of reasons, your right hand side is not exactly zero but equal to some $\tilde{\mathbf{f}}$ which is almost but not quite zero.

In your code, you can test that by computing

norm(rhs(0,[pi/2 0 0 0]))

which gives 6.191e-16 - so almost but not exactly zero. How does that impact the dynamics of your system?

Under some assumptions, the effect of $\mathbf{f}$ not being precisely zero is the same is if you would not start from the initial value $\mathbf{y}_0$ you prescribe but from a value that is very slightly different, let us call it $\tilde{\mathbf{y}}_0$.

Furthermore, over a very short time, the solution to your system looks like the solution of the linearised system

\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{f}(0, \mathbf{y}_0) + \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) = \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) \end{equation}

where $\mathbf{f}'$ is the Jacobian of your function $\mathbf{f}$ or rhs in your code. Since $\mathbf{y}_0$ is constant, we can recast this into an equation for $\mathbf{d}(t) := \mathbf{y}(t) - \mathbf{y}_0$, where $\mathbf{d}$ says how far away from the initial value we are:

\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{f}'(0, \mathbf{y}_0) \mathbf{d}(t). \end{equation}

I couldn't be bothered to compute the Jacobian by hand so I used automatic differentiation to get a good approximation:

\begin{equation} \mathbf{J} := \mathbf{f}'(0, \mathbf{y}_0) = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 9.81 & 0 & 2.4525 & 0 \\ 0 & 0 & 0 & 1 \\ 2.4525 & 0 & 2.4525 & 0 \end{bmatrix} \end{equation}

so that your equation becomes

\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{J} \mathbf{d}(t), \mathbf{d}(0) = \mathbf{\tilde{y}}_0 - \mathbf{y}_0 \end{equation}

Now we need one final step: we can compute an eigenvalue decomposition of the Jacobian such that

\begin{equation} \mathbf{J} = \mathbf{Q} \mathbf{D} \mathbf{Q}^{-1} \end{equation}

where $\mathbf{D}$ is a diagonal matrix containing the eigenvalues of $\mathbf{J}$ and $\mathbf{Q}$ are orthogonal matrices, representing coordinate transformations. We can then transform the equation for $\mathbf{d}$ into an equation for $\mathbf{e}(t) := \mathbf{Q}^{-1} \mathbf{d}(t)$ which reads

\begin{equation} \dot{\mathbf{e}}(t) = \mathbf{D} \mathbf{e}(t), \mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}_0. \end{equation}

Because $\mathbf{D}$ is diagonal, these are effectively four independent equations

\begin{equation} \dot{e}_i(t) = \lambda_i e_i(t), e_i(0) = i-\textrm{th component of}~\mathbf{Q}^{-1} \mathbf{d}_0 \end{equation}

with $i=1,2,3,4$. If you compute the eigenvalues, you'll find that the largest one is $\lambda_1 = 3.2485$. Therefore,

\begin{equation} e_1(t) = e_1(0) e^{3.2485 t}. \end{equation}

Now if arithmetic in your computer were exact, you'd have $\mathbf{d}(0) = 0$, thus $\mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}(0) = 0$ and thus $e_1(0) = 0$ and nothing would happen. But since this is not the case, you have a small but finite $e_1(0)$ which gets exponentially amplified. Hence the rapid deviation from the equilibrium in your solution.

Answer 2

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.

Answer 3

Look at the components of the forces calculated in your functions.

You will probably find they are never exactly zero, because as other answers have said, you can't represent the value of $\pi$ exactly in computer arithmetic, and the routines that calculate trig functions are not exact either.

Eventually, the tiny forces (probably of order $10^{-16}$ at the start) will move the system away from its unstable equilibrium position.

While the displacement of the system is still very small, all the calculations will lose a lot of precision through rounding errors (you are doing calculations similar to $a = 1.0$; $ a = a + 10^{-16}$) so the amount of time before the system "topples over" in the simulation will depend on exactly what integration method you used, what time steps you requested, etc.

Answer 4

The initial assumption was that the initial position was at a stable equilibrium (i.e., a minimum of the potential energy) with zero kinetic energy and the system started moving away from the equilibrium.
Since physically it can't happen (if we consider classical mechanics), two things came to my mind:

1. The first one is that maybe the initial position is: both pendulums pointing upwards ($\pi/2$ instead of $-\pi/2$?), which is a point of unstable equilibrium;

2. The second one is that perhaps there is something wrong with the equations of movement (perhaps a typo somewhere?). Can you please write the equations explicitly? Perhaps you could plot the angular acceleration as a function of the initial position of each pendulum, assuming zero angular velocity to check if there is something weird.

Answer 5

You should search more about double pendulums: they are what we call "chaotic systems". Even though they behave following simple rules, starting from slightly different initial conditions, solutions diverge quite fast. Doing numerical simulations for this kind of systems is not easy. Take a look at the following video to get more insight into the problem.

For the simple or double pendulum you could write a formula for the total energy of the system. Assuming that friction forces are neglected, this total energy is preserved by the analytic system. Numerically this is a whole other issue.

Before trying the double pendulum, try the simple pendulum. You will notice that for Runge-Kutta methods of small order the energy of the system will grow in the numerical simulations, instead of remaining constant (this is what happens in your simulations: you get movement out of nothing). In order to prevent this, higher order RK methods could be used (ode45 is of order 4; RK of order 8 would work better). There are other methods called "symplectic methods" which are designed such that the numerical simulations conserve the energy. In general you should stop the simulation as soon as the energy increases significantly compared to your initialization.

And try to understand the simple pendulum before going to the double one.