Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?


share|improve this question
Dissipation is to be expected for long times. Try with or some other symplectic integrator. – Juan M. Bello-Rivas 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. – Rhino Sep 4 '14 at 13:36
up vote 12 down vote accepted

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.

Conclusion: this is a serious issue with the method, so use a symplectic integrator that solves it.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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