I have the following two step Nyström method: $u_{k+1}=u_{k-1}+2h\cdot f(u_k)$. I want to know the stability region, so I wrote this as $w_{k+1} = A\cdot w_k$ with $A=\left(\begin{matrix} 2h\lambda & 1 \\ 1 & 0\end{matrix}\right)$ and $w_k = (u_{k},u_{k-1})^t$. I use the sample/test equation $u'=\lambda\cdot u = f(u)$ with $\lambda\in\mathbb{C}$. To find the stability region I must set some limitations on the matrix A.

Do I want the whole vector $w_{k+1}$ to go to the zero vector as $k$ goes to infinity, or just the first component that contains $u_{k+1}$? I don't think I can let the second component go to zero, since the matrix just has a $1$ in the bottomleft corner. So I think I want $2h\lambda$ to be negative in the reals (like Forward Euler) , or at least less than 1? I'm unfortunately just guessing at this point.

For (Forward) Euler we have $u_{k+1}=u_k + h\cdot f(u_k) = (1 + h\lambda)u_k$ and so we require that $h\lambda$ does not have a positive real part.

So I guess I also want some anologue for matrices being less than 1? I do not know how to deal with the matrices since the left hand vector not only contains $u_{k+1}$ but also the old $u_k$.


The eigenvalues of $A$ are $hλ\pm\sqrt{1+(hλ)^2}\approx \pm\exp(\pm hλ)$. The method is stable if for $Re(λ)<0$ the iteration converges to zero. As the product of the eigenvalues is always $-1$, you will always have one eigenvalue greater than $1$ in absolute value, and thus the component of the solution corresponding to it growing, not falling to $0$.

In conclusion, the method is nowhere stable.

Following Hairer/Wanner: "Solving ODE II", ch. V.1 "Stability of multi-step methods", the characteristic polynomials here are $\rho(\zeta)=\zeta^2-1$ and $\sigma(\zeta)=2\zeta$ with characteristic equation $\rho(\zeta)-\mu\sigma(\zeta)=0$, where $\mu=hλ$, and $\zeta=e^{hλ}$ is the step factor of the exact solution.

Then the stability domain is defined as the set of all $\mu$ where the characteristic equation only has solutions in the closed unit disk. But again, here the solutions have product $-1$, so they can not all be inside the unit disk, and both roots are on the boundary for $\mu$ on the segment between $-i$ and $i$. This is called "weak instability" and causes the oscillating behavior (slowly increasing from floating point noise) shown in my previous answer resp. vol. I, chapter III.9 of the cited book.

  • $\begingroup$ Are you assuming $\lambda$ has to be a real number? It can be complex, I will update my question. $\endgroup$ Jan 15 '19 at 21:16
  • $\begingroup$ If I impose conditions on the real part of $λ$, it obviously can be complex. $\endgroup$ Jan 15 '19 at 21:21
  • $\begingroup$ Okay, so you're saying one component of the solution isn't falling to 0, but the second component can stay at 1, because that would just mean that the second component equals 1 times the same term because of the 1 in the bottom left? $\endgroup$ Jan 15 '19 at 21:40
  • 1
    $\begingroup$ For $h$ small enough, $q=hλ+\sqrt{1+(hλ)^2}=hλ+(1+\frac12(hλ)^2-\frac18(hλ)^4+...)=e^{hλ}+O((hλ)^3)$. For stability you need $\|A^k\|\le 1$ for large $k$. The eigenvalues of $A$ tell you if that is the case. $\endgroup$ Jan 15 '19 at 22:39
  • 1
    $\begingroup$ Yes it can, for the differential equation $y'=0$, that is, for $λ=0$. Hopefully any implementation of a solver will return the constant solution in this case. Also for $y''=-y$ you get $λ=\pm i$ so that for $h<1$ the points $\pm ih\pm\sqrt{1-h^2}$ are on the unit circle. But that is the best you get. $\endgroup$ Jan 16 '19 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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