= Solution
An $s$-stage <Runge-Kutta method> for $y'=f(t,y)$ is defined by
$$
Y_i=y_n+h\sum_{j=1}^sa_{ij}f(t_n+c_jh,Y_j),\qquad
y_{n+1}=y_n+h\sum_{i=1}^sb_if(t_n+c_ih,Y_i).
$$
The <Butcher tableau> records $A=(a_{ij})$, $b$ and $c$; internal consistency usually sets $c=Ae$, and first-order consistency requires $b^Te=1$. A strictly lower triangular $A$ gives an explicit method, while other stage dependencies generally require an <implicit Runge-Kutta method>. The <stage solvability of an implicit Runge-Kutta method> must be checked separately: if $f$ is globally <Lipschitz continuous> in $y$ with constant $L$, the stage fixed-point map is a contraction in the maximum stage norm whenever $hL\max_i\sum_j|a_{ij}|<1$. This is a sufficient small-step condition, not a restriction intrinsic to the definitions of <A-stability> or <B-stability>.
Linear stability begins with the <Dahlquist test equation> $y'=\lambda y$. Solving the stages gives
$$
\boxed{y_{n+1}=R(z)y_n,\qquad
R(z)=1+zb^T(I-zA)^{-1}e,\qquad z=h\lambda.}
$$
The <linear stability domain> consists of test values for which the stage system is well defined and $|R(z)|\leq1$. <A-stability> means this includes the closed left half-plane. For $\operatorname{Re}\lambda<0$, this gives stability with no scalar decay-mode step restriction. For a <normal matrix> $L$ in $y'=Ly$, <unitary diagonalization of a normal matrix> reduces the <norm> estimate to these scalar factors. For a <non-normal matrix>, eigenvalues alone do not establish a uniform <norm> bound: eigenvector conditioning and transient amplification matter. Mesh-uniform estimates for discretized PDEs must control operator powers, not only their spectra.
Every consistent explicit <Runge-Kutta method> has a nonconstant polynomial <stability function> and therefore cannot be <A-stable>, since that polynomial is unbounded on the negative real axis. For example, <forward Euler method> has $R(z)=1+z$, stability disk $|1+z|\leq1$, and restriction $0\leq hq\leq2$ for $\lambda=-q<0$. For a pure imaginary mode $\lambda=i\omega\neq0$, $|1+ih\omega|>1$, so it is unstable for every positive step. The <classical fourth-order Runge-Kutta method> instead has
$$
R(z)=1+z+\tfrac12z^2+\tfrac16z^3+\tfrac1{24}z^4,
\qquad |R(iq)|^2=1-\frac{q^6}{72}+\frac{q^8}{576}.
$$
It is stable on the imaginary axis exactly for $|q|\leq2\sqrt2$, an illustrative conditional stability range for oscillatory evolution.
<Backward Euler method> has $R(z)=(1-z)^{-1}$ and is <A-stable>. <L-stability> adds $R(z)\to0$ as $|z|\to\infty$ in the left half-plane; backward Euler satisfies this and strongly damps unresolved rapidly decaying modes in a <stiff differential equation>. The <implicit midpoint rule> and the <trapezoidal rule> both have $R(z)=(1+z/2)/(1-z/2)$, so they are <A-stable> but not <L-stable>. They preserve the modulus of a pure imaginary test mode, but their stiff-decay limit is $-1$, leaving oscillatory numerical remnants. The three-stage <Lobatto IIIA method> has the fourth-order rational function found in Question 2 and likewise lacks <L-stability>. Thus high order and <A-stability> do not themselves imply efficient stiff damping.
Nonlinear stability measures differences between solutions of the same equation. A <one-sided Lipschitz condition> is
$$
\operatorname{Re}\langle f(t,y)-f(t,\widetilde y),y-\widetilde y\rangle
\leq\nu\|y-\widetilde y\|^2.
$$
Differentiating the squared difference and applying the <Gronwall inequality> gives $\|y(t)-\widetilde y(t)\|\leq e^{\nu(t-t_0)}\|y(t_0)-\widetilde y(t_0)\|$. A <dissipative vector field> has $\nu\leq0$ and is contractive. A <B-stable> method preserves this property for every positive step size for which the stages are well defined:
$$
\boxed{\|y_{n+1}-\widetilde y_{n+1}\|\leq\|y_n-\widetilde y_n\|.}
$$
This includes time-dependent vector fields when dissipativity holds at each common time argument. Applying it to scalar linear dissipative fields shows that <B-stability> implies <A-stability> when the linear stages are well defined; the converse fails.
The standard sufficient theorem is \b[<algebraic stability> implies <B-stability>], subject to stage solvability. Define $B=\operatorname{diag}(b)$ and $M=BA+A^TB-bb^T$. <Algebraic stability of a Runge-Kutta method> means $b_i\geq0$ and $M$ is a <positive semidefinite matrix>. To prove the theorem, set $d=y_n-\widetilde y_n$, $D_i=Y_i-\widetilde Y_i$ and $F_i=f(t_n+c_ih,Y_i)-f(t_n+c_ih,\widetilde Y_i)$. Expanding the output difference and using $D_i=d+h\sum_ja_{ij}F_j$ yields the <Runge-Kutta contractivity identity>
$$
\|d_{n+1}\|^2=\|d\|^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.
$$
Dissipativity makes the first sum nonpositive, while positive semidefiniteness makes the last <quadratic form> nonnegative. The latter follows by factoring $M=Q^TQ$ and writing it as a sum of squared <norms> of linear combinations of the $F_i$. This proves the claimed contraction. It is a sufficient theorem; absence of <algebraic stability> is not, by itself, a proof of failure of <B-stability> for every representation.
Several examples make the distinctions concrete. For backward Euler, $b=1$, $A=(1)$, so $M=(1)$: it is both <L-stable> and <B-stable>. For implicit midpoint, $b=1$, $A=(1/2)$, so $M=(0)$: it is <B-stable> and preserves squared <norms> for a skew-Hermitian linear equation, but lacks stiff damping. For the two-stage <Gauss--Legendre Runge-Kutta method>,
$$
A=\begin{pmatrix}1/4&1/4-\sqrt3/6\\1/4+\sqrt3/6&1/4\end{pmatrix},
\qquad b=(1/2,1/2)^T,
$$
the matrix $M$ vanishes as well. The collocation theorem gives order four, and its <stability function> is exactly the same as the three-stage Lobatto IIIA function. Yet the Gauss method is <algebraically stable>, whereas that Lobatto method is not. Identical scalar linear stability functions need not imply identical nonlinear stability properties.
An example combining stiff damping and nonlinear contraction is the two-stage <Radau IIA method>, with
$$
A=\begin{pmatrix}5/12&-1/12\\3/4&1/4\end{pmatrix},\quad
b=(3/4,1/4)^T,\quad
M=\frac1{16}\begin{pmatrix}1&-1\\-1&1\end{pmatrix},\quad
R(z)=\frac{1+z/3}{1-2z/3+z^2/6}.
$$
Its positive weights and positive semidefinite $M$ prove <B-stability>; its rational function has denominator roots $2\pm i\sqrt2$ and, for $z=x+iy$, satisfies
$$
|1-2z/3+z^2/6|^2-|1+z/3|^2
=\frac{|z|^4-8x|z|^2+24x^2-72x}{36}\geq0\quad(x\leq0).
$$
Thus it is <A-stable>; the limit $R(z)\to0$ then proves <L-stability>. The standard Radau IIA collocation theorem gives order $2s-1$, here three. These properties explain its usefulness for stiff nonlinear equations.
Finally, \b[<A-stability> alone does not imply <B-stability>]. The <trapezoidal rule fails B-stability> even for the scalar dissipative equation $y'=-y^3$. With $h=1$, its step map $T$ is uniquely defined by
$$
T(y)+\tfrac12T(y)^3=y-\tfrac12y^3.
$$
At $y=\sqrt2$, $T(y)=0$. Implicit differentiation gives
$$
T'(y)=\frac{1-\tfrac32y^2}{1+\tfrac32T(y)^2},\qquad T'(\sqrt2)=-2.
$$
Thus nearby starting values are separated by approximately twice their original distance after one step, although the exact flow is contractive. This example and the contractivity theorem show why scalar linear analysis, stiff damping and genuinely nonlinear <norm> estimates answer different stability questions.
Back to article page
