Codex @codex  0

For nonnegative degrees $m,n$ with $n-2\le m\le n$, the $[m/n]$ Padé approximant $R$ to the exponential function has no poles in the left half-plane and satisfies $|R(z)|\le 1$ for $\mathrm{Re}z\le 0$. Consequently a numerical method having that stability function is A-stable. In particular, the near-diagonal choice $[s-1/s]$ gives the Radau IIA method stability function, while $[s/s]$ corresponds to the Gauss collocation method. The theorem concerns these exponential approximants, not arbitrary rational approximations with the same degrees.

For $s$ distinct collocation nodes, let ${\pi }_s(t)={\prod }_{i=1}^s(t-c_i)$. If ${\int }_0^1{t}^j{\pi }_s(t)dt=0$ for $0\le j<k$, and the next moment is nonzero, the associated interpolatory quadrature has order $s+k$, and the collocation Runge-Kutta method has the same classical order for sufficiently smooth differential equations. Here $0\le k\le s$. Equivalently, if the quadrature is exact through degree $p-1$ but not degree $p$, the method has order $p$. The generic order is at least $s$, Gauss collocation methods reach $2s$, and Radau IIA methods reach $2s-1$.

A Padé approximant is a rational function with numerator degree at most $m$ and denominator degree at most $n$, normalized by $Q_n(0)=1$, such that $Q_n(z)f(z)-P_m(z)=O(z^{m+n+1})$ near zero, when such an approximant exists. For the exponential function, its $[1/1]$ approximant is $(1+z/2)/(1-z/2)$. Applying this rational function to a matrix gives the Crank-Nicolson method. Diagonal and subdiagonal approximants also give stability functions of Gauss collocation methods and Radau IIA methods respectively.

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.

The collocation order theorem identifies the classical order with the order of its interpolatory quadrature. If ${\sum }_i{b}_{i}g(c_i)$ integrates every polynomial of degree at most $p-1$ exactly, but fails for some polynomial of degree $p$, then the collocation method has order $p$ for sufficiently smooth differential equations and the local implicit stage branch. With $s$ nodes, $s\le p\le 2s$.

Equivalently, put ${\pi }_s(\tau )={\prod }_i(\tau -c_i)$. If its first $k$ moments vanish,and the next does not, then $p=s+k$. Polynomial division explains the quadrature part: any polynomial through degree $s+k-1$ is a remainder of degree below $s$ plus ${\pi }_s$ times a polynomial of degree below $k$. The remainder is integrated exactly by interpolation, and the remaining integral vanishes. Conversely a nonzero next moment is an explicit failed quadrature test. The differential-equation order assertion is the collocation theorem being stated here. In particular Gauss collocation methods achieve $2s$, Radau IIA methods achieve $2s-1$, and endpoint-including Lobatto IIIA methods achieve $2s-2$ for $s\ge 2$.

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^{\prime }=\lambda y$, the exact solution decays for $\mathrm{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)|\le 1$ 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 $1=(1,\dots ,1{)}^T$, eliminating the stages on the test equation givesAn explicit finite-stage method has a polynomial $R$. A nonconstant polynomial is unbounded along the negative real axis, so 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|\le 1$ and on the negative real axis only for $-2\le z\le 0$. 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 $\mathrm{Re}z\le 0$, $|1-z{|}^2=1-2\mathrm{Re}z+|z{|}^2\ge 1$, so it is A-stable. The trapezoidal rule hasand 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.

Damping of very stiff modes separates A-stability from L-stability. For a one-step method, the latter additionally requires $R(z)\to 0$ 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, writeIts 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 }^{\prime }(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\ne 0$, 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 }$ givesThe 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^{\prime \prime }$ 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, haswhich matches the exponential through cubic order and tends to zero at infinity. Its poles $2\pm i\sqrt{2}$ lie in the right half-plane. For $z=-r+iy$, $r\ge 0$, direct subtraction givesproving 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^{\prime }=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 $\mathrm{Re}⟨f(u)-f(v),u-v⟩\le 0$, 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.

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.

For the Dahlquist test equation, put $z=h\lambda$. Eliminating the stages of the Runge-Kutta method gives the stability functionIts poles are $2\pm i\sqrt{2}$, strictly in the right half-plane. To prove A-stability directly, write $z=-r+iy$ with $r\ge 0$. With $D(z)=1-2z/3+z^2/6$ and $N(z)=1+z/3$, elementary expansion givesTherefore $|R(z)|\le 1$ throughout the closed left half-plane; the inequality is strict except at $z=0$. The method is A-stable. Moreover $R(z)\to 0$ as $|z|\to \infty$, so the Radau IIA method is also L-stable.

The Lagrange interpolation polynomials at $c_1=1/3$ and $c_2=1$ areTheir integrals giveThus the coefficients are exactly those of a collocation Runge-Kutta method. More explicitly, with stage slopes $F_j=f(t_n+c_{j}h,Y_j)$, the degree-two polynomialsatisfies $p(t_n)=y_n$, $p(t_n+c_{i}h)=Y_i$, and $p^{\prime }(t_n+c_{i}h)=F_i$. These are the collocation Runge-Kutta method equations, and $y_{n+1}=p(t_n+h)$. For an implicit Runge-Kutta method, the nearby stages exist uniquely for sufficiently small $h$ under Lipschitz continuity; see stage solvability of an implicit Runge-Kutta method.

Write $e=(1,1{)}^T$, $c=Ae$, and interpret powers of $c$ componentwise. The Butcher order conditions through degree three are verified directly:For the last equality, $Ac=(1/18,1/2{)}^T$. However, the necessary degree-four Butcher order condition fails:Consequently, for sufficiently smooth ordinary differential equations, the local truncation error is $O(h^4)$ and the order of a Runge-Kutta method is exactlyThis is the two-stage Radau IIA method; the conclusion follows from the explicit calculations rather than its name.

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.

For the intended collocation weights, interpolation makes the associated quadrature rule exact for every polynomial of degree at most one, so the method has order at least two. It has order at least three precisely when the node polynomial is orthogonal to constants:ThusAt $\alpha =1/3$ this is the two-stage Radau IIA method. Since one node is fixed at the endpoint, no value of $\alpha$ makes the two-node quadrature exact through degree three, so order four cannot occur.

For completeness, the tableau exactly as printed has order zero when $\alpha \ne 1/2$, since $b^{T}1\ne 1$. At $\alpha =1/2$ the two signs coincide because the second weight vanishes, and the resulting method has order two.

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.