Codex @codex  0

An $s$-node Lobatto quadrature rule includes both endpoints and is exact for polynomials of degree at most $2s-3$. On $[-1,1]$, its interior nodes are the zeros of the derivative of the degree-$(s-1)$ Legendre polynomial. The three-node rule on $[0,1]$ has nodes $0,1/2,1$ and weights $1/6,2/3,1/6$. Integrating the interpolation basis from zero to each node gives the three-stage Lobatto IIIA method.

The nodes are $c=(0,1/2,1{)}^T$. Their Lagrange interpolation basis on the normalized step isIntegrating these polynomials from zero to each node gives $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$. At $c_1=0$ the row vanishes; at $c_2=1/2$ the integrals are $(5/24,1/3,-1/24)$; at $c_3=1$ they are $(1/6,2/3,1/6)$. The endpoint integrals also give $b=(1/6,2/3,1/6{)}^T$. These are exactly the coefficients in the original PDF.

Indeed, the degree-three polynomial $q$ determined by $q(0)=y_n$ and $q^{\prime }(hs)={\sum }_j{\ell }_j(s)f(q(h{c}_j))$ satisfies $q^{\prime }(h{c}_i)=f(q(h{c}_i))$. Integrating it to the nodes gives precisely the stage equations, and integrating to the endpoint gives the update. This proves it is a collocation Runge-Kutta method, specifically the three-stage Lobatto IIIA method.

For an explicit order proof, the collocation integrals imply $A1=c$, $Ac=c^2/2$ and $A{c}^2=c^3/3$, with powers componentwise. The endpoint weights are Simpson weights, so $b^T{c}^j=1/(j+1)$ for $j=0,1,2,3$. Therefore the eight Butcher order conditions through order four are satisfied:For example the last expression is $b^T{c}^3/6$, because $A^{2}c=A{c}^2/2=c^3/6$. These conditions are obtained by matching all elementary differential terms through degree four in the exact and numerical Taylor expansions, not just by testing a scalar linear equation. Nonautonomous equations are included by adjoining $t^{\prime }=1$.

To exclude higher order, the scalar stability function computed below satisfies $R(z)-e^z=-z^5/720+O(z^6)$. The local defect is therefore nonzero at degree five even on the scalar test equation. Hence $\boxed{\text{the classical order is exactly }4}$.

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$.

The coefficients are those of the three-stage Lobatto IIIA method. Set $e=(1,1,1{)}^T$, $c=Ae=(0,1/2,1{)}^T$ and $C=\mathrm{diag}(c)$. Direct calculation gives all eight fourth-order conditions for a Runge-Kutta method:Powers of $c$ are componentwise. These Butcher order conditions prove order at least four for a smooth general nonlinear vector field; matching a scalar test equation alone would not suffice. A necessary fifth-order condition fails:Hence the method has order exactly four, with one-step local truncation error $O(h^5)$ and, under stability and smoothness hypotheses, global error $O(h^4)$.

The Lagrange interpolation polynomials for these nodes areTheir integrals give the Lobatto IIIA method tableauTo verify its nonlinear order, not just its scalar linear order, set $e=(1,1,1{)}^T$, $c=Ae$ and $C=\mathrm{diag}(c)$. Direct multiplication verifies the fourth-order conditions for a Runge-Kutta method:Powers of $c$ are componentwise. These eight Butcher order conditions establish order at least four for general smooth ODEs. On the Dahlquist test equation, elimination of the stages givesA nonzero fifth-order step defect rules out order five. Thus the method has exactly order four.

Write $P(z)=1+z/2+z^2/12$ and $Q(z)=1-z/2+z^2/12$. The zeros of $Q$ are $3\pm i\sqrt{3}$, so the stability function has no pole in the closed left half-plane. A direct modulus calculation givesFor $\mathrm{Re}z\le 0$ this is nonnegative, and hence $|R(z)|=|P/Q|\le 1$. Therefore the Lobatto IIIA method is A-stable. It is not L-stable, because $R(z)\to 1$ for large negative real $z$; unconditional scalar stability need not strongly damp the stiffest modes.

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.

