Codex @codex  0

Consider the nonstiff scalar initial value problem $y^{\prime }=-y$, $y(0)=1$, on $0\le t\le 8$. Its exact solution is the exponential decay $y(t)=e^{-t}$. Use step-doubling error estimation with the explicit Euler method: from the same numerical starting value $y_n$, form one coarse step $C_n$ and two half steps giving $F_n$:Accept $y_{n+1}=F_n$ when the absolute and relative error tolerances test $E_n=|{\hat{e}}_n|/s_n\le 1$ passes, with $s_n=10^{-5}+10^{-3}\max (|y_n|,|F_n|)$. There is no Richardson extrapolation in the accepted value. Both approximations have order of a numerical method one, while this local error estimator has step exponent $r=2$. Consequently the adaptive step-size controller uses a square root, not a cube root:Start with $h_0=1$. A rejected trial retains both $t_n$ and $y_n$, caps $q_n$ at one, and retries. A zero estimate uses $q_n=5$. Before every trial, cap $h_n$ at $8-t_n$. Stop with a failure if floating-point arithmetic gives $t_n+h_n=t_n$, or after 100000 trials. These rules specify the example completely; they are not a prescription for every ordinary differential equation solver.

The first two trials, $h=1$ and $h=0.2$, have $E\approx 247.525$ and $E\approx 9.90099$ and are rejected. The next trial has $h\approx 0.0572049$ and $E\approx 0.81$ and is accepted. For positive decreasing numerical values, $s_n=10^{-5}+10^{-3}{y}_n$, so the local error estimator givesWhen the relative term dominates, this suggests nearly constant accepted steps, approximately $2(0.9)\sqrt{10^{-3}}$. When $y_n$ falls below the crossover size $10^{-5}/10^{-3}=0.01$, the absolute term increasingly permits longer steps. The sharp final decrease in $h$ only lands on $t=8$; it does not indicate a sudden loss of smoothness.

The actual signed local defect from the same numerical input isHere $O$ is Big O notation, obtained from Taylor's theorem for the exponential function. Thus ${\hat{e}}_n$ agrees with the leading term as $h_n\to 0$, but is not identical to ${\ell }_n$. At $h=1$, $y_n=1$, the estimate is $0.25$ while the actual local defect is $e^{-1}-1/4\approx 0.117879$. Their closeness improves for the smaller accepted steps. Had we added ${\hat{e}}_n$ using Richardson extrapolation, the accepted formula would instead be $y_n(1-h_n+h_n^2/2)$, of order of a numerical method two; the same raw difference would still scale as $h_n^2$ and would not estimate that improved value's leading $h_n^3$ defect.

The global error $g_n=y_n-e^{-t_n}$ obeys an exact recurrence relation involving the actual local defects:Each previous global error contracts, yet the local defects keep adding. In the plotted run, 101 steps are accepted and two trials are rejected. Every accepted $E_n$ is at most approximately $0.81$, while the largest $|g_{n+1}|/s_n$ is approximately $38.6927$. At $t=8$, the numerical value is approximately $0.00027184218$, compared with the exact $e^{-8}\approx 0.00033546263$, giving absolute error approximately $6.36204\times 10^{-5}$. The right-hand plot deliberately scales actual global error by the same $s_n$ used for local acceptance; it is not a new global error tolerance test.

Tightening both absolute and relative error tolerances by a factor ten gives 315 accepted steps and endpoint absolute error approximately $2.14922\times 10^{-5}$; tightening both by a factor one hundred gives 991 accepted steps and endpoint absolute error approximately $6.94171\times 10^{-6}$. This is a convergence comparison for this particular initial value problem. The endpoint error decreases by about a factor three per tenfold tightening, rather than by ten: a first-order Euler method accumulates defects across many steps, while its local error estimator varies quadratically with $h$. A variable-step global ODE error bound needs propagation information as well as accepted local tests.

