Using Implicit Euler with second order differential equations

Question

We can numerically integrate first order differential equations using Euler method like this: $$y_{n+1} = y_n + hf(t_n, y_n)$$

And with Implicit Euler like this: $$y_{n+1} = y_n + hf(t_{n+1},y _{n+1})$$

If I have a differential equation $y' - ky = 0$, I can integrate $y$ numerically using Implicit Euler: $$y_{n+1} = y_n + hky_{n+1}$$ $$y_{n+1} = y_n\frac{1}{1-hk}$$

But how I do use Implicit Euler for second order differential equations, like for instance the equation for simple harmonic motion? $$y'' + w^2y = 0$$

We have to integrate with respect to $y$ and $y'$. For explicit Euler the numerical integration would look like this (?): $$y_{n+1} = y_n + hf(t_n, y'_n)$$ $$y'_{n+1} = y'_n + hg(t_n, y_n)$$

How would we do integrate using Implicit Euler instead?

Answer 1

You have to write your second order equation as a system of two first order equations. Let $y' = v$, then your equation

$$ y'' + \omega^2 y = 0 $$

becomes

$$ \begin{pmatrix} y' \\ v' \end{pmatrix} = \begin{pmatrix} v \\ -\omega^2 y \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \begin{pmatrix} y \\ v \end{pmatrix} $$

If you denote $u = (y,v)$ and the matrix as $A$, this system can be written as

$$ u' = A u $$

Now implicit Euler reads

$$ u_{n+1} = u_n + h A u_{n+1} $$

or

$$ \left( I - h A \right) u_{n+1} = u_n $$

Translated back into components you get

$$ \begin{pmatrix} 1 & -h \\ h \omega^2 & 1 \end{pmatrix} \begin{pmatrix} y_{n+1} \\ v_{n+1} \end{pmatrix} = \begin{pmatrix} y_{n} \\ v_n \end{pmatrix} $$

Addendum. It might be worth pointing out that implicit Euler is not a very good integrator for this type of problem as it will lead to artificial energy dissipation. You might be better of with what is called symplectic Euler method.