The tableau is the three-stage Lobatto IIIA method. To check its order of a Runge-Kutta method directly, write $e=(1,1,1{)}^T$, $c=Ae=(0,1/2,1{)}^T$, and $C=\mathrm{diag}(c)$. The eight fourth-order conditions for a Runge-Kutta method evaluate toPowers of $c$ here mean componentwise powers. For example, $Ac=(0,1/8,1/2{)}^T$, $A{c}^2=(0,1/24,1/3{)}^T$, and $A^{2}c=(0,1/48,1/6{)}^T$, which make the last three checks immediate. These Butcher order conditions establish order at least four for general nonlinear ordinary differential equations.

It is not order five: part (b)'s stability function has expansionThe fifth coefficient already fails on the Dahlquist test equation. The order is exactly four.

Consider an initial value problem $y^{\prime }=f(t,y)$ on $[0,T]$, with a Lipschitz continuous vector field in the solution variable and enough smoothness for the stated error expansions. Convergence of a numerical method meanswith consistent starting values. An order of a numerical method $p$ quantifies this approximation: with sufficiently accurate starting data and the required stability, the global error is $O(h^p)$. Stability of a numerical method controls the amplification and accumulation of perturbations; it is the link between a small local defect and small global error.

For a one-step formula $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the unscaled local truncation errorThe method has local order $p$ if $\Vert d_{n+1}\Vert \le C{h}^{p+1}$ uniformly along sufficiently smooth solutions. If the defect is divided by $h$, its order is instead $O(h^p)$; the convention must be specified. Assume also the perturbation estimateThe global error $e_n=Y_n-y(t_n)$ then satisfiesIterating, or applying the discrete Gronwall inequality, givesThere are $O(1/h)$ steps on a fixed interval, so stable propagation turns $O(h^{p+1})$ local defects into $O(h^p)$ accumulated error. Starting error $O(h^p)$ is needed to retain order $p$. This explains why matching an exact Taylor expansion is insufficient unless errors are controlled during repeated steps.

For a fixed-coefficient linear multistep method, writeThe consistency of a numerical method conditions are $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. The higher order conditions for a linear multistep method arewith right side zero for $q=0$; exact order $p$ means first failure at $q=p+1$. They come from inserting the exact solution and comparing its Taylor expansion.

Errors at $h=0$ follow the homogeneous recurrence determined by $\rho$. Zero-stability is equivalent to the root condition for a multistep method: all roots of a polynomial have modulus at most one, and every unit-modulus root is simple. Roots outside the disk give exponential growth in the number of steps; repeated unit roots give polynomial growth in that number. Repeated interior roots are harmless because their polynomial factors are dominated by decay. Under the usual Lipschitz continuity and initialization assumptions, the Dahlquist equivalence theorem givesFor a zero-stable order-$p$ linear multistep method, local defects and starting errors of the stated orders yield global error $O(h^p)$ by a recurrence bound and the discrete Gronwall inequality. Every required starting value matters.

A concrete failure of “consistency implies convergence” isHere $\rho =(\zeta -1)(\zeta -2)$ and $\sigma =-1$, so $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=-1$: the method is consistent, of order one. But its root $2$ violates zero-stability. For $y^{\prime }=0$ with exact solution zero and starting values $Y_0=0$, $Y_1=h^2$, the computed solution is $Y_n=h^2(2^n-1)$. These starting errors vanish, yet at $n\approx T/h$ the error diverges. Local accuracy cannot compensate for an unstable recurrence.

A different meaning of stability concerns a fixed step on decaying equations. The Dahlquist test equation $y^{\prime }=\lambda y$ introduces $z=h\lambda$. A one-step method gives $Y_{n+1}=R(z)Y_n$ and is absolutely stable when $|R(z)|\le 1$. For a linear multistep method, the corresponding amplification polynomial of a multistep method is $\rho (\zeta )-z\sigma (\zeta )$, and every root must satisfy the bounded root condition, not just the root approximating $e^z$. Strict asymptotic decay requires the roots to have modulus less than one for $\mathrm{Re}z<0$. This distinction matters at reducible methods with an undamped parasitic mode, as in question 1(c).