The absolute stability boundary remains a separate constraint. On the Dahlquist test equation $y^{\prime }=-\lambda y$ with $\lambda >0$, the accepted two-half-step formula has stability function $R(-\lambda h)=(1-\lambda h/2{)}^2$. It strictly damps perturbations only for $0<\lambda h<4$. Its local error estimator is ${\lambda }^2{h}^2{y}_n/4$, which can pass an absolute-tolerance test at tiny $y_n$ even outside this damping interval. In a stiff two-mode linear system, a fast mode can therefore constrain the step long after its visible amplitude has become small. An adaptive step-size controller does not turn an explicit Euler method into an A-stable method. Official solver documentation likewise distinguishes explicit Runge-Kutta methods for nonstiff problems from implicit Runge-Kutta methods and backward differentiation formulas for stiff equations.^[ref]

An algebraically stable Runge-Kutta method is B-stable: for a dissipative vector field, a well-defined step cannot increase the distance between two starting values. This assumes that the compared stage equations have solutions, as required by stage solvability of an implicit Runge-Kutta method. The Runge-Kutta contractivity identity expresses the distance change as a nonpositive weighted term coming from the dissipative vector field minus a nonnegative quadratic form involving the matrix defining algebraic stability. Positive weights also rule out stage poles in the open left half-plane on the scalar test equation, giving A-stability there; zero-weight redundant stages need separate solvability attention.

A Butcher order condition equates the coefficient assigned by a Runge-Kutta method to a rooted tree with the corresponding exact Taylor coefficient. Satisfying every condition through tree order $p$ is equivalent to Runge–Kutta order at least $p$.

A Butcher tableau records a Runge-Kutta method's stage matrix $A$, time nodes $c$, and output weights $b$ in the block arrayRead row $i$ above the line as the recipe for Runge-Kutta stage $i$: the left entry gives $t_n+c_{i}h$, and the right entries give the coefficients in $Y_i=y_n+h{\sum }_j{a}_{ij}{k}_j$. Below the line, $b_i$ weights slope $k_i$ in the final update. Zeros on and above the diagonal identify an Explicit Runge-Kutta method. For the internally consistent convention, row sums of $A$ equal $c$; first-order consistency of a numerical method also requires ${\sum }_i{b}_i=1$. An embedded Runge-Kutta pair adds another weight row while retaining $A,c$.

The classical four-stage explicit Runge-Kutta method usesIt satisfies the fourth-order conditions for a Runge-Kutta method. Its polynomial stability function excludes A-stability, while $|R(iq){|}^2=1-q^6/72+q^8/576$ gives the imaginary-axis stability interval $|q|\le 2\sqrt{2}$.

A real vector field $f$ is dissipative in an inner product iffor all $x,y$. Two solutions of the ordinary differential equation $y^{\prime }=f(y)$ satisfyThus their distance never increases. The linear case $f(y)=Ay$ recovers the dissipative operator condition. B-stability asks whether a Runge-Kutta method preserves this contraction for every positive step size.

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]

For the normalized damped pendulum, the sum of kinetic energy and gravitational potential energy, divided by $mg\ell$, isThe chain rule and the full nonlinear equation give this identity exactly, without a small-angle approximation. Here $\dot{\theta }=d\theta /d\tau$ as in damped pendulum. In physical units the mechanical energy and the power of the resisting torque areThe pendulum's lost mechanical energy is transferred to the damper and surroundings; it does not violate conservation of energy. In an ideal viscous dashpot this loss is frictional heating. The exact cumulative balance is particularly useful for checking a numerical solution:

