Adaptive Runge-Kutta (Dormand-Prince)

Deep dive

How it works

In: f, t₀, y₀, t_end, rtol, atol, initial h · Out: y at t_end (accepted steps, dense output on request) · State: t, y, h, k₁ (reused from the last stage), previous error for the PI controller · Cost: 6 f-evaluations per accepted step; steps ∝ tol^(-1/5)

1. Stages k₁..k₇ from the Dormand-Prince tableau (code below). y₅ = y + h·Σ bᵢkᵢ, where b equals the 7th row of A, so k₇ = f(t + h, y₅).
2. Error vector e = h·Σ (bᵢ − b̂ᵢ)kᵢ, where b̂ are the 4th-order weights.
3. err = √(mean((eᵢ/scᵢ)²)) with scᵢ = atol + rtol·max(|yᵢ|, |y₅ᵢ|).
4. err ≤ 1: accept (t += h, y = y₅, k₁ ← k₇). Otherwise reject and retry with smaller h.
5. Next h from the controller. I-controller: h·0.9·err^(−1/5). PI-controller (smoother, fewer rejections): h·0.9·err^(−α)·err_prev^(β), with Hairer's DOPRI5 defaults β = 0.04, α = 0.2 − 0.75β = 0.17. After a rejection, never grow h on the next step.
6. Initial h: Hairer's heuristic from ‖y₀‖, ‖f₀‖ and one extra evaluation estimating ‖y''‖.

Code

