What is the flaw in my stability analysis?

Question

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?

Answer 1

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)[0])
                x2.append(np.matmul(C,X0)[1])

        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:

Answer 2

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

Answer 3

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.