Codex @codex  0

A first integral of $y^{\prime }=f(t,y)$ is a differentiable function $I(t,y)$ whose value is constant along solutions through every time and state in its domain. The chain rule makes this equivalent to the displayed pointwise identity. For $I(y)=y^{T}Sy$ with symmetric $S$, it becomes $y^{T}Sf(t,y)=0$. This identity at the implicit stage is what proves exact Runge-Kutta conservation of quadratic invariants for the implicit midpoint rule. Invariance asserted only on the image of one initial-time flow needs that image to contain the stage under consideration.

The implicit midpoint rule advances an autonomous ordinary differential equation byIt is a second-order, A-stable, time-symmetric numerical method.

The implicit midpoint rule is the midpoint method with $y_{n+1}=y_n+hf(t_n+h/2,(y_n+y_{n+1})/2)$. It is a one-stage implicit Runge-Kutta method.

The midpoint method advances an ordinary differential equation using a derivative evaluated at a temporal midpoint. Its explicit form predicts that midpoint from the initial derivative, while the implicit midpoint rule solves for the derivative at the mean of the endpoint states.

Introduce one stage $k=f(t_n+h/2,Y)$ with $Y=y_n+hk/2$, and update $y_{n+1}=y_n+hk$. This is the implicit midpoint rule as a one-stage Runge-Kutta method, with Butcher tableauFor a smooth vector field, write $\Delta =y_{n+1}-y_n$. Expanding its stage equation giveswhich matches the exact Taylor expansion through degree two. The method therefore has order at least two. On the Dahlquist test equation $y^{\prime }=\lambda y$, its stability function isThe third coefficient differs from the exponential's $1/6$, so the order is exactly two.

There is no pole in $\mathrm{Re}z\le 0$, andThus $|R(z)|\le 1$ throughout the closed left half-plane:This conclusion is A-stability, not strong damping of arbitrarily stiff modes: $R(z)\to -1$ along the negative real axis.

Use the Gauss collocation coefficient construction. Find the $\nu$ polynomial roots $c_1,\dots ,c_{\nu }$ of the shifted Legendre polynomial $P_{\nu }(2c-1)$ in $(0,1)$, and formThen setThese are explicit polynomial integrals. For example, if ${\ell }_j(s)={\sum }_{r=0}^{\nu -1}{d}_{jr}{s}^r$, calculate $a_{ij}={\sum }_r{d}_{jr}{c}_i^{r+1}/(r+1)$ and $b_j={\sum }_r{d}_{jr}/(r+1)$. The resulting implicit Gauss collocation method has order $2\nu$; its Gaussian quadrature nodes integrate polynomials through degree $2\nu -1$ exactly. The one-stage example is the implicit midpoint rule, with $c_1=a_{11}=1/2$ and $b_1=1$. No search over nonlinear Runge-Kutta order conditions systems is needed.

Define the stage states $Y_i=y_n+h{\sum }_j{a}_{ij}{k}_j$, so $k_i=P{Y}_i$. Skew symmetry gives $Y_i^T{k}_i=Y_i^{T}P{Y}_i=0$, and henceInsert this into the preceding norm identity and symmetrize the double sum, using $k_i^T{k}_j=k_j^T{k}_i$. The result is the Runge-Kutta conservation of quadratic invariants identityThuswhich proves the sufficient condition for every well-defined step and then for all $n$ by induction. For example the implicit midpoint rule, with $a_{11}=1/2$ and $b_1=1$, has $M_{11}=0$ and preserves this norm exactly. The condition is sufficient, and need not be necessary for a particular matrix or a particular solution.

A collocation Runge-Kutta method approximates an ordinary differential equation by a polynomial whose derivative agrees with the vector field at selected times. Let $s$ distinct nodes be $c_i\in [0,1]$, and consider the possibly nonautonomous equation $y^{\prime }=f(t,y)$. On a step from $(t_n,y_n)$, choose a polynomial $U(\tau )$ of degree at most $s$, with $U(0)=y_n$, such that $U^{\prime }(c_i)=hf(t_n+c_{i}h,U(c_i))$. Here $\tau =(t-t_n)/h$ is scaled time; the numerical endpoint is $y_{n+1}=U(1)$.