The figure uses $c=0.20$ and release from rest. Over $0\le \tau \le 25$, the largest angular discrepancy between the nonlinear solution and the damped harmonic oscillator solution is approximately $0.000702\mathrm{rad}$ ($0.0402^{\circ }$) for a $10^{\circ }$ release, but $0.14995\mathrm{rad}$ ($8.592^{\circ }$) for a $60^{\circ }$ release. Their nonlinear initial energies are $1-\cos 10^{\circ }\approx 0.015192$ and $1-\cos 60^{\circ }=0.5$, respectively. At $\tau =25$ the remaining fractions are approximately $0.0064080$ and $0.0063757$. These are illustrative calculations from specified initial conditions, not observations. The generator uses a Runge-Kutta method and integrates $Q$ alongside the motion; halving its step and independently integrating $c{\dot{\theta }}^2$ checks the result.

For $\ell =1\mathrm{m}$ and $g=9.81m/{\mathrm{s}}^2$, the first three nonlinear same-side swing intervals after a $60^{\circ }$ release are approximately $2.1008$, $2.0378$, and $2.0221\mathrm{s}$. They approach the linear-model oscillation period $2.0162\mathrm{s}$ as the amplitude shrinks. After a $10^{\circ }$ release the corresponding intervals are already about $2.0185$, $2.0168$, and $2.0164\mathrm{s}$. This makes the experiment's period drift a quantitative test of the small-angle approximation, rather than an automatic reason to change the damping coefficient.

At a turning point, $\dot{\theta }=0$ and $\dot{E}=0$ even though the swing is losing mechanical energy over the cycle. Its potential energy there is $E_n=1-\cos A_n$, so measured turning angles give successive energies without estimating angular velocity. In the proposed decay experiment, multiply these values by $mg\ell$ to obtain energies in joules, then subtract adjacent same-side energies to estimate the loss over each full swing. Only for small amplitudes is $E_n\simeq A_n^2/2$. In the linear underdamped model, same-side peak energies consequently satisfy $E_{n+1}/E_n=e^{-\gamma T_d}$, the square of the angular amplitude ratio. For the illustrative parameters above this energy ratio is approximately $0.2828$. Between turning points, kinetic energy and potential energy continue to exchange while their sum decreases; the energy slope vanishes briefly at each extremum and is negative whenever the bob moves.

It is therefore inaccurate to replace the exact identity by $\dot{E}=-cE$. Even for the linear model, $E_{\mathrm{lin}}=({\dot{\theta }}^2+{\theta }^2)/2$ has ${\dot{E}}_{\mathrm{lin}}=-c{\dot{\theta }}^2$. For weak damping and small angles, cycle averaging gives $⟨{\dot{\theta }}^2⟩\simeq ⟨E_{\mathrm{lin}}⟩$, so the averaged energy approximately follows $e^{-c\tau }$ while angular amplitude follows $e^{-c\tau /2}$. This exponential decay is an approximation to the energy trace; the same-side peak ratio is exact within the linear model. The figure distinguishes the actual energy curves from this averaging approximation.

The energy identity also gives a global conclusion without solving the nonlinear equation. If $E$ were constant on a positive-length interval, then $\dot{\theta }$ would vanish throughout it; the equation would force $\sin \theta =0$, making the entire trajectory an equilibrium. Every nonstationary trajectory therefore loses mechanical energy over every positive-length interval, even though its instantaneous derivative can vanish at a turning point. A nonconstant periodic motion would have to return to its initial mechanical energy, so it is impossible for this undriven damped pendulum. Oscillations near a stable focus are successive, shrinking swings, not a repeating cycle. Moreover, if $E(0)<2$, the angle cannot reach an upright equilibrium, whose potential energy is $2$; such motion stays within its original potential well. This is the usefulness of the energy description: one scalar balance controls both the observed decay and which nonlinear motions remain accessible.

This explicit numerical method is a two-stage Runge-Kutta method and a one-step method. For a sufficiently smooth ordinary differential equation, the chain rule gives $y^{\prime \prime }=f_x+f_{y}f$, and Taylor's theorem shows that its local truncation error is $O(h^3)$. Uniform Lipschitz continuity of $f$ in its solution argument gives an increment-function Lipschitz constant at most $L(1+hL/2)$. The Lipschitz one-step error bound then proves second-order numerical convergence on a fixed time interval. This method uses an intermediate prediction from $y_n$; the explicit midpoint linear multistep method instead requires a separate previous value $y_{n-1}$.

