Adaptive Step Size Control with Embedded Runge-Kutta Pairs

An integrator with a fixed step has to be told the step before it starts, which means the user must know how the solution behaves along the trajectory. The alternative is a controller: estimate the error of the step just taken, compare it to a requested tolerance, and choose the next step from that ratio. The mechanism is small enough to state exactly. The details, namely how the estimate is obtained for free, why the exponent in the update is not 1, and what happens at the ends of the range, decide whether the controller is stable or oscillates.

Two approximations from one set of stages

The error estimate must not double the cost of the step. The embedded pair achieves this by computing the stage values once and combining them with two different weight vectors, one of order $p$ and one of order $\hat p = p - 1$. The tableau’s stage values $k_i$ are already available, so the second combination costs $s$ extra multiply-adds per component, not $s$ extra function evaluations.

The local error estimate is the difference of the two combinations:

$$ \hat e = h \sum_{i=1}^{s} \left(b_i - \hat b_i\right) k_i $$

Its leading term is $O(h^{p})$, one order above the local truncation error of the lower-order method and one below that of the higher, in the usual arrangement where the higher order is propagated. This matters for the controller’s exponent, since the estimate’s scaling with step is the quantity that appears in the step update rule.

The difference-of-weights form has a consequence worth stating: any stage whose two weights coincide contributes nothing to the estimate. Pairs are constructed so that the free parameters are used to make the estimate good rather than to save arithmetic. The two classical families are distinguished by what they optimize. Fehlberg’s pairs minimize the error of the lower-order combination, on the theory that the estimate should be a good approximation to the true error. Dormand-Prince pairs optimize the propagated higher-order combination, on the theory that what the integrator propagates is what should be accurate. The choice matters in long integrations where the accumulated error is dominated by the propagated solution, and it is why one family is favored for problems with modest tolerances and the other for problems where the propagated solution’s quality is the goal.

Implementing the estimate and the reject loop is short:

/* One adaptive step with an embedded RK pair; weights in b_hi, b_lo. */
double dt = dt0, t = t0;
while (t < t_end) {
    rk_stages(f, t, y, dt, k);                    /* fills s stage vectors */
    for (int c = 0; c < n; c++) {
        double hi = 0.0, lo = 0.0;
        for (int i = 0; i < s; i++) {
            hi += b_hi[i] * k[i][c];
            lo += b_lo[i] * k[i][c];
        }
        y_hi[c] = y[c] + dt * hi;
        err[c]  = dt * (hi - lo);
    }
    double e = wrms_norm(err, y, y_hi, n, atol, rtol);
    if (e <= 1.0) {                               /* step accepted */
        memcpy(y, y_hi, n * sizeof *y);
        t += dt;
    }
    dt *= fmin(FACMAX, fmax(FACMIN, SAFETY * pow(e, -1.0 / (order + 1))));
    if (dt > dt_max) dt = dt_max;
}

fmin and fmax clamp the growth and shrinkage factors; without the clamps a step that lands far inside the tolerance produces a wild increase, and a step that fails badly produces a wild decrease that saturates the arithmetic. Typical values are a safety factor a little below 1, a maximum growth of about 5 per step and a maximum shrinkage of about 0.1 per rejected step.

Why the exponent is not one

The error estimate scales with the step as $O(h^{q})$ with $q = \hat p + 1$, the order of the estimate’s leading term. If the current step gives an error $e$ and the target is a tolerance normalized to 1, the step that would give exactly the target satisfies

$$ h_{\text{new}} = h \left(\frac{1}{e}\right)^{1/q} $$

so the exponent is $1/q$, not 1. For a fourth-order propagated solution with a third-order error estimate, $q = 4$ and the exponent is $1/4$. This is the source of a common early bug: applying a linear scaling makes the controller overshoot, and the overshoot shows up as an alternating pattern of accepted and rejected steps rather than a smooth adaptation.