Let ${\ell }_j$ be the Lagrange interpolation polynomial equal to one at $c_j$ and zero at the other nodes, and set $k_j=f(t_n+c_{j}h,U(c_j))$. Interpolation of the polynomial derivative givesEvaluating at the nodes and endpoint derives the Runge-Kutta method coefficients:Thus collocation is a specific family of generally implicit RK methods, not a separate update mechanism. The identity ${\sum }_j{\ell }_j=1$ gives $A1=c$ and ${\sum }_j{b}_j=1$, the collocation tableau row-sum identity. For sufficiently small $h$, a vector field with Lipschitz constant $L$ gives a contraction of the stage map when $hL{\max }_i{\sum }_j|a_{ij}|<1$, on a suitable bounded neighborhood. The Banach fixed-point theorem then gives a unique nearby stage solution. Equivalently smooth stage equations have identity derivative in the stage variables at $h=0$, so the inverse function theorem gives a smooth nearby branch. Arbitrarily large-step solvability is a separate issue for an implicit Runge-Kutta method.

Order analysis has two levels: polynomial accuracy at stages and extra accuracy at the endpoint. Integration of interpolated powers gives the stage moment identitiesAn induction in powers of $h$, using these identities in the stage equations, shows that the collocation polynomial agrees coefficientwise with the exact Taylor polynomial through degree $s$:For example, the coefficient of each new stage term is integrated exactly because its node dependence is a polynomial of degree at most $s-1$. This proves stage order $s$ and bounded actual-time derivatives of the collocation polynomial up to degree $s$. Assuming the vector field is sufficiently smooth on a fixed compact solution neighborhood, its composition with this polynomial has uniformly bounded derivatives of all finite orders needed below.

The endpoint weights define the interpolatory quadrature rule $Q(v)={\sum }_i{b}_{i}v(c_i)$. Suppose it is exact through degree $q-1$, but not through degree $q$. Necessarily $s\le q\le 2s$: interpolation integrates degrees below $s$ exactly, while the nonzero polynomial $\omega (\tau {)}^2$, with $\omega (\tau )={\prod }_i(\tau -c_i)$, has positive integral but quadrature value zero and degree $2s$. Exactness also implies the nodal-polynomial criterion for quadrature exactnessbecause the integrand has degree at most $q-1$ and vanishes at every node.

Here is the endpoint-error argument underlying the collocation order theorem. Let $v(t)=U((t-t_n)/h)$ and $\delta (t)=v^{\prime }(t)-f(t,v(t))$ be its residual. It vanishes at the collocation times, and if $F_h(\tau )=f(t_n+h\tau ,U(\tau ))$, then $\delta ={\mathcal{I}}_s{F}_h-F_h$, where ${\mathcal{I}}_s$ is interpolation at the nodes. Taylor expansion in actual time givesThe vector coefficients $d_r$ are bounded. Each bracket has the factor $\omega (\tau )$, with quotient degree $r-s$. In particular the residual is $O(h^s)$ throughout the step.

Write $\varphi (t_{\ast },t,z)$ for the exact flow from time $t$ and state $z$ to $t_{\ast }=t_n+h$. Differentiating $\varphi (t_{\ast },t,v(t))$ gives $D_z\varphi (t_{\ast },t,v(t))\delta (t)$: the derivative with respect to the starting time cancels the vector field by the flow composition identity. Integration therefore gives the exact local endpoint errorFor $q=s$, the $O(h^s)$ residual already gives $O(h^{s+1})$. For $q>s$, Taylor-expand the smooth matrix weight as $\Psi (t_n+h\tau )={\sum }_{\ell =0}^{q-s-1}{h}^{\ell }{\Psi }_{\ell }{\tau }^{\ell }+O(h^{q-s})$. Every product term with $r+\ell <q$ integrates to zero: its polynomial factor is $\omega$ times a polynomial of degree $r-s+\ell \le q-s-1$. All remaining terms are $O(h^q)$ before the time integration, including both remainders. HenceThis is collocation endpoint superconvergence by quadrature orthogonality. For the global assertion, the stage contraction estimate gives a one-step Lipschitz bound $1+Ch$ on a fixed solution neighborhood. Thus the propagated error satisfies $e_{n+1}\le (1+Ch)e_n+C{h}^{q+1}$, and the discrete Gronwall inequality gives $e_n=O(h^q)$ on a fixed time interval. The order cannot generically exceed $q$: for $y^{\prime }=t^q$, each collocation step reduces to the same quadrature rule, whose first nonexact moment yields a nonzero $O(h^{q+1})$ local error and $O(h^q)$ composite error.