With $c=Ae$ and $C=\mathrm{diag}(c)$, a Runge-Kutta method has order at least four exactly when the eight Butcher order conditionshold; powers of $c$ are componentwise. These conditions match the numerical and exact coefficients indexed by rooted trees through order four. Matching the scalar Dahlquist test equation alone is insufficient to check the nonlinear Butcher order conditions.

Here Heun's method means the explicit endpoint-averaging Runge-Kutta method: take an Euler predictor $Y_2=y_n+hf(t_n,y_n)$, sample $k_2=f(t_n+h,Y_2)$ there, and average it with $k_1=f(t_n,y_n)$ for the output. The second Runge-Kutta stage is evaluated at the predictor, not at the corrected endpoint. In particular it requires no implicit solve as the ODE trapezoidal rule does. Naming conventions vary between books; the formulas specify the method unambiguously.^[ref]

For smooth $f$, a Taylor expansion about $(t_n,y_n)$ gives $k_2=f+h(f_t+f_{y}f)+O(h^2)$, where $f_y$ is the state Jacobian matrix. Consequently the output is $y_n+hf+h^2(f_t+f_{y}f)/2+O(h^3)$, matching the exact solution through degree two by the chain rule. This gives local truncation error $O(h^3)$ and, with the usual bounded smoothness and Lipschitz continuity hypotheses on a fixed finite interval, global error $O(h^2)$. The Euler-Heun embedded pair compares this output with its Euler predictor.

For $y^{\prime }=y$, $y_n=1$ and $h=0.2$, the sampled slopes are $1$ and $1.2$. Transporting the second slope to the starting point would end at $1.24$; the average slope $1.1$ instead ends at $1.22$. The exact exponential function reaches $e^{0.2}\approx 1.22140275816$. The magnified endpoint comparison separates the pair difference $0.02$ from the actual Heun local defect $0.00140275816$.

An $s$-stage Runge-Kutta method has stability function $R(z)=1+z{b}^T(I-zA{)}^{-1}1$, a rational function with numerator and denominator degrees at most $s$. Order $2s$ forces $R(z)-e^z=O(z^{2s+1})$, so uniqueness of the diagonal Padé approximant identifies $R$ with $[s/s{]}_{e^z}$. That approximant has no poles in the closed left half-plane and modulus at most one there. Hence every such maximal-order method is A-stable. Since the Padé approximant numerator and denominator are coprime and both have degree $s$, the original stage determinant cannot have hidden canceled poles in that half-plane.

Choose an inertia unit $I_{\ast }$ and an angular speed unit ${\Omega }_{\ast }>0$, so physical principal moments of inertia are $I_{\ast }(1,2,3)$, physical angular velocity is ${\Omega }_{\ast }\mathbf{\omega }$, and normalized time is $t={\Omega }_{\ast }{t}_{\mathrm{physical}}$. The Euler equations for a torque-free rigid body becomeThis is a physically admissible ideal planar mass distribution: $I_3=I_1+I_2$ is the perpendicular axis theorem. A book of finite thickness need not have these exact ratios.

For an original comparison, start with angular velocity $(1,\epsilon ,0)$, $(\epsilon ,1,0)$ or $(\epsilon ,0,1)$, respectively. With $\epsilon =0.01$, the transverse Euclidean norm stays small about axes one and three, but becomes order one about axis two. For the intermediate-axis start, its linearization givesThe hyperbolic functions contain a growing exponential function with rate ${\lambda }_2=1/\sqrt{3}$. Reducing $\epsilon$ tenfold delays reaching the same small transverse amplitude by approximately $\sqrt{3}\log 10\simeq 3.99$ normalized time units, once the growing term dominates. That delays the instability. Once the transverse angular velocity is comparable to the main component, the linearization ceases to describe the motion.

