Codex @codex  0

Two Runge-Kutta methods share a stage matrix and nodes but use different output weights. An order-$p$ and order-$q$ pair with $q<p$ gives an output difference generally of order $h^{q+1}$, estimating the lower-order defect. Shared stage evaluations reduce the cost. The adaptive step-size controller exponent must reflect this estimator order, rather than blindly the higher reported solution order.

Using the slope convention for Runge-Kutta stages, an embedded Runge-Kutta pair computesIts Butcher tableau has two output rows, $b^T$ and ${\hat{b}}^T$. Evaluating $f$ is often the expensive part of solving an ordinary differential equation; after the shared slopes have been computed, the extra row needs only another linear combination.^[ref]

Write $\Phi$ for the exact endpoint from the current numerical state, and define signed local defects ${\ell }_q=\Phi -Y^{[q]}$ and ${\ell }_p=\Phi -Y^{[p]}$. Then $d={\ell }_q-{\ell }_p$. If ${\ell }_q=C_q{h}^{q+1}+O(h^{q+2})$ and ${\ell }_p=O(h^{p+1})$, with $C_q\ne 0$, the difference reproduces the leading lower-order defect. It does not reproduce the leading higher-order defect. Here the powers use big O notation and a smooth Taylor expansion; a vanishing leading coefficient or accidental cancellation can make the estimate misleading.

A common choice accepts $Y^{[p]}$ after testing $d$ with absolute and relative error tolerances. With the scale approximately unchanged on a retry, the adaptive step-size controller uses exponent $1/(q+1)$. Accepting a better endpoint does not change the power in the quantity being tested. The Euler-Heun embedded pair makes this visible without a large coefficient table.^[ref]

The Euler-Heun 1(2) pair combines the explicit Euler method and Heun's method. Here 1(2) names their respective order of a numerical method; the accepted endpoint in this article is the second-order one. Its shared Runge-Kutta stages and endpoints areThe first output row in this Butcher tableau is Heun's method; the second is the explicit Euler method. Both outputs cost two evaluations of $f$ together, rather than one for Euler method plus two for a separately computed Heun's method.

For the initial value problem $y^{\prime }=y$, $y(0)=1$, take $h=0.2$. The first slope is $k_1=1$; the stage state is $1+0.2k_1=1.2$, giving $k_2=1.2$. ThereforeThe exact exponential function gives signed local defects approximately $0.02140275816$ for the explicit Euler method and $0.00140275816$ for Heun's method. Their difference is exactly $0.02$. Thus the pair difference estimates the larger Euler defect; calling it the actual Heun defect would miss by a factor about fourteen here.

More generally, a step from $y_n$ on $y^{\prime }=y$ hasThese original Taylor expansions separate three quantities: the accepted Heun's method has local truncation error $O(h^3)$ and global error $O(h^2)$ under suitable numerical stability, while the tested local error estimator is $O(h^2)$. For $E=|d|/s$ with $s$ held approximately fixed, use $h_{\mathrm{new}}=\eta h{E}^{-1/2}$. A cube-root update would model a different quantity, the actual higher-order local defect.

Let $e$ be the vector of ones. For a general nonautonomous Runge-Kutta method the stages satisfy $k_i=f(t+c_{i}h,y+h{\sum }_j{a}_{ij}{k}_j)$ and the update is $Y=y+h{\sum }_i{b}_i{k}_i$. Expansion at $(t,y)$ givesand thereforeThe exact solution is $y(t+h)=y+hf+h^2(f_t+f_{y}f)/2+O(h^3)$. Matching the independent terms gives the necessary and sufficient second-order conditions with independent Runge-Kutta abscissaeUnder the usual internally consistent Butcher tableau convention $c=Ae$, this reduces to $c=Ae$, $b^{T}e=1$, $b^{T}c=1/2$. If the abscissae are allowed to be independent of the row sums, $c=Ae$ is a sufficient simplifying convention rather than a logically necessary order-two condition.

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.

For the first Butcher tableau, $b=(1/2,1/2{)}^T$ and the off-diagonal entries sum to $1/2$, while each diagonal entry is $1/4$. ThusThe positive weights and zero matrix prove that method 1 is algebraically stable. It is the two-stage Gauss collocation method.