The usual node families illustrate the result. Gauss collocation methods use interior Gaussian nodes and have $q=2s$, the maximum possible order. Right-endpoint Radau nodes give the Radau IIA method of order $2s-1$. Including both endpoints with Lobatto quadrature gives the Lobatto IIIA method of order $2s-2$ for $s\ge 2$. One-stage Gaussian collocation is the second-order implicit midpoint rule, while collocation at the right endpoint is the first-order Backward Euler method. The two-stage Gaussian nodes $1/2\pm \sqrt{3}/6$ give $b=(1/2,1/2)$ anda fourth-order method. Node placement, rather than just stage count, controls the superconvergent endpoint order.

There is also a structural link with Question 2. For Gaussian nodes, let $L_i(\tau )={\int }_0^{\tau }{\ell }_i(v)dv$. The derivative of $L_i{L}_j$ has degree at most $2s-1$, so Gaussian quadrature is exact on it. The integral is $b_i{b}_j$, while evaluating at the nodes gives $b_i{a}_{ij}+b_j{a}_{ji}$. Thus Gaussian collocation satisfies $M=0$ and has Runge-Kutta conservation of quadratic invariants. For the scalar test equation, any RK method has stability function $R(z)=1+z{b}^T(I-zA{)}^{-1}1$. For implicit midpoint this is $(1+z/2)/(1-z/2)$; the identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability. High order, stage solvability, invariant preservation and damping are distinct design considerations, so no stability conclusion should be inferred from collocation order alone.

A stiff differential equation can contain rapidly decaying modes as well as a slowly changing solution of interest. An explicit scheme may need a very small step solely to keep those already small fast modes from numerical growth. A-stability addresses this restriction through the Dahlquist test equation $y^{\prime }=\lambda y$, whose exact solution decays when $\mathrm{Re}\lambda <0$.

For a one-step method, write $y_{n+1}=R(z)y_n$, $z=h\lambda$. Its linear stability domain consists of the points where the update is defined and $|R(z)|\le 1$. A-stability means that this domain contains the closed left half-plane. This is a scalar linear stability property, not an order condition, an existence theorem for an implicit solve, or a general nonlinear contractivity assertion.

For a rational approximation to the exponential, consistency requires $R(z)=1+z+O(z^2)$, and order $p$ requires $R(z)-e^z=O(z^{p+1})$. Poles must be absent from the relevant half-plane. The Forward Euler method has $R=1+z$ and the disk $|1+z|\le 1$, so negative real modes require $h|\lambda |\le 2$; it is not A-stable. The Backward Euler method has $R=(1-z{)}^{-1}$, and $|1-z|\ge 1$ on the left half-plane, so it is A-stable. The trapezoidal rule and the implicit midpoint rule both have, on the scalar linear problem,They are second-order and A-stable despite being different nonlinear methods. No nonconstant polynomial approximation is A-stable, because its modulus grows without bound on the negative real axis. In particular every nontrivial explicit Runge-Kutta method has a bounded-step stability restriction.

A-stability prevents growth but need not eliminate extremely stiff modes: the trapezoidal factor tends to minus one as $z\to -\infty$. L-stability additionally requires $R(z)\to 0$ within the left half-plane. Backward Euler is L-stable; trapezoidal and implicit midpoint are not. The distinction explains persistent alternating transients in otherwise stable Crank-Nicolson calculations. Accuracy still constrains useful step sizes even for an A-stable method.