The estimate is asymptotic, so its predicted scaling is only valid near the tolerance. The safety factor below 1 exists for the same reason. A controller that requests exactly the target error has a roughly even chance of exceeding it on the next step because the local error varies along the solution, and at its $O(h)$ accuracy the estimate itself carries error. A factor of 0.8 to 0.9 is the conventional buffer.

Normalizing the error: relative, absolute, and combined

The tolerance must be on a scale-free number, or a step in a solution of magnitude $10^{6}$ and a step in a solution of magnitude $10^{-6}$ cannot both be controlled by the same request. The standard normalization combines relative and absolute components:

scale_i = atol + rtol * max(|y_i|, |y_new_i|)
err_i   = e_i / scale_i
e       = sqrt( (1/n) * sum_i err_i^2 )

The maximum of the old and new values in the denominator is deliberate: using the new value alone makes the scale drop when a component decays toward zero, which inflates the normalized error and produces a spurious rejection. The root-mean-square over components means a single component out of tolerance does not force a rejection on its own unless it is far outside, a choice that differs from the infinity norm used in some solvers. The two norms give different behavior in practice: the RMS allows a component to sit slightly outside the tolerance while the worst-case component stays within, and the maximum norm does not.

This norm belongs with the notion of tolerance as a per-component scale rather than a global error bound. The request names a magnitude for each component, and the controller treats a solution whose component $i$ is large and a solution whose component $i$ is near zero with the same normalization; the same tolerance therefore means different things in different regions of the solution’s range. Other integrators normalize by the initial value or by the maximum over the integration, which produces a constant denominator and a controller that does not react to a component decaying through zero. The choice is part of the solver’s contract rather than an internal detail, and it changes which steps get rejected.

Saturation and the end of the interval

Two boundary effects deserve handling. The first is the end of the integration interval: a step that would carry $t$ past $t_{\text{end}}$ must be truncated to land exactly there, and a truncation by a tiny amount can produce a step far smaller than the controller’s choice, with a correspondingly smaller error and a large growth request that is immediately discarded. The usual approach is a final step that satisfies the tolerance without being artificially cut.

The second is tolerance saturation. The local error estimate has a floor set by roundoff, roughly the machine epsilon times the solution scale, and a tolerance below that floor cannot be met. The controller then reduces the step, the rejection repeats, and the integration stalls. A guard is to check for a step below a minimum fraction of the elapsed interval or of the machine epsilon’s square root times the solution scale, and to report the stall rather than reduce further. The floor exists for the reasons given in local truncation error: truncation and roundoff contributions move in opposite directions with the step and cross at a nonzero optimum.

What the controller does not guarantee

The local error estimate bounds the error of a single step, not the accumulated error. The global error can be much larger than the tolerance times the number of steps, and it can be dominated by an amplification factor tied to the problem’s sensitivity, as the derivation in global error accumulation makes explicit. Tightening the tolerance tightens the local estimate and, for a well-behaved problem, tightens the global error proportionally, but a problem with exponential divergence has global error that grows independently.

The controller also does not handle stiffness. Its response to a stiff mode is to reduce the step until the stability requirement is met, which is exactly the collapse described in stiffness; the estimate becomes tiny and the step grows again, and the cycle repeats at the stability limit rather than at the accuracy limit. Detectable symptom: the accepted step tracks the inverse of the Jacobian spectral radius rather than the solution’s variation.

When the estimate can be wrong

An embedded estimate compares two different methods, and it is possible for their errors to cancel to leading order on a particular problem, making the estimate far smaller than the true error. This is not hypothetical: it happens systematically when the right-hand side is a low-degree polynomial in the state along a stretch of the trajectory, and it can cause a step to grow without bound until the higher-order terms dominate. The standard guard is a maximum step increase per step, which catches the runaway before it produces an inaccurate accepted step, and, for long integrations, periodic step-size reductions that force the method out of the degenerate regime.