# Stability region for two step Nyström method

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.

• Are you assuming $\lambda$ has to be a real number? It can be complex, I will update my question. Jan 15 '19 at 21:16
• If I impose conditions on the real part of $λ$, it obviously can be complex. Jan 15 '19 at 21:21
• 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? Jan 15 '19 at 21:40
• 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. Jan 15 '19 at 22:39
• 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. Jan 16 '19 at 21:51