# What is the flaw in my stability analysis?

The ODE $${d^2x\over dt^2}=-kx; k>0$$can be converted in the system of linear equations as \begin{align} {dx\over dt} & =v\\ {dv\over dt} &= -kx\\ \end{align}

Using Euler’s method, given $$x_n$$ and $$y_n$$ and for the time step $$\Delta t$$, the next values can be determined as $$\left[ \begin{matrix} x_{n+1}\\ v_{n+1}\\ \end{matrix}\right] = \left[\begin{matrix} 1&\Delta t\\ -k\Delta t&1 \end{matrix}\right] \left[\begin{matrix} x_n\\ v_n\\ \end{matrix}\right].$$

Now the absolute value of the (possibly complex) eigenvalues should be less than $$1$$ for this algorithm to be stable. But the eigenvalues turn out to be $$1\pm i\sqrt{k}\Delta t$$ whose absolute values are strictly greater than $$1$$ for any nonzero time-step $$\Delta t$$.

So the algorithm should not work for any value of $$\Delta t$$, however small. But clearly, this is not the case as my programs do come up with (an approximate) solution though.

So where is the flaw in my reasoning?

• Its hard to tell when you have a pretty obvious typo in the statement. Is your ODE in fact $$\frac{d^2x}{dt^2} = -kx$$. please clarify so we can help.
– EMP
Oct 17 '19 at 16:53
• Yep, sorry! I’ll edit.
– Atom
Oct 17 '19 at 16:56
• For future readers: the system is purely oscillatory (imaginary eigenvalues) and conserves its energy. No fixed-step explicit schemes allow for a correct long term solution (they either explode or dissipate). You can write the Hamiltonian and use a sympletic integrator so that the numerical solution respect the energy conservation. More simply, you can use the Crank-Nicolson scheme which has the property that it does not dissipate or amplify pure complex eigenvalues, therefore the oscillation will keep a constant amplitude (but its phase may gradually degrade if the time step is too large). Nov 22 '20 at 15:29

But clearly, this is not the case as my programs do come up with (an approximate) solution though.

I believe you did not continue the integration until you see that your integration is not convergent and is not bounded.

I could rewrite your system of ODEs as:

$$\dot{x_{1}} = x_{2}$$

$$\dot{x_{2}} = -kx_{1}$$

Or in matrix form:

$$\dot{X} = AX$$

Where: $$X = \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix}$$ and $$A = \begin{bmatrix} 0 & 1 \\ -k & 0 \end{bmatrix}$$

The equilibrium point of your system of ODEs is: $$X^{*} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$, but this equilibrium point is unstable cause the real part of eigenvalues of $$A$$ are not all negative: $$\lambda_{1} = i\sqrt{k}$$ and $$\lambda_{2} = -i\sqrt{k}$$. In fact, the real part of eigenvalues are zero for these two eigenvalues. So, the conclusion is: no matter how you choose a small $$\Delta t$$, the forward Euler integration will not remain bounded.

Let's look at your discretization. I could discretize this system of ODEs in matrix form as:

$$X_{n+1} = (I+\Delta t A) X_{n}$$

Where $$X_{n+1}$$ and $$X_{n}$$ are $$X$$ vectors at times $$n+1$$ and $$n$$ respectively. The general formula is:

$$X_{n} = (I+\Delta t A)^{n} X_{0}$$

Where $$X_{0}$$ is initial condition for vector $$X$$. In order to have a bounded solution, I need to make sure the Frobenius norm of $$||I+\Delta t A||_{F} < 1$$. But we have:

$$||I+\Delta t A||_{F} = \sqrt{2+(1+k^{2})\Delta t^{2}} > 1$$

Which shows that no matter what you choose for $$\Delta t$$, if you continue the integration long enough, finally $$(I+\Delta t A)^{n}$$ will be blown up at some point.

This is the implementation with Python:

import numpy as np
import matplotlib.pyplot as plt

k = 1

deltats = np.linspace(0.01,0.1,5)

A = [[0,1],[-k,0]]
I = [[1,0],[0,1]]

A = np.array(A)
I = np.array(I)

X0 = [0,np.sqrt(k)]

X0 = np.array(X0)

for deltat in deltats:
x1 = []
x2 = []
B = I + deltat * A
ts = np.linspace(0,100,int(100/deltat))
for i,t in enumerate(ts):
C = np.linalg.matrix_power(B,i)
x1.append(np.matmul(C,X0))
x2.append(np.matmul(C,X0))

plt.plot(ts,x1,label=r'$$x_{1}$$, $$\Delta t$$ = '+str(deltat))
#plt.plot(ts,x2,label=r'$$x_{2}$$, $$\Delta t$$ = '+str(deltat))

plt.xlabel('t')
plt.ylabel('X')
plt.legend(loc='best')
plt.show()


And you see, when we expect the solution of this system of ODEs with initial condition of $$X_{0} = \begin{bmatrix} 0 \\ \sqrt{k} \end{bmatrix}$$ to be $$X(t) = \begin{bmatrix} \sin(\sqrt{k}t) \\ \sqrt{k}\cos(\sqrt{k}t) \end{bmatrix}$$ and clearly the solution should be bounded smaller than 1, but you see it's not bounded when you continue the integration long enough: This problem has an invariant which is the total energy $$E(t) = \frac{1}{2}(\dot{x}^2 + k x^2) = \textrm{constant}$$ As done by AloneProgrammer, write as first order system $$\dot{x}_1 = x_2, \qquad \dot{x}_2 = - k x_1$$ In the phase space $$(x_1,x_2)$$, the solution must stay on an ellipse whose size is determined by the initial energy.

Applying forward Euler to this, you can show that $$E^{n+1} = E^n + \frac{1}{2}k (\Delta t)^2[ k (x_1^n)^2 + (x_2^n)^2] > E^n$$ no matter what $$\Delta t$$ you choose. The solution spirals out in phase space.

Using backward Euler, you can show that $$E^{n+1} < E^n$$ for any $$\Delta t$$. The solution spirals inward in phase space.

Trapezoidal method would conserve energy, see

https://nbviewer.jupyter.org/github/cpraveen/na/blob/master/ode/periodic1.ipynb

For such problems, look for symplectic methods, e.g.

https://www.unige.ch/~hairer/poly_geoint/week2.pdf

So I'm not too familiar with this form of analysis, but it looks to me like the issue is that your $$x_{n+1}$$ update should be using $$v_{n+1}$$ rather than $$v_n$$ but I could be wrong.

• No, in forward-Euler method, as I’ve used, it is as I’ve mentioned.
– Atom
Oct 17 '19 at 17:04