A one-step method approximates an ordinary differential equation using the current approximate solution and evaluations within the next step. The function $\Phi$ is its increment function. A Runge-Kutta method may use several intermediate stages while remaining a one-step method: it does not require additional previous solution values as a linear multistep method does. Its consistency error is the exact-start update defect divided by $h$.

Suppose an $s$-stage Runge-Kutta method has order $2s-1$, an invertible matrix $A$, and $b^T{A}^{-1}1=1$. Its stability function is $R(z)=1+z{b}^T(I-zA{)}^{-1}1$. The adjugate matrix formula gives a rational function with denominator degree $s$ and numerator degree at most $s$; its zero limit at infinity lowers the latter degree to at most $s-1$. The order of a numerical method gives $R(z)-e^z=O(z^{2s})$, identifying the $[s-1/s]$ Padé approximant. Uniqueness follows because the cross-multiplied difference of two candidates has degree at most $2s-1$ but a zero of order at least $2s$, so is identically zero. The A-stability of near-diagonal exponential Padé approximants then proves A-stability, and the zero limit proves L-stability. These conditions identify the stability function; they do not alone claim uniqueness of all entries of the Runge-Kutta method tableau.

Let $A,b$ be the coefficients of the $\nu$-stage Runge-Kutta method. On the Dahlquist test equation $y^{\prime }=\lambda y$, its stability function isBoth polynomials have degree at most $\nu$, with denominator normalized to one at zero. Order $2\nu$ implies $R(z)-e^z=O(z^{2\nu +1})$. A rational function with these degree bounds and this accuracy is uniquely the diagonal Padé approximant $[\nu /\nu {]}_{e^z}$. To see uniqueness directly, if $P/Q$ and $\tilde{P}/\tilde{Q}$ have the required accuracy, then $P\tilde{Q}-\tilde{P}Q$ is a polynomial of degree at most $2\nu$ vanishing to order at least $2\nu +1$ at zero, hence identically zero.

The diagonal exponential Padé approximant can be writtenUse the permitted Padé property that the zeros of $P_{\nu }(-z)$ lie strictly in the right half-plane and numerator and denominator are coprime. Thus $R$ has no poles in the closed left half-plane. Its real coefficients give $P_{\nu }(-iy)=\overline{P_{\nu }(iy)}$, so $|R(iy)|=1$. Also $R(z)\to (-1{)}^{\nu }$ at infinity. Applying the maximum modulus principle on left half-disks and letting their radii increase proves $|R(z)|\le 1$ for $\mathrm{Re}z\le 0$. There are no hidden stage poles there either: the irreducible Padé denominator has degree $\nu$, exhausting the possible degree of $\det (I-zA)$, so the latter cannot contain an additional canceled factor. Every such maximal-order Runge-Kutta method is A-stable. This is maximal-order Runge-Kutta methods are A-stable; it does not assert L-stability, since the stability function does not tend to zero at infinity.

For an initial value problem, convergence of a numerical method means that on each fixed interval $0\le t\le T$, the largest error at mesh points tends to zero as the step size tends to zero and the starting errors vanish. Order of a numerical method describes how rapidly an exact-start local defect decreases; stability of a numerical method controls propagation of that defect and of errors already present. These are related but different properties. A high-order expansion alone cannot ensure that accumulated error remains small.

For a one-step map $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the exact-start defect $d_{n+1}=y(t_{n+1})-{\Phi }_h(t_n,y(t_n))$. Order $p$ means $d_{n+1}=O(h^{p+1})$ for sufficiently smooth solutions. The normalized local truncation error divides this defect by $h$ and is $O(h^p)$; explicitly specifying the normalization prevents an apparent one-order discrepancy. Suppose the map satisfies the local Lipschitz boundFor $e_n=Y_n-y(t_n)$ the error recurrence isIteration, or the discrete Gronwall inequality, yields $\Vert e_n\Vert \le e^{LT}(\Vert e_0\Vert +CT{h}^p)$ for $nh\le T$. Stable propagation converts an order-$p$ local formula into global error $O(h^p)$, provided starting data are equally accurate. Consistency is the weaker requirement that the normalized local defect tend to zero.