const C: [f64; 7] = [0.0, 1.0 / 5.0, 3.0 / 10.0, 4.0 / 5.0, 8.0 / 9.0, 1.0, 1.0];
const A: [[f64; 6]; 7] = [
    [0.0; 6],
    [1.0 / 5.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [3.0 / 40.0, 9.0 / 40.0, 0.0, 0.0, 0.0, 0.0],
    [44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0, 0.0, 0.0, 0.0],
    [19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0, 0.0, 0.0],
    [9017.0 / 3168.0, -355.0 / 33.0, 46732.0 / 5247.0, 49.0 / 176.0, -5103.0 / 18656.0, 0.0],
    [35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, 11.0 / 84.0], // = b (5th order)
];
/// b(5th) - b(4th): the embedded error estimate.
const E: [f64; 7] = [71.0 / 57600.0, 0.0, -71.0 / 16695.0, 71.0 / 1920.0, -17253.0 / 339200.0, 22.0 / 525.0, -1.0 / 40.0];

pub struct Opts { pub rtol: f64, pub atol: f64, pub h_max: f64, pub max_steps: usize }
#[derive(Default, Debug)]
pub struct Stats { pub accepted: usize, pub rejected: usize, pub evals: usize }

fn rms_scaled(v: &[f64], y: &[f64], o: &Opts) -> f64 {
    (v.iter().zip(y).map(|(a, b)| (a / (o.atol + o.rtol * b.abs())).powi(2)).sum::<f64>() / v.len() as f64).sqrt()
}

pub fn dopri5<F: FnMut(f64, &[f64], &mut [f64])>(
    mut f: F, t0: f64, y0: &[f64], t1: f64, o: &Opts, mut on_step: impl FnMut(f64, &[f64]),
) -> Result<(Vec<f64>, Stats), String> {
    let n = y0.len();
    let (mut t, mut y) = (t0, y0.to_vec());
    let mut k = vec![vec![0.0; n]; 7];
    let mut ys = vec![0.0; n];
    let mut st = Stats::default();
    f(t, &y, &mut k[0]);
    // Initial step (Hairer, Norsett, Wanner, Sec. II.4).
    let (d0, d1) = (rms_scaled(&y, &y, o), rms_scaled(&k[0], &y, o));
    let h0 = if d0 < 1e-5 || d1 < 1e-5 { 1e-6 } else { 0.01 * d0 / d1 };
    for i in 0..n { ys[i] = y[i] + h0 * k[0][i]; }
    f(t + h0, &ys, &mut k[1]);
    let diff: Vec<f64> = (0..n).map(|i| (k[1][i] - k[0][i]) / h0).collect();
    let d2 = rms_scaled(&diff, &y, o);
    let h1 = if d1.max(d2) <= 1e-15 { (h0 * 1e-3).max(1e-6) } else { (0.01 / d1.max(d2)).powf(0.2) };
    let mut h = (100.0 * h0).min(h1);
    st.evals = 2;
    let (mut err_prev, mut rejected_last) = (1e-4f64, false);
    let (beta, alpha) = (0.04, 0.2 - 0.75 * 0.04);
    while t < t1 {
        if st.accepted + st.rejected >= o.max_steps { return Err(format!("max_steps reached at t = {t}")); }
        h = h.min(o.h_max).min(t1 - t);
        if h <= 10.0 * f64::EPSILON * t.abs().max(1.0) { return Err(format!("step size underflow at t = {t}")); }
        for s in 1..7 {
            for i in 0..n { ys[i] = y[i] + h * (0..s).map(|j| A[s][j] * k[j][i]).sum::<f64>(); }
            f(t + C[s] * h, &ys, &mut k[s]); // after s = 6, ys holds the 5th-order solution
        }
        st.evals += 6;
        let err = ((0..n)
            .map(|i| {
                let e = h * (0..7).map(|j| E[j] * k[j][i]).sum::<f64>();
                (e / (o.atol + o.rtol * y[i].abs().max(ys[i].abs()))).powi(2)
            })
            .sum::<f64>() / n as f64).sqrt();
        if err <= 1.0 {
            t = if t1 - t <= h { t1 } else { t + h }; // land exactly on t1
            y.copy_from_slice(&ys);
            k.swap(0, 6); // first same as last
            st.accepted += 1;
            on_step(t, &y);
            let fac = if err == 0.0 { 10.0 } else { 0.9 * err.powf(-alpha) * err_prev.powf(beta) };
            h *= if rejected_last { fac.clamp(0.2, 1.0) } else { fac.clamp(0.2, 10.0) };
            err_prev = err.max(1e-4);
            rejected_last = false;
        } else {
            h *= (0.9 * err.powf(-alpha)).max(0.2);
            st.rejected += 1;
            rejected_last = true;
        }
    }
    Ok((y, st))
}

Complexity and constants

• 6 evaluations per accepted step plus 6 per rejected one. Measured with the code above: y' = −y on [0, 10] at rtol 10⁻⁸ takes 91 steps, 0 rejections; the Arenstorf orbit (one period) at rtol = atol = 10⁻⁸ takes 338 accepted + 23 rejected steps (2,168 evaluations) and returns to its start within 5·10⁻⁷.
• Error per step ∝ h⁵ (local error of the 4th-order estimate), so steps scale like tol^(1/5): 10× tighter tolerance costs about 10^(1/5) ≈ 1.6× more steps.
• Overhead per step besides f: about 7·7·n multiply-adds for the stage sums and the error norm. For n in the millions this is memory-bound; fuse loops.
• Global error is not bounded by tol; it is usually within 1-2 orders of magnitude of it for non-chaotic problems. Check by halving tol and comparing.

Efficiency

• Cost: 6 f-evaluations per accepted step and 6 per rejected one (FSAL: 7 stages, first reused). Measured with the code above: y' = −y on [0, 10] at rtol 10⁻⁸ takes 91 steps; the Arenstorf orbit at rtol = atol = 10⁻⁸ takes 2,168 evaluations. Besides f, about 7·7·n multiply-adds per step and order of 9n doubles of state. Bound: cost of f; memory bandwidth when n is in the millions.
• Levers:
  • first-same-as-last stage — 1 of 7 evaluations saved per step (Dormand-Prince 1980) — always
  • PI step controller — fewer rejected steps than the I-controller (Gustafsson 1991; Hairer et al. 1993) — oscillating or noisy error estimates
  • looser tolerance — steps ∝ tol^(-1/5), so 10× looser is about 1.6× fewer steps (derived above) — when output accuracy allows
  • higher-order pair (DOP853) — fewer steps at tight tolerances (Hairer et al. 1993) — tol of order 10⁻⁷ and below
  • dense output instead of shrinking h to hit output times — no extra evaluations — fixed output grids and events
• Naive → standard → state of the art:
  • fixed-step RK4 with h set by the hardest region · 4 evaluations per step at the smallest h needed anywhere
  • Dormand-Prince 5(4), I-controller, scalar tolerance · 6 evaluations per step; 2,168 for the Arenstorf orbit (above)
  • PI controller, per-component atol, dense output, pair chosen by tolerance (Tsit5, Vern, DOP853) · same cost per step, fewer steps

Implementation notes

• Coefficients: verify A[s] row sums equal C[s] and E sums to 0 in a unit test; one wrong digit gives a silently first-order method.
• Mixed scales: atol is per component in production solvers (concentrations of 10⁻⁹ next to temperatures of 10³); a scalar atol here keeps the code short.
• Dense output: Dormand-Prince has a free 4th-order interpolant from the stored stages (as in Hairer's CONTD5); use it for output at fixed times and for event location instead of forcing small steps.
• Events (ground contact, zero crossing): detect a sign change of g(t, y) between steps, then root-find on the interpolant (Newton or Brent).
• Stiffness detection: Hairer's DOPRI5 estimates h·ρ(J) from the last two stages and warns when that estimate repeatedly exceeds 3.25, near this pair's stability boundary.
• Determinism: adaptive step sequences depend on every floating-point bit; different compilers or FMA contraction yield different step counts.

Pitfalls

• atol = 0 with a solution crossing zero: the scale becomes 0, err = ∞, step size collapses.
• Using the 4th-order solution to propagate: legal (Fehlberg did), but DP's coefficients are optimized for propagating the 5th.
• Growing h right after a rejection: oscillating accept/reject cycles.
• Non-smooth right-hand sides (switches, if on state): the controller shrinks h to tiny values at each discontinuity. Integrate to the event and restart.
• Tolerances near machine precision (rtol < 100·ε): err is dominated by rounding, the controller thrashes.
• Not reporting Stats: rejections and evaluation counts are the first diagnostic for stiffness or bad tolerances.

Variants and alternatives

In production

Engineering and numerical computing: MATLAB ode45, scipy solve_ivp(method='RK45') (default) and 'DOP853', Boost.Odeint runge_kutta_dopri5, DifferentialEquations.jl DP5, Hairer's Fortran DOPRI5/DOP853, GSL rk8pd/rkf45. Statistics: Stan's ode_rk45 is the recommended first solver for non-stiff ODE models (Dormand-Prince 1980 via Boost.Odeint). Machine learning: torchdiffeq (neural ODEs) uses dopri5 as its default solver.

Implementations

Sources

• John R. Dormand, Peter J. Prince, A Family of Embedded Runge-Kutta Formulae, Journal of Computational and Applied Mathematics 6(1), 1980. paper
• Ernst Hairer, Syvert P. Nørsett, Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993 (Sections II.4-II.6). book
• Lawrence F. Shampine, Mark W. Reichelt, The MATLAB ODE Suite, SIAM Journal on Scientific Computing 18(1), 1997. paper
• Kjell Gustafsson, Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods, ACM TOMS 17(4), 1991. paper
• Stan Development Team, Stan User's Guide: Ordinary Differential Equations (rk45 solver). docs
• Ricky T. Q. Chen et al., Neural Ordinary Differential Equations, NeurIPS 2018. paper · code