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$.
Best Answer
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.