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?

  • $\begingroup$ 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. $\endgroup$
    – EMP
    Oct 17, 2019 at 16:53
  • $\begingroup$ Yep, sorry! I’ll edit. $\endgroup$
    – Atom
    Oct 17, 2019 at 16:56
  • $\begingroup$ 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). $\endgroup$
    – Laurent90
    Nov 22, 2020 at 15:29

3 Answers 3


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)

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


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:

enter image description here


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


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



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.

  • 3
    $\begingroup$ No, in forward-Euler method, as I’ve used, it is as I’ve mentioned. $\endgroup$
    – Atom
    Oct 17, 2019 at 17:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.