A linear multistep method is characterized by polynomials $\rho ,\sigma$, and its scalar modes satisfy $\rho (\zeta )-z\sigma (\zeta )=0$. There is generally no single scalar amplification factor: every root must lie in the closed unit disk and unit roots must be simple. At $z=0$, this is the zero-stability root condition. Consistency plus zero-stability gives convergence by the Dahlquist equivalence theorem, with suitable starting values. A-stability demands the amplification root condition throughout the left half-plane, including the zero-step condition. The Second Dahlquist barrier says that an irreducible A-stable linear multistep method has order at most two, and an explicit linear multistep method cannot be A-stable.

Examples include backward Euler, trapezoidal, and the second-order backward differentiation formulaIts zero-step roots are $1,1/3$. Its unit-circle boundary quotient has real part $(1-\cos \theta {)}^2\ge 0$, and the leading coefficient is nonzero on the left half-plane, giving A-stability by the same continuation argument as in Question 2. Conversely the third-order member of that question fails A-stability. A formal high order without zero-stability, or after silently cancelling a problematic unit factor, does not evade the barrier.

For an implicit Runge-Kutta method, the stages on the test equation satisfy $Y=e{y}_n+zAY$. Whenever these stages are solvable, its stability function isThis rational function is used to test scalar A-stability; the Butcher order conditions separately establish nonlinear order. The Lobatto IIIA method in Question 3 is an A-stable fourth-order example with $R$ the diagonal degree-two Padé approximant. More generally the $s$-stage Gauss--Legendre Runge-Kutta method has order $2s$ and is A-stable, with the diagonal Padé approximation to $e^z$; its stiff-limit factor is nonzero. The $s$-stage Radau IIA method has order $2s-1$ and is L-stable. Thus Runge-Kutta stages permit arbitrarily high A-stable order, unlike irreducible linear multistep formulas.

For a normal matrix, scalar amplification bounds transfer through unitary diagonalization. For a non-normal matrix, eigenvalue information alone does not assert Euclidean contraction: transient growth and eigenvector conditioning can matter. Direct energy method arguments, such as the dissipative Cayley-transform proof in Question 1, use the symmetric part and are more informative in that setting.

For nonlinear dissipative vector fields, the stronger B-stability property asks that distances between two numerical solutions not increase. A standard sufficient condition for Runge-Kutta methods is algebraic stability: $b_i\ge 0$ andis positive semidefinite. It follows from the Runge-Kutta contractivity identity when the implicit stages exist. The three-stage Lobatto IIIA method is A-stable but has the first diagonal entry of that matrix equal to $-1/36$, so scalar A-stability should not be confused with that sufficient nonlinear condition. Stage solvability of an implicit Runge-Kutta method remains a separate issue and computational cost of the implicit equations must be considered. A-stability removes a scalar decay-mode stability restriction; order, stiff damping, nonlinear stability and solvability remain distinct properties.

A Runge-Kutta method advances $y^{\prime }=f(t,y)$ through stagesThe Butcher tableau records $A,b,c$, normally with $c=Ae$. Linear absolute stability tests the amplification of a linear decay mode; nonlinear B-stability compares distances between two solutions of a dissipative vector field. Both describe propagation of perturbations, but they impose different conditions.

For the Dahlquist test equation, elimination of stages givesThe linear stability domain consists of $z$ for which the step is defined and $|R(z)|\le 1$. Poles or singular stage systems must be excluded, even if an output formula has a formal cancellation. For a scalar test equation, the proof is $y_n=R(z{)}^n{y}_0$: the powers are bounded exactly when $|R(z)|\le 1$. For $y^{\prime }=Ly$, the update is $R(hL)$. If $L$ is a normal matrix, unitary diagonalization of a normal matrix makes the operator norm equal to the largest scalar amplification modulus. For a diagonalizable matrix, powers are bounded by the condition number of the eigenvector matrix times the maximal scalar power. These constants must remain uniform for mesh-dependent systems. A defective unit-modulus eigenvalue can instead produce a growing Jordan block; the scalar spectral test alone does not control a general matrix.