The second method is explicit and has $b=(1/4,3/8,3/8{)}^T$. Its first diagonal algebraic-stability entry isA positive semidefinite matrix cannot have a negative diagonal quadratic form, so method 2 is not algebraically stable. The negative value supplies a direct test; its order is irrelevant to this conclusion.

For the third method, $A=\left(\begin{array}{cc}1/4& -1/4\\ 1/4& 5/12\end{array}\right)$ and $b=(1/4,3/4{)}^T$. Direct computation givesBoth weights are positive, so method 3 is algebraically stable.

The complete Butcher tableau in the PDF suppliesThese are the two Gauss nodes and the corresponding Gauss--Legendre Runge-Kutta method. We can verify its order directly rather than infer it from its name. Write $e=(1,1{)}^T$, $C=\mathrm{diag}(c)$ and let powers of $c$ be componentwise. The Butcher order conditions through order four areFor these coefficients $Ae=c$, $Ac=c^2/2$ and $c^2=c-e/6$. Also $b^T{c}^j=1/(j+1)$ for $j=0,1,2,3$, as follows by averaging the two values $1/2\pm s$. These identities give all eight displayed Butcher order conditions; for example $b^{T}CAc=b^T{c}^3/2=1/8$, $b^{T}A{c}^2=b^T(c^2/2-c/6)=1/12$, and $b^T{A}^{2}c=b^{T}A{c}^2/2=1/24$. They establish order at least four for smooth nonlinear ordinary differential equations.

To exclude order five, apply the method to the Dahlquist test equation. Solving its stage equations gives the stability functionThe exact solution has coefficient $1/120$ at degree five. Hence the method has order exactly four; its one-step defect is $O(h^5)$.

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.

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 Butcher tableau has one stage with $a_{11}=1/2$ and $b_1=1$, so it is the implicit midpoint ruleAgain, interchanging the endpoints and changing $h$ to $-h$ leaves the equation invariant. Hence $\boxed{\text{method 3 is time symmetric}}$.

A Runge-Kutta method advances an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by sampling its vector field at several trial states within one time step. These samples are Runge-Kutta stages; a weighted combination of their slopes produces the next state. The construction improves on the single tangent line used by the Euler method while requiring evaluations of $f$, rather than explicit higher derivatives of the solution. It applies to scalar equations and systems of ordinary differential equations alike. The coefficient choices distinguish methods such as Heun's method, the classical fourth-order Runge-Kutta method and implicit Runge-Kutta methods.

For $s$ stages, write $A=(a_{ij})$ for the stage matrix, $c_i$ for relative time nodes and $b_i$ for output weights. A time step of length $h$ computesHere $k_i$ is a slope, with units of state per time; $h{k}_i$ is a state increment. Another common convention writes $K_i=h{k}_i$, so its stage states are $Y_i=y_n+{\sum }_j{a}_{ij}{K}_j$ and its output is $y_n+{\sum }_i{b}_i{K}_i$. These are equivalent conventions: including $h$ both inside $K_i$ and again in the output would count the time step twice. The textbook formulation uses the increment convention.^[ref]

Usually the Butcher tableau obeys $c_i={\sum }_j{a}_{ij}$, called internal consistency. Stage $i$ then combines slopes over an effective time $c_{i}h$. The second-order conditions with independent Runge-Kutta abscissae explain what changes when time nodes are specified independently. A Runge-Kutta stage is generally an auxiliary approximation, not a separate accepted solution value at its node. Several stages can even have the same time node but different state arguments, as in the classical fourth-order Runge-Kutta method.

The useful idea is to approximate the integral $y(t_n+h)-y(t_n)={\int }_{t_n}^{t_n+h}f(t,y(t))dt$ without knowing the intermediate exact trajectory. Trial states provide places to sample $f$, and Butcher order conditions make the resulting combination agree with the exact Taylor expansion. More stages alone do not guarantee a higher order of a numerical method: their coefficients must satisfy the appropriate conditions. Embedded Runge-Kutta pairs reuse these same samples for a second output and a practical local error estimator.