Solution (source code)

= Solution

\b[A-stability removes the scalar linear decay restriction on the step size; it does not replace accuracy, nonlinear solvability or a suitable system norm.] Its importance is clearest for a <stiff differential equation>, where rapidly decaying modes coexist with slowly varying quantities of interest. On the <Dahlquist test equation> $y'=\lambda y$, the exact solution decays for $\operatorname{Re}\lambda<0$. A stable time integrator should not turn that decay into numerical growth.

For a one-step method write $y_{n+1}=R(z)y_n$, where $z=h\lambda$ and $R$ is its <stability function>. Its <linear stability domain> consists of parameters with $|R(z)|\le1$ and a well-defined update. The method is <A-stable> when this domain contains the closed left half-plane. A consistent order-$p$ one-step method has $R(z)=e^z+O(z^{p+1})$ near zero, but that local agreement alone says nothing about stability at large negative $z$.

For a <Runge-Kutta method> with stage <matrix> $A$, weights $b$ and $\mathbf1=(1,\ldots,1)^T$, eliminating the stages on the test equation gives
$$
R(z)=1+zb^T(I-zA)^{-1}\mathbf1.
$$
An explicit finite-stage method has a <polynomial> $R$. A nonconstant <polynomial> is unbounded along the negative real axis, so \b[no consistent explicit <Runge-Kutta method> is A-stable]. The <explicit Euler method>, with $R(z)=1+z$, is stable in the disk $|1+z|\le1$ and on the negative real axis only for $-2\le z\le0$. A stiff mode with a large negative <eigenvalue> therefore forces an unnecessarily small step even after its exact contribution has almost disappeared.

The <Backward Euler method> has $R(z)=1/(1-z)$. For $\operatorname{Re}z\le0$, $|1-z|^2=1-2\operatorname{Re}z+|z|^2\ge1$, so it is <A-stable>. The <trapezoidal rule> has
$$
R(z)=\frac{1+z/2}{1-z/2},\qquad
|1-z/2|^2-|1+z/2|^2=-2\operatorname{Re}z\ge0,
$$
and is also <A-stable>. The <implicit midpoint method> has the same scalar <stability function>, although its nonlinear update differs from the trapezoidal rule. The first method has order one and the latter two order two. Their implicit solves trade greater work per step for freedom from the scalar stiff stability restriction.

\Image[/past-exam-of-the-mathematics-course-of-the-university-of-cambridge/2010/iii/paper-63-stability-domains.png]
{title=Absolute stability domains of explicit Euler, backward Euler and the trapezoidal rule; shaded regions satisfy modulus of the amplification factor at most one}

Damping of very stiff modes separates <A-stability> from <L-stability>. For a one-step method, the latter additionally requires $R(z)\to0$ as $z$ tends to infinity in the left half-plane. <Backward Euler method> satisfies this. The trapezoidal and midpoint rules instead have $R(z)\to-1$ along the negative real axis: a very stiff mode can alternate in sign with little damping even though its amplitude remains bounded. <L-stability> is often preferable for diffusion and rapid relaxation, while retaining oscillatory modes can be desirable in conservative dynamics.

For a <linear multistep method>, write
$$
\sum_{j=0}^k\alpha_jy_{n+j}=h\sum_{j=0}^k\beta_jf_{n+j},\qquad
\rho(\zeta)=\sum_j\alpha_j\zeta^j,\quad \sigma(\zeta)=\sum_j\beta_j\zeta^j.
$$
Its test-equation amplification roots satisfy $\rho(\zeta)-z\sigma(\zeta)=0$. <Absolute stability> requires every root to lie in the <unit disk> or on its boundary with all unit-modulus roots simple; all parasitic modes matter. At $z=0$ this is the <zero-stability> root condition. Consistency requires $\rho(1)=0$ and $\rho'(1)=\sigma(1)$, and the <Dahlquist equivalence theorem> combines consistency and <zero-stability> to characterize convergence. Thus a claim of stability obtained by canceling a common factor is not a proof about the original recurrence and its arbitrary consistent starting data.

Explicit <linear multistep methods> cannot be <A-stable> either. If $\beta_k=0$ and $\alpha_k\ne0$, some lower coefficient of $\rho-z\sigma$ becomes unbounded as $z\to-\infty$, while the leading coefficient stays fixed. If all roots remained in the <unit disk>, every coefficient-to-leading-coefficient ratio, an elementary symmetric function of those roots, would remain bounded. This contradiction forces an unstable root.

An important implicit example is <BDF2 method>, whose characteristic equation is $(3-2z)\zeta^2-4\zeta+1=0$. A unit root $\zeta=e^{i\theta}$ gives
$$
z(\theta)=\tfrac12(3-4e^{-i\theta}+e^{-2i\theta}),\qquad
\operatorname{Re}z(\theta)=(1-\cos\theta)^2\ge0.
$$
The leading coefficient vanishes only at $z=3/2$, outside the left half-plane. Near a small negative $z$, both roots are inside the disk; no unit-circle crossing can occur in the connected open left half-plane. On the imaginary axis the boundary equality can occur only at $z=0$, where the roots are $1,1/3$. Hence BDF2 is <A-stable> and its two roots tend to zero for very large negative $z$, giving strong stiff-mode damping. This is a complete boundary-locus argument, not just a boundary plot.

The <Second Dahlquist barrier> says that a convergent <A-stable> ordinary <linear multistep method> has order at most two. This theorem does not limit <implicit Runge-Kutta methods> or <multiderivative multistep methods>. In particular the <A-stable third-order two-step multiderivative method> in question 1 uses $y''$ and satisfies different order conditions. Among <collocation Runge-Kutta methods>, the standard family results give <Gauss collocation methods> order $2s$ with diagonal $[s/s]$ <Padé approximants> to $e^z$ as <stability functions>, and <Radau IIA methods> order $2s-1$ with subdiagonal $[s-1/s]$ approximants. Both families are <A-stable>; the Gauss functions have a nonzero limit at infinity, whereas Radau IIA is <L-stable>. The corrected two-stage right-Radau method, for example, has
$$
R(z)=\frac{1+z/3}{1-2z/3+z^2/6},
$$
which matches the exponential through cubic order and tends to zero at infinity. Its poles $2\pm i\sqrt2$ lie in the right half-plane. For $z=-r+iy$, $r\ge0$, direct subtraction gives
$$
|1-2z/3+z^2/6|^2-|1+z/3|^2
=\frac{r(r+6)(r^2+2r+12)+2r(r+4)y^2+y^4}{36}\ge0,
$$
proving its <A-stability> directly.

Sectorial <A-alpha stability> relaxes the requirement to a sector around the negative real axis, useful when the spectrum of a stiff problem remains in that sector; higher-order backward differentiation formulas illustrate the usefulness of this weaker requirement. For a normal linear system $y'=Ly$, unitary diagonalization reduces the numerical update to its scalar factors $R(h\lambda_j)$. For a nonnormal system, a scalar <eigenvalue> test alone does not give a uniform bound in the chosen norm. Nonlinear stability is different again: <B-stability> asks for contractivity when $\operatorname{Re}\langle f(u)-f(v),u-v\rangle\le0$, and <algebraic stability of a Runge-Kutta method> is a standard sufficient criterion. Thus <A-stability> addresses a precise and valuable linear test, while step selection still has to meet accuracy requirements and the actual implicit equations must be solved on the relevant branch.