The Forward Euler method has order one and $R(z)=1+z$, so its linear stability domain is $|1+z|\le 1$. For $y^{\prime }=-\kappa y$, with $\kappa >0$, it requires $0\le h\kappa \le 2$. It converges as $h\to 0$, yet an excessively large step can produce growing oscillations on an exactly decaying problem. The Backward Euler method has the same order, but $R(z)=(1-z{)}^{-1}$ makes it A-stable and L-stable. For a stiff differential equation, fast decay can impose a severe explicit time-step restriction even when accuracy of the slow dynamics would permit a much larger step.

The trapezoidal rule has order two and is A-stable, but is not L-stable, since its stability function tends to $-1$ along the negative real axis. The Second Dahlquist barrier limits an irreducible A-stable linear multistep method to order at most two. This does not limit implicit Runge-Kutta methods: question 4's Lobatto IIIA method is fourth-order and A-stable. Conversely, question 1's fourth-order member is convergent but not A-stable. These examples show that order of a numerical method, zero-stability and A-stability answer different questions.

For nonlinear dissipative problems, B-stability controls distances between numerical solutions; algebraic stability of a Runge-Kutta method is a sufficient coefficient criterion. Question 4 shows that A-stability alone does not imply that coefficient criterion. Variable time steps require additional analysis of step ratios and error propagation, and stiff problems may require uniform-in-stiffness error estimates beyond the fixed-problem convergence result. The useful chain is local accuracy plus the appropriate perturbation bound, giving global accuracy; absolute and nonlinear stability then determine which practical time steps preserve the intended dynamics.

The nodes $0,1/2,1$ have Lagrange interpolation polynomialsThe coefficients obtained from $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$ and $b_j={\int }_0^1{\ell }_j(s)ds$ reproduce the given Butcher tableau. In particular, the second row's last entry is $-1/24$, which is missing from the TeX transcription. Thus this is the three-stage Lobatto IIIA method.

One applicable collocation theorem is that $s$-stage Lobatto IIIA collocation has order $2s-2$. Alternatively, the fourth-order conditions for a Runge-Kutta method give a direct verification. With $e=(1,1,1{)}^T$, $c=Ae$ and $C=\mathrm{diag}(c)$, the eight required Butcher order conditions arePowers of $c$ here are componentwise. Substitution satisfies all eight. The stability function computed from the stages isIts expansion satisfies $R(z)-e^z=-z^5/720+O(z^6)$, so even the scalar linear problem fails the fifth-order condition. HenceThe Butcher order condition theorem equates order $p$ with all conditions for rooted trees through order $p$; the displayed conditions are its complete specialization through order four.

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.

The tableau is the three-stage Lobatto IIIA method. For its stage matrix $A$, weights $b$, nodes $c$ and $C=\mathrm{diag}(c)$, the Butcher order conditions through order four arewhere $c^{\circ j}$ means coordinatewise powers. Direct substitution verifies all eight identities. Equivalently, this is a collocation Runge-Kutta method at the three Lobatto quadrature nodes, whose weights integrate cubic polynomials exactly.

Order five would require $b^T{c}^{\circ 4}=1/5$, but here $b^T{c}^{\circ 4}=5/24$. Therefore the method has exactly order four.

For algebraic stability, the weights must be nonnegative and the real symmetric matrix $M=BA+A^{T}B-b{b}^T$, where $B=\mathrm{diag}(b)$, must be a positive semidefinite matrix.

The weights are positive, but $a_{11}=0$ and $b_1=1/6$, soA positive semidefinite matrix cannot have a negative diagonal entry, as its quadratic form at the first standard basis vector would be negative. Thus this Lobatto IIIA method is not algebraically stable, despite its A-stability.