An A-stable method accepts the entire closed left half-plane. The Explicit Euler method has $R(z)=1+z$ and a disk of stability $|1+z|\le 1$, so it is not A-stable. No nonconstant polynomial stability function can be A-stable, since it is unbounded on the negative real axis; this excludes all consistent explicit Runge-Kutta methods. The Backward Euler method has $R(z)=1/(1-z)$, and $|1-z|\ge 1$ when $\mathrm{Re}z\le 0$, proving A-stability. It is also L-stable, because $R(z)\to 0$ for large left-half-plane $z$.

The implicit midpoint rule and the trapezoidal rule have the same scalar stability function,The identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability. On the imaginary axis their amplification has unit modulus, and along the negative real axis it tends to $-1$. Thus neither is L-stable; very stiff decay can persist as alternating numerical values. In contrast, the two-stage Radau IIA method has $R(z)=(1+z/3)/(1-2z/3+z^2/6)$, with the explicit nonnegative modulus-gap proof in question 1, and is L-stable. This illustrates why L-stability matters beyond A-stability for a stiff differential equation.

For nonlinear problems, assume a real or complex Euclidean inner product and the dissipativity hypothesisfor all relevant states at each common time. Exact solutions satisfyA B-stable Runge-Kutta method reproduces this nonexpansion of distances whenever its stages exist. Put $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$, and $G_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,{\tilde{Y}}_i)$. Then $D_i=d+h{\sum }_j{a}_{ij}{G}_j$. Expansion of the output squared norm, followed by substituting $d=D_i-h{\sum }_j{a}_{ij}{G}_j$ in its linear terms, proves the Runge-Kutta contractivity identityIf $b_i\ge 0$ and $M=(m_{ij})$ is a positive semidefinite matrix, the first sum is nonpositive by dissipativity, and the second is nonnegative: expand each coordinate of the stage vectors to express it as a sum of nonnegative quadratic forms. These are exactly the conditions for algebraic stability of a Runge-Kutta method. Hence algebraic stability implies B-stability. Applying B-stability to the linear dissipative problem $y^{\prime }=\lambda y$, or its two-dimensional real form, also proves A-stability whenever stages are well defined.

The Backward Euler method has $b=A=1$, hence $M=1$; the implicit midpoint rule has $b=1$, $A=1/2$, hence $M=0$. Both are algebraically stable and B-stable. For the two-stage Radau IIA method,a positive semidefinite matrix, while both weights are positive. It therefore combines third-order accuracy, B-stability, and L-stability. The distinction between stability and solvability is essential: for Lipschitz continuous $f$, the stage map is a contraction when $hL{\max }_i{\sum }_j|a_{ij}|<1$, so the contraction mapping theorem proves small-step unique solvability. The preceding algebraic stability proof is an estimate for existing stages; it is not a theorem asserting stage existence for arbitrary step sizes and domains.

Linear A-stability is insufficient for nonlinear B-stability. The trapezoidal rule fails B-stability, even though its scalar stability function is identical to that of the B-stable implicit midpoint rule. Its Runge-Kutta method coefficients areFailure of this sufficient algebraic stability test alone would not prove failure of B-stability. An actual counterexample does: take $f(y)=-y^3$ and $h=2$. Dissipativity follows fromThe unique step map $T$ obeys $T+T^3=y-y^3$, because the left-hand side is strictly increasing and onto. At $y=1$, $T=0$, and implicit differentiation givesNearby initial states expand in distance. This proves a true nonlinear failure rather than merely failure of a coefficient criterion.

Finally, stability of a numerical method on finite intervals connects these contractivity properties to accuracy. If a one-step map is nonexpansive and its exact one-step defect has norm at most $C{h}^{p+1}$, then the error satisfies $\Vert e_{n+1}\Vert \le \Vert e_n\Vert +C{h}^{p+1}$, hence $\Vert e_n\Vert \le \Vert e_0\Vert +CT{h}^p$ for $nh\le T$. A Lipschitz step bound $1+Kh$ gives the same order using the discrete Gronwall inequality. For a stiff differential equation, the value of B-stability is that the propagation estimate need not grow with a large negative dissipative rate. The accuracy constant still needs appropriate smoothness and uniform derivative bounds; stability alone does not establish a mesh-uniform error order.

