Stability regions and stiff equations

Let $d\ge 1$ be an integer and let $G\in {\mathbb{R}}^{d\times d}$ be symmetric positive definite. Define the inner product and its vector norm on ${\mathbb{R}}^d$ by

$⟨x,y{⟩}_G=x^{\mathsf{T}}Gy,\mathrm{\parallel }x{\mathrm{\parallel }}_G=\sqrt{x^{\mathsf{T}}Gx}\mathrm{.}$

For $M\in {\mathbb{R}}^{d\times d}$, its induced operator norm is $\mathrm{\parallel }M{\mathrm{\parallel }}_G={\sup }_{x\mathrm{\ne }0}\mathrm{\parallel }Mx{\mathrm{\parallel }}_G\mathrm{/}\mathrm{\parallel }x{\mathrm{\parallel }}_G$. A matrix step $M$ is contractive in this specified inner product when $\mathrm{\parallel }M{\mathrm{\parallel }}_G\le 1$.

For $A\in {\mathbb{R}}^{d\times d}$, its logarithmic norm in the same inner product is

${\mu }_G(A)={\sup }_{x\mathrm{\ne }0}\frac{x^{\mathsf{T}}GAx}{x^{\mathsf{T}}Gx}\mathrm{.}$

In the Euclidean inner product $G=I$, write $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_2$ and ${\mu }_2$. A real matrix is normal when $A{A}^{\mathsf{T}}=A^{\mathsf{T}}A$ and nonnormal otherwise.

For a rational stability function $R=P\mathrm{/}Q$, the matrix step $R(hA)$ is defined only when $Q(hA)$ is invertible; its contractivity is the separate condition $\mathrm{\parallel }R(hA){\mathrm{\parallel }}_G\le 1$.

For the scalar test equation $y^{\mathrm{\prime }}=\lambda y$, a step size $h>0$ gives $z=h\lambda$. A one-step scalar update has stability function $R$ when one step sends $y_n$ to $y_{n+1}=R(z)y_n$. If $R$ is rational, write it in reduced form; it is defined at $z$ exactly when its denominator is nonzero there. Its absolute-stability region is

${\mathcal{S}}_R=\{z\in \mathbb{C}:R(z)\text{ is defined and }\mathrm{\mid }R(z)\mathrm{\mid }\le 1\}\mathrm{.}$

The update is A-stable when the closed left half-plane $\{z:\mathrm{Re}z\le 0\}$ is contained in ${\mathcal{S}}_R$. A rational one-step update is L-stable when it is A-stable and

${\lim }_{\begin{array}{c}\mathrm{\mid }z\mathrm{\mid }\to \mathrm{\infty }\\ \mathrm{Re}z\le 0\end{array}}R(z)=0,$

where the limit means that the values tend to zero along every sequence in the closed left half-plane whose modulus tends to infinity.

For a linear multistep method with the coefficient and polynomial convention of Linear multistep methods, consistency, order, and the root condition, substitution of the scalar test equation gives the characteristic polynomial

${\pi }_z(\zeta )=\rho (\zeta )-z\sigma (\zeta )\mathrm{.}$

The multistep absolute-stability region consists of those $z$ for which the leading coefficient ${\alpha }_k-z{\beta }_k$ is nonzero and every root of ${\pi }_z$ satisfies the root condition: all roots have modulus at most one, and every root of modulus one is simple. This criterion concerns constant-step scalar test recurrences.

Absolute stability is a scalar test-equation property. It does not itself assert an accuracy bound for an ODE solution, contractivity of a matrix step, positivity preservation, or conservation.

Backward Euler on $y^{\mathrm{\prime }}=\lambda y$ is the scalar update $y_{n+1}=y_n+z{y}_{n+1}$, where $z=h\lambda$. Its stability function is

$R_{\mathrm{B}\mathrm{E}}(z)=\frac{1}{1-z}\mathrm{.}$

It is A-stable and L-stable.

The trapezoidal rule on $y^{\mathrm{\prime }}=\lambda y$ is

$y_{n+1}=y_n+\frac{z}{2}(y_n+y_{n+1}),z=h\lambda \mathrm{.}$

Its stability function is $R_{\mathrm{T}}(z)=(1+z\mathrm{/}2)\mathrm{/}(1-z\mathrm{/}2)$. It is A-stable but not L-stable.

On a uniform grid, the variable-step BDF2 algorithm of Variable-step BDF2 with an explicit local implicit-branch rule has the scalar test-equation polynomial

$P_z(\zeta )=(3-2z){\zeta }^2-4\zeta +1,z=h\lambda \mathrm{.}$

If $\zeta$ is a root and $\mathrm{\mid }\zeta \mathrm{\mid }\ge 1$, then $\mathrm{Re}z\ge 0$. Equality is possible only for $z=0$ and $\zeta =1$.

For forward Euler applied to $y^{\mathrm{\prime }}=\lambda y$, let $z=h\lambda$. Its stability function is $R_{\mathrm{E}}(z)=1+z$, and its absolute-stability region is exactly the closed disk

$\{z\in \mathbb{C}:\mathrm{\mid }z+1\mathrm{\mid }\le 1\}\mathrm{.}$

On the negative real axis, for $y^{\mathrm{\prime }}=-\rho y$ with $\rho >0$, the method is absolutely stable exactly when $0<h\rho \le 2$.

For each $z=h\lambda$ with $\mathrm{Re}z\le 0$, the constant-step BDF2 test polynomial $P_z(\zeta )=(3-2z){\zeta }^2-4\zeta +1$ has degree two, all its roots have modulus at most one, and every root of modulus one is simple. Thus its multistep absolute-stability region contains the closed left half-plane, so BDF2 is A-stable.

For $\mathrm{Re}z\le 0$, let $\zeta$ be any root of the constant-step BDF2 test polynomial $P_z(\zeta )=(3-2z){\zeta }^2-4\zeta +1$. As $\mathrm{\mid }z\mathrm{\mid }\to \mathrm{\infty }$ within the closed left half-plane, every such root satisfies $\zeta \to 0$.

This is a multistep damping statement. It does not extend the one-step definition of L-stability to multistep methods.