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

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?

• Besides the equations, I think the a schematic would also help a lot to understand b the question. – nicoguaro Jun 21 '19 at 3:06
• If you think it's appropriate, you may accept one of the answers (there is a green button). – Ertxiem - reinstate Monica Jun 21 '19 at 7:57
• You don't say, but you seem to be plotting x1 and x3. (Insert dry comment about graphs without legends or descriptions.) Try plotting the logarithms of (the absolute values of) x2 and x4. – Eric Towers Jun 21 '19 at 8:33

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

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

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

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

$$$$\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)$$$$

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:

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

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

$$$$\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}$$$$

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

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

$$$$\mathbf{J} = \mathbf{Q} \mathbf{D} \mathbf{Q}^{-1}$$$$

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

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

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

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

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

$$$$e_1(t) = e_1(0) e^{3.2485 t}.$$$$

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.

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.

• If you monitor the state variables carefully (by looking at the values printed out in scientific notation), you should be able to see the initial very slow movement away from equilibrium. – Brian Borchers Jun 21 '19 at 4:15
• This makes sense and indeed, when I start the system in a downwards vertical position (being a point of stable equilibrium), the system does not move at all, at least for a simulation of 1000 seconds which I consider a very long period of time. – Meclassic Jun 21 '19 at 4:50
• If you want to solve the problem, change your parametrization so that $x_1$ is the angle with respect to the vertical axis. Then in your equations you'll need to compute sin(0) and cos(0) which are exact, rather than sin(pi/2) and cos(pi/2). In this way I think you can guarantee that rhs(t,[0,0,0 0] == [0,0,0,0], which means your system won't move at all. – Federico Poloni Jun 21 '19 at 8:30
• @jrojasqu just to make explicitly clear what has been said: your variables are represented by double precision floats, which can only represent a finite set of numbers, and $\pi/2$ is not one of them for several reasons, the simplest being that it's irrational and so has an infinite number of digits, but floats have a finite number of digits. Zero is a member of this set and so it can be exactly represented. – llama Jun 21 '19 at 19:34
• Note that if the straight down direction is not represented by $\theta = 0$, you should also see some movement where you would mathematically expect none, but in this case it will just be oscillating about the stable equilibrium with an amplitude on the order of the float precision ($~10^{-16}$ for doubles) – llama Jun 21 '19 at 19:41

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.

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.

• Indeed, I started the system in an upwards vertical position. Therefore, it is a point of unstable equilibrium. Brian Borcher's comment completes your answer by explaining the issue with $\pi$ approximation which leads to the system eventually moving from that position. – Meclassic Jun 21 '19 at 4:45
• By the way, just for fun, if you wanted to keep the system at the unstable vertical position, you could change your origin of coordinates to have the angle equal to zero pointing upwards. – Ertxiem - reinstate Monica Jun 21 '19 at 7:56
• @Ertxiem another option is to introduce small friction in the pins that would eat numerical errors. – svavil Jun 22 '19 at 1:04
• @Ertxiem For the sake of fun, I tried changing the coordinate system so that the zero angle makes the system point upwards. It is really the best parameterization here. Obviously the system stays in the upwards position indefinitely. However oscillations still arise (minimal for 1000 seconds of simulation but there) at the stable equilibrium position (straight downwards) because then, there is a $\sin(\pi)$ to be calculated in the force deriving from the potential energy. So I insist in the fact that if I simulate this long enough, the system will start deviating also from that position. – Meclassic Jun 22 '19 at 5:52
• Since physically it can't happen – Given the insight that we are at an unstable equilibrium, I somewhat challenge this. Physical systems (without too much friction) do not stay in unstable equilibria. More generally, if you simulate real systems, you would want to avoid that it accidentally gets stuck in an unstable equilibrium (however it got there) – it’s a feature, not a bug. (There are some rare exceptions to this, such as the uninfected state in immunology, which is an unstable equilibrium that can be maintained.) – Wrzlprmft Jun 22 '19 at 18:38

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.

• This is not a matter of the system being chaotic, though. You can have an unstable equilibrium in non-chaotic systems as well, e.g. single pendulum being "on its head", and it will exhibit the same behaviour described in the question. – Daniel Jun 21 '19 at 12:07
• It is also not true that the energy increases for RKM of small order: implicit Euler ist first order and shows exactly the opposite behaviour. – Daniel Jun 21 '19 at 12:08
• @BeniBogosel You mention symplectic methods which caught my attention because obviously, in my example, the energy is not conserved. However, could you indicate a specific symplectic method that could be implemented here? – Meclassic Jun 21 '19 at 15:15
• @jrojasqu why do you say that energy is not conserved on your system? – Ertxiem - reinstate Monica Jun 22 '19 at 0:18
• @Ertxiem When I calculate the total Mechanical Energy of the system (Kinetic+Potential energies) with the outputs provided by ode45, I get a value that starts with zero, then grows with time. The value is very very small at the first seconds, but still, it consistently grows away from zero. I believe it is because of the issues that were addressed in the answers above (approximation of $\pi$, etc.). – Meclassic Jun 22 '19 at 2:21