An $s$-stage Runge-Kutta method for $y^{\prime }=f(t,y)$ is defined byThe Butcher tableau records $A=(a_{ij})$, $b$ and $c$; internal consistency usually sets $c=Ae$, and first-order consistency requires $b^{T}e=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^{\prime }=\lambda y$. Solving the stages givesThe linear stability domain consists of test values for which the stage system is well defined and $|R(z)|\le 1$. A-stability means this includes the closed left half-plane. For $\mathrm{Re}\lambda <0$, this gives stability with no scalar decay-mode step restriction. For a normal matrix $L$ in $y^{\prime }=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|\le 1$, and restriction $0\le hq\le 2$ for $\lambda =-q<0$. For a pure imaginary mode $\lambda =i\omega \ne 0$, $|1+ih\omega |>1$, so it is unstable for every positive step. The classical fourth-order Runge-Kutta method instead hasIt is stable on the imaginary axis exactly for $|q|\le 2\sqrt{2}$, 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)\to 0$ 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 isDifferentiating the squared difference and applying the Gronwall inequality gives $\Vert y(t)-\tilde{y}(t)\Vert \le e^{\nu (t-t_0)}\Vert y(t_0)-\tilde{y}(t_0)\Vert$. A dissipative vector field has $\nu \le 0$ and is contractive. A B-stable method preserves this property for every positive step size for which the stages are well defined: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 algebraic stability implies B-stability, subject to stage solvability. Define $B=\mathrm{diag}(b)$ and $M=BA+A^{T}B-b{b}^T$. Algebraic stability of a Runge-Kutta method means $b_i\ge 0$ and $M$ is a positive semidefinite matrix. To prove the theorem, set $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$ and $F_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,{\tilde{Y}}_i)$. Expanding the output difference and using $D_i=d+h{\sum }_j{a}_{ij}{F}_j$ yields the Runge-Kutta contractivity identityDissipativity makes the first sum nonpositive, while positive semidefiniteness makes the last quadratic form nonnegative. The latter follows by factoring $M=Q^{T}Q$ 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,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, withIts positive weights and positive semidefinite $M$ prove B-stability; its rational function has denominator roots $2\pm i\sqrt{2}$ and, for $z=x+iy$, satisfiesThus it is A-stable; the limit $R(z)\to 0$ 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, A-stability alone does not imply B-stability. The trapezoidal rule fails B-stability even for the scalar dissipative equation $y^{\prime }=-y^3$. With $h=1$, its step map $T$ is uniquely defined byAt $y=\sqrt{2}$, $T(y)=0$. Implicit differentiation givesThus 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.

A Runge-Kutta method with stage matrix $A$, weights $b$ and step size $h$ isIts order measures accuracy as $h\to 0$, while stability concerns how errors or stiff components propagate. Linear A-stability, strong stiff damping through L-stability, and nonlinear B-stability address different questions.

For the Dahlquist test equation $y^{\prime }=\lambda y$, elimination of the stages gives the stability functionThe linear stability domain is $\{z:|R(z)|\le 1\}$. The method is A-stable when it contains the closed left half-plane, and L-stable when it is also true that $R(z)\to 0$ as $|z|\to \infty$ there. These conditions describe, respectively, stability for every non-growing scalar linear mode and damping of very rapidly decaying modes.

The Forward Euler method has $R(z)=1+z$, so it fails A-stability, for example at $z=-3$. In fact no nontrivial explicit Runge-Kutta method is A-stable: its stability function is a nonconstant polynomial, which is unbounded along the negative real axis. The Backward Euler method has $R(z)=1/(1-z)$. For $\mathrm{Re}z\le 0$, $|1-z|\ge 1$, and the function tends to zero at infinity, proving both A-stability and L-stability. The implicit midpoint rule has $R(z)=(1+z/2)/(1-z/2)$. The identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability, but its limit $-1$ rules out L-stability.

For nonlinear equations, a dissipative vector field satisfies $⟨f(x)-f(y),x-y⟩\le 0$. The corresponding exact solutions contract because the time derivative of their squared distance is twice that inner product. A B-stable method preserves this contraction for every $h>0$, whenever its stages are well defined.

