I am trying to control that the following differentialequation is stable for the numerical time-integration method: Euler Forward with $\Delta t = 1$.
$$
y''(t)+y'(t)+\sin(y(t)) = 0, \ \ \ \ y(0 ) = 1, y'(1) = 0.
$$
By using the following, we can transform this into a system of differential equations. Since we know that our ODE is non-linear, we will use the Jacobi-matrix of $\bf{y'}= \begin{bmatrix} y_1'(t) \\ y_2'(t) \end{bmatrix}$.
\begin{align}
y_1&=y &\Rightarrow y_1' &= y_2,\\
y_2&=y' &\Rightarrow y_2'& = -y_2-\sin(y_1).
\end{align}
The Jacobi-matrix is:
\begin{equation}
J(y_1,y_2)= \begin{bmatrix} 0 & 1\\
-\cos(y_1) & -1 \end{bmatrix}.
\end{equation}
We will calculate the eigen-values of this matrix, these are:
\begin{align}
\lambda_1 &= \frac{-1 + \sqrt{1-4 \cos(y_1)} }{2},\\
\lambda_2 &= \frac{-1 – \sqrt{1-4 \cos(y_1)} }{2}.
\end{align}
Now, my question is: These eigenvalues can be complex, so how is it possible for Euler Forward to be stable, I recall that Euler Forward is not stable if the eigenvalues have complex parts.
Thanks for your time,
K. Kamal
Best Answer
The condition for stability is that for eigenvalues $λ$ with negative real part you have $$|1+hλ|<1.$$ For your problem these eigenvalue are $λ=\frac12(-1\pm i\sqrt{3})$ and $λ=-\frac12(1+\sqrt5)$ so that the condition is $$ (1-\tfrac12h)^2+\tfrac34h^2<1\iff h^2-h<0~~\text{ and }~~h(1+\sqrt5)<4 $$ and as $h>0$, the bound for stability is $h<\min(1, \frac{4}{1+\sqrt5})=1$.
However this only means that you avoid visibly, blatantly wrong behavior of the numerical solution, that is, that if the exact solution moves towards a fixed point, that the numerical solution were to explode away from it. The overall global error of size $O(h)$ or for variable time, $O((e^{Lt}-1)h)$, still remains. It may be advisable to take $h$ much smaller than the stability bound to get a sensible global error.
While the exact solution converges very rapidly towards the origin, with a bit more than one visible rotation about zero, the numerical solution with step size $h=0.5$ still somewhat shows that rapidity while $h=0.75$ has a rather slow convergence and $h=1$ looks undecided if it approaches the origin at all.