1

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 Coding Wombat
  • 657
  • 3
  • 20

1 Answers1

2

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.

Lutz Lehmann
  • 126,666
  • Are you assuming $\lambda$ has to be a real number? It can be complex, I will update my question. – The Coding Wombat Jan 15 '19 at 21:16
  • If I impose conditions on the real part of $λ$, it obviously can be complex. – Lutz Lehmann 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? – The Coding Wombat Jan 15 '19 at 21:40
  • No, it means that $u_k=aq^k+b(-q)^{-k}$ where $q=hλ+\sqrt{1+(hλ)^2}$ is the root closer to $1$. Assume $Re(λ)<0$ so that $|q|<1$, at least for small $h$. You can construct $u_1$ in such a way that $b=0$ initially. But during the computation floating point errors accumulate. Their contribution is distributed to both $a$ and $b$. With time, the first term will fall to zero, but in the second term $|q|^{-k}$ will be large enough that this term becomes noticeable, even if $b$ stays a small multiple of the machine number. – Lutz Lehmann Jan 15 '19 at 21:50
  • How do the eigenvalues of $A$ say anything about the stability? I know that if $\tau$ is an eigenvalue of $A$ we have that $A\cdot w_k=\tau\cdot w_k$. And where can I find more on how you get to the approximation $e^{h\lambda}$? – The Coding Wombat Jan 15 '19 at 22:29
  • 1
    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. – Lutz Lehmann Jan 15 '19 at 22:39
  • How is that norm defined? – The Coding Wombat Jan 16 '19 at 21:06
  • This is true for any operator norm, $\limsup \sqrt[n,]{|A^n|}\le 1 $ if and only if all eigenvalues are $\le 1$ in abs. value and when $=1$ they have to be simple. – Lutz Lehmann Jan 16 '19 at 21:15
  • But if the product of the eigenvalues is -1, can't one eigenvalue be 1 and the other -1? So that neither of them grow? Wouldn't that suffice? – The Coding Wombat Jan 16 '19 at 21:38
  • 1
    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. – Lutz Lehmann Jan 16 '19 at 21:51