Solution (source code)

= Solution

An $s$-stage <Runge-Kutta method> has stages and update
$$
Y_i=y_n+h\sum_{j=1}^sa_{ij}f(Y_j),
\qquad
y_{n+1}=y_n+h\sum_{i=1}^sb_if(Y_i).
$$
On the <Dahlquist test equation> $y'=\lambda y$, elimination of the stages gives the <stability function>
$$
R(z)=1+z b^T(I-zA)^{-1}\mathbf1,
\qquad z=h\lambda.
$$
The <linear stability domain> is the set where $|R(z)|\leq1$. A method is <A-stable> when this domain contains $\operatorname{Re}z\leq0$, so every exactly decaying scalar linear mode remains bounded for every step size. It is <L-stable> when it is A-stable and $R(z)\to0$ as $|z|\to\infty$ in the left half-plane; this extra limit strongly damps unresolved stiff modes.

The rational function $R$ makes several useful conclusions immediate. No explicit Runge--Kutta method is A-stable because its stability function is a nonconstant <polynomial> and is therefore unbounded on the negative real axis. The <implicit midpoint rule> has $R(z)=(1+z/2)/(1-z/2)$ and is A-stable, but $R(z)\to-1$, so it is not L-stable. The <Backward Euler method> has $R(z)=(1-z)^{-1}$ and is L-stable. More generally, a rational $R$ with no pole in the closed left half-plane is A-stable if and only if $|R(iy)|\leq1$ for every real $y$; this follows by applying the <maximum modulus principle> on expanding left half-disks.

Scalar linear stability does not by itself control nonlinear perturbations. Suppose the vector field is dissipative in the sense that
$$
\operatorname{Re}\langle f(u)-f(v),u-v\rangle\leq0.
$$
A method is <B-stable> if it preserves the resulting contractivity: two numerical solutions satisfy $\|y_{n+1}-\widetilde y_{n+1}\|\leq\|y_n-\widetilde y_n\|$. A practical sufficient condition is <algebraic stability of a Runge-Kutta method>: $b_i\geq0$ and
$$
M=BA+A^TB-bb^T\succeq0,
\qquad B=\operatorname{diag}(b_i).
$$
To prove the implication, let $D_i=Y_i-\widetilde Y_i$ and $F_i=f(Y_i)-f(\widetilde Y_i)$. Expanding the squared distance and substituting the stage equations gives the <Runge-Kutta contractivity identity>
$$
\|y_{n+1}-\widetilde y_{n+1}\|^2
=\|y_n-\widetilde y_n\|^2
+2h\sum_i b_i\operatorname{Re}\langle D_i,F_i\rangle
-h^2\sum_{i,j}m_{ij}\operatorname{Re}\langle F_i,F_j\rangle.
$$
The dissipativity inequalities make the middle sum nonpositive, and positive semidefiniteness of $M$ makes the final quadratic form nonnegative before its minus sign. The distance therefore cannot increase. In particular, algebraic stability implies B-stability and, by applying contractivity to the real two-dimensional form of $y'=\lambda y$, implies A-stability.

Important collocation families illustrate these notions. <Gauss--Legendre Runge-Kutta method>[Gauss methods] are A-stable, symmetric, and have order $2s$, but they do not damp infinitely stiff modes. <Radau IIA method>[Radau IIA methods] have order $2s-1$, are algebraically stable, and are L-stable. These properties explain why A-stability controls unrestricted linear decay, L-stability is useful for stiff transients, and algebraic or B-stability is the stronger tool for nonlinear dissipative equations.