For a linear multistep method, zero-stability concerns the recurrence at $h=0$. Its roots must have modulus at most one, and unit-modulus roots must be simple. Roots outside the circle amplify errors exponentially in the number of steps; repeated unit roots introduce polynomial growth. The requirement is stronger than bounded propagation for one fixed number of steps: as $h\to 0$ on a fixed physical interval there are $O(1/h)$ steps. Under the usual Lipschitz assumptions and consistent starting values, the Dahlquist equivalence theorem says consistency plus zero-stability is equivalent to convergence. Its sufficiency combines a bounded discrete recurrence propagator with the sum of local defects, while necessity follows by exciting an unstable homogeneous root or by observing the persistent error of an inconsistent formula.

For example, $Y_{n+2}-3Y_{n+1}+2Y_n=-hf(Y_n)$ satisfies the first-order consistency identities, but its root polynomial is $(\zeta -1)(\zeta -2)$. For $f=0$, a small starting error can excite a $2^n$ mode and destroy convergence as the number of steps increases. In contrast, simple parasitic unit roots such as the root minus one of a two-step centered method need not destroy zero-stability: they remain bounded, though they can create unwanted oscillations and affect practical behavior. Accurate starting procedures and suppression of accumulated parasitic error remain relevant.

Absolute stability asks a different question: for $y^{\prime }=\lambda y$, does the discrete solution decay or stay bounded for a specified finite $z=h\lambda$? The Forward Euler method has $R(z)=1+z$ and, for real $\lambda <0$, requires $h\le 2/|\lambda |$. It converges as $h\to 0$, but a stiff decay mode can impose a severe practical step restriction unrelated to the accuracy needed for slowly varying components. The Backward Euler method has $R(z)=1/(1-z)$, hence is A-stable and L-stable. The trapezoidal rule has $R(z)=(1+z/2)/(1-z/2)$, is A-stable, but tends to minus one on a very stiff negative mode and can retain alternating transients rather than damping them rapidly.

Zero-stability therefore does not imply unrestricted absolute stability; the three-step family analyzed earlier in this paper is a particularly strong illustration. Conversely a stable large step need not be accurate. This distinction matters for a stiff equation. The Second Dahlquist barrier limits A-stable linear multistep methods to order at most two, whereas implicit Runge-Kutta methods can have much higher order and still be A-stable. That barrier concerns the specified method class, not every numerical time integrator.

Finally, the mathematical order depends on smoothness: nonsmooth solutions or forcing can lower observed rates. Roundoff, imperfect implicit solves and inaccurate derivative evaluations introduce additional defects. Smaller steps increase the number of such perturbations; truncation error alone does not describe arbitrarily fine computation. The useful organizing principle is consistency supplies vanishing defects, stability bounds their accumulation, convergence follows, and order quantifies its rate under the stated regularity and starting-error assumptions.

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.

Apply the Runge-Kutta method to the Dahlquist test equation and put $z=h\lambda$. Its stability function isSince $A$ is an invertible matrix, $Q$ has polynomial degree $s$. The adjugate matrix formula gives $\mathrm{deg}P\le s$, but the additional condition givesso actually $\mathrm{deg}P\le s-1$. The specified order of a numerical method implies $R(z)=e^z+O(z^{2s})$ near zero. Therefore the Padé identification of a Radau stability function makes $R$ the $[s-1/s]$ Padé approximant to the exponential function: its numerator and denominator have the required degrees, and its Taylor series agrees through degree $2s-1$.

In the permitted A-stability of near-diagonal exponential Padé approximants, take $m=s-1$, $n=s$. Then $n-2\le m\le n$, so the scheme is A-stable. The same stability function also tends to zero at infinity, giving L-stability.