# Non-conservative implementation implicit Euler

In Matlab R2013a I have implemented the Implicit Euler (time) integration scheme. To find the $x^{n+1}$ value I use fixed point iterations:

$x^{n+1} = \Delta t f(x^{n+1}) + x^n$

To test this, I use an undamped mass-spring system, which should give a simple harmonic as result. However, I seem to loose a considerable amount of energy over time.

Is this due to the fixed point iterations finding a 'wrong' fixed point?

Thanks!

• Dissipation is to be expected for long times. Try with en.wikipedia.org/wiki/Semi-implicit_Euler_method or some other symplectic integrator. Sep 4 '14 at 13:24
• I can agree with some dissipation, but after two periods the amplitude is nearly halved. This seems quite extreme to me. A smaller timestep seems to dissipate less. Sep 4 '14 at 13:36

This might seem extreme, but this can be analysed exactly. Take the system $$\dot x_1 = x_2, \qquad \dot x_2=-x_1, \qquad x_1(0) = 1, \qquad x_2(0)=0.$$ Let $X=(x_1,x_2)$ be the state vector, $\delta t$ the time step, and $X^+$ the state vector for the next time step. Then the implicit Euler scheme is $$X^+ = \delta t\left(\begin{array}{cc}0&1\\-1&0\end{array}\right) X^+ + X,$$ or $$X^+ = \frac{1}{1+\delta t^2}\begin{pmatrix}1&\delta t\\-\delta t&1\end{pmatrix}X = AX.$$ Now, the conserved quantity is $\|X\|_2$, so take the matrix, call it $A$, that multiplies $X$ every time step, and consider how it changes the norm of $X$. Its 2-norm is $$\|A\|_2 = \frac{1}{(1+\delta t^2)^{1/2}},$$ so if there are $n$ steps in a period of length $2\pi$, the reduction in $\|X\|_2$ is on the order of $$\left(1+\frac{4\pi^2}{n^2}\right)^{-n/2} \approx e^{-2\pi^2/n},$$ where we would like this to be $1$. For $n=100$, this is $0.82$, which is indeed drastic, but shouldn't be surprising in light of the mathematics.
If you want correctness to $d$ digits after one period, the condition on $n$ is that $e^{-2\pi^2/n}\sim 1-10^{-d}$, or $n\sim 2\pi^210^d$. So to get even three digits of accuracy you need $20,000$ steps or so.