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!
Non-conservative implementation implicit Euler
ode;time integration;implicit methods
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.