Let $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$, and $F_i=f(Y_i)-f({\tilde{Y}}_i)$. Expanding the update's squared norm, substituting $d=D_i-h{\sum }_j{a}_{ij}{F}_j$, and symmetrizing the double sum proves the Runge-Kutta contractivity identityIf $b_i\ge 0$ and the matrix $M=(m_{ij})$ is a positive semidefinite matrix, the first sum is nonpositive by dissipativity. The second double sum is nonnegative: sum the quadratic form of $M$ over each coordinate of the $F_i$. Thus algebraic stability implies B-stability. Applied to complex scalar linear equations with the real part of the Hermitian inner product, the same argument also implies A-stability, where the stages are defined.

The Backward Euler method has $b_1=a_{11}=1$, hence $M=(1)$, and is algebraically stable. One can prove its nonlinear contraction directly: $d_{n+1}=d+h(f(y_{n+1})-f({\tilde{y}}_{n+1}))$ gives $\Vert d_{n+1}{\Vert }^2\le ⟨d,d_{n+1}⟩\le \Vert d\Vert \Vert d_{n+1}\Vert$. The implicit midpoint rule has $a_{11}=1/2$, $b_1=1$ and $M=0$, so it is B-stable even though it is not L-stable.

Linear A-stability does not by itself imply nonlinear B-stability. The ODE trapezoidal rule has the same stability function as the implicit midpoint rule, but take the dissipative scalar field $f(y)=-\max (y,0)$ and $h=6$. For any positive starting value, its trapezoidal update is $y_{n+1}=y_n+3(f(y_n)+f(y_{n+1}))=-2y_n<0$. Thus two positive starting values have their distance doubled, violating B-stability. Its algebraic stability matrix also has a negative first diagonal entry. The field is globally Lipschitz continuous, so this is a well-defined nonlinear counterexample.

These examples distinguish accuracy, scalar stiff-mode stability, nonlinear contraction, and rapid stiff-mode damping. For implicit methods, solvability of the stage equations remains an additional requirement; a formal stability function alone does not establish it for arbitrary nonlinear problems.

Apply the implicit midpoint rule, equivalently the Crank-Nicolson method, to the semidiscrete equation:It has order two. Its amplification matrix is the Cayley transformBecause $H_h$ is Hermitian, $Q$ is unitary. Hence $\Vert {\mathbf{u}}^{n+1}{\Vert }_2=\Vert {\mathbf{u}}^n{\Vert }_2$ exactly.

An $s$-stage Runge-Kutta method has stages and updateOn the Dahlquist test equation $y^{\prime }=\lambda y$, elimination of the stages gives the stability functionThe linear stability domain is the set where $|R(z)|\le 1$. A method is A-stable when this domain contains $\mathrm{Re}z\le 0$, 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)\to 0$ 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)|\le 1$ 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 thatA method is B-stable if it preserves the resulting contractivity: two numerical solutions satisfy $\Vert y_{n+1}-{\tilde{y}}_{n+1}\Vert \le \Vert y_n-{\tilde{y}}_n\Vert$. A practical sufficient condition is algebraic stability of a Runge-Kutta method: $b_i\ge 0$ andTo prove the implication, let $D_i=Y_i-{\tilde{Y}}_i$ and $F_i=f(Y_i)-f({\tilde{Y}}_i)$. Expanding the squared distance and substituting the stage equations gives the Runge-Kutta contractivity identityThe 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^{\prime }=\lambda y$, implies A-stability.

Important collocation families illustrate these notions. Gauss methods are A-stable, symmetric, and have order $2s$, but they do not damp infinitely stiff modes. 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.

Set $\overline{\mathbf{y}}=({\mathbf{y}}_n+{\mathbf{y}}_{n+1})/2$. The implicit midpoint rule readsTaking the inner product with ${\mathbf{y}}_{n+1}+{\mathbf{y}}_n=2\overline{\mathbf{y}}$ givesThe last equality again follows from skew-symmetry. Mathematical induction therefore yields