Embedded Runge-Kutta pair

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.

Continue the initial value problem $y^{\prime }=y$, $y(0)=1$, to $t=4$ using the Euler-Heun embedded pair. Choose absolute and relative error tolerances $\mathrm{atol}=10^{-3}$ and $\mathrm{rtol}=0$, so the scalar scale is $s=10^{-3}$ throughout. This intentionally absolute-only test makes the local and accumulated errors directly comparable. Start with $h=0.2$, accept the Heun's method endpoint precisely when $E=|d|/s\le 1$, and after every trial use the capped adaptive step-size controllerFor $E=0$ use the growth factor two. Cap every trial at $4-t_n$; on rejection retain the old time and state. These choices specify an original computation, rather than a claim about the defaults of a particular software solver.

The initial trial has $E=20$ and is rejected. Retrying from $(0,1)$ with $h=0.9(0.2)/\sqrt{20}\approx 0.0402492236$ gives $E=0.81$ and is accepted. The run takes 315 accepted steps and one rejected trial. Increasing exponential growth makes the local error estimator $|d|=|y_n|h^2/2$ grow at a fixed $h$, so subsequent steps shrink roughly like $|y_n{|}^{-1/2}$. The final step is additionally shortened to land at $t=4$.

The right panel plots errors against the accepted-step number, all divided by the same $s$. The actual Heun local defect uses the restarted exact endpoint $y_n{e}^{h_n}$, while actual global error uses $e^{t_{n+1}}$. All accepted $E$ are at most one, and the actual Heun local defects here are even smaller, at most about $0.01143s$. Nevertheless the endpoint is $y_N\approx 54.5833912205$, versus $e^4\approx 54.5981500331$, with absolute error about $0.0147588=14.7588s$. Passing local tests therefore does not certify the original trajectory's global error.

This failure of a global guarantee has an exact explanation for this example. Put $R(h)=1+h+h^2/2$ and $g_n=e^{t_n}-y_n$. ThenEach new local defect is added to the previous error amplified by $e^{h_n}$. Iterating this identity, rather than merely inspecting $E$, accounts for accumulation. The variable-step global ODE error bound describes the corresponding mechanism for more general ordinary differential equations.

Repeating with $\mathrm{atol}=10^{-2}$ gives 98 accepted steps, one rejection and endpoint absolute error about $0.153252$; with $\mathrm{atol}=10^{-4}$ it gives 1001 accepted steps, two rejections and endpoint absolute error about $0.00145642$. Accuracy improves on this smooth problem, but the endpoint error remains about fifteen times the requested local scale. A different initial value problem can amplify or damp defects differently.