2

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

ArianJ
  • 690

1 Answers1

1

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.

numerical solution for 3 step sizes

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.

Lutz Lehmann
  • 131,652
  • How do you get those eigenvalues? I don't understand how you 'lost' the $\cos( y_1 )$. – ArianJ Jun 27 '18 at 19:24
  • I used the extremes of $\cos(y)\in [-1,1]$. Actually relevant is only the behavior at $(y,y')\approx(0,0)$ where $\cos(y)\approx 1$. Stability at the saddle points is somewhat less relevant. – Lutz Lehmann Jun 27 '18 at 19:33
  • But then again, we get complex eigenvalues; so therefore Euler Forward should not give stability for any value of $h$ right? – ArianJ Jun 27 '18 at 20:13
  • Please give your definition of stability which excludes complex eigenvalues. As per the arguments above, $h=Δt=1$ in not stable in my definition which is the one of Dahlquist A-stability. Which is the one where you take $y'=λy$ so that the exact solution converges to $0$ and demand that also the numerical solution for the given method with stepsize $h$ converges to $0$, – Lutz Lehmann Jun 27 '18 at 21:22
  • I see what my problem was: I read somewhere that when your system has \bf{pure imaginary} eigenvalues Euler Forward can't be used (we can see this in the stability region http://www.s3datascience.com/_/rsrc/1481650194998/z-s3-data/educating-the-user/a-brief-history-of-calculation/StabilityRegion.png). And here is where my mistake comes in: You can have complex eigenvalues and Euler Forward can give a stable method according to Dahlquist A-stability if $| 1+h \cdot \lambda| < 1$. Thanks for your explanation. – ArianJ Jun 28 '18 at 08:32