Runge-Kutta 4

Deep dive

How it works

In: f, t, y, fixed step h · Out: y(t + h) with local error O(h⁵) · State: none between steps (four stage vectors as scratch) · Cost: 4 f-evaluations and about 10n flops per step

For y' = f(t, y), step h:

k1 = f(t,       y)
k2 = f(t + h/2, y + h/2·k1)
k3 = f(t + h/2, y + h/2·k2)
k4 = f(t + h,   y + h·k3)
y_next = y + h/6·(k1 + 2k2 + 2k3 + k4)

Butcher tableau: c = (0, ½, ½, 1), a₂₁ = ½, a₃₂ = ½, a₄₃ = 1, b = (⅙, ⅓, ⅓, ⅙). For f depending only on t it reduces to Simpson's rule.

Error: local O(h⁵), global O(h⁴). Stability: applied to y' = λy one step multiplies y by R(hλ), the degree-4 Taylor polynomial of e^z. Every eigenvalue of the Jacobian times h must land where |R| ≤ 1.

Code

/// Fixed-size systems: stack arrays, no allocation; the compiler unrolls small N.
pub fn rk4_step<const N: usize>(f: impl Fn(f64, &[f64; N]) -> [f64; N], t: f64, y: &[f64; N], h: f64) -> [f64; N] {
    let axpy = |a: &[f64; N], b: &[f64; N], s: f64| -> [f64; N] { std::array::from_fn(|i| a[i] + s * b[i]) };
    let k1 = f(t, y);
    let k2 = f(t + 0.5 * h, &axpy(y, &k1, 0.5 * h));
    let k3 = f(t + 0.5 * h, &axpy(y, &k2, 0.5 * h));
    let k4 = f(t + h, &axpy(y, &k3, h));
    std::array::from_fn(|i| y[i] + h / 6.0 * (k1[i] + 2.0 * (k2[i] + k3[i]) + k4[i]))
}

/// Large systems (method of lines, particle sets): reuse four stage buffers.
pub struct Rk4 { k: [Vec<f64>; 4], tmp: Vec<f64> }

impl Rk4 {
    pub fn new(n: usize) -> Self { Rk4 { k: std::array::from_fn(|_| vec![0.0; n]), tmp: vec![0.0; n] } }

    pub fn step<F: FnMut(f64, &[f64], &mut [f64])>(&mut self, f: &mut F, t: f64, y: &mut [f64], h: f64) {
        let [k1, k2, k3, k4] = &mut self.k;
        let tmp = &mut self.tmp;
        f(t, y, k1);
        for i in 0..y.len() { tmp[i] = y[i] + 0.5 * h * k1[i]; }
        f(t + 0.5 * h, tmp, k2);
        for i in 0..y.len() { tmp[i] = y[i] + 0.5 * h * k2[i]; }
        f(t + 0.5 * h, tmp, k3);
        for i in 0..y.len() { tmp[i] = y[i] + h * k3[i]; }
        f(t + h, tmp, k4);
        for i in 0..y.len() { y[i] += h / 6.0 * (k1[i] + 2.0 * (k2[i] + k3[i]) + k4[i]); }
    }
}

/// Cubic Hermite interpolation between two accepted steps (y0, f0 at t0; y1, f1 at t0 + h).
/// Third-order accurate: good for plotting and event location, one order below the step.
pub fn hermite(y0: f64, f0: f64, y1: f64, f1: f64, h: f64, theta: f64) -> f64 {
    let (t2, t3) = (theta * theta, theta * theta * theta);
    (2.0 * t3 - 3.0 * t2 + 1.0) * y0 + (t3 - 2.0 * t2 + theta) * h * f0
        + (-2.0 * t3 + 3.0 * t2) * y1 + (t3 - t2) * h * f1
}

/// Largest stable |h·lambda| along the negative real axis, found by bisection on |R(z)| = 1.
pub fn rk4_real_stability_limit() -> f64 {
    let r = |z: f64| 1.0 + z + z * z / 2.0 + z.powi(3) / 6.0 + z.powi(4) / 24.0;
    let (mut lo, mut hi) = (-2.0, -3.0); // |R(-2)| < 1 < |R(-3)|
    for _ in 0..60 {
        let mid = 0.5 * (lo + hi);
        if r(mid).abs() <= 1.0 { lo = mid } else { hi = mid }
    }
    lo
}

Complexity and constants

• 4 f-evaluations + about 10n flops of vector updates per step; for n large and f cheap it is memory-bound (5 vectors touched per stage).
• Harmonic oscillator y'' = −y on [0, 10]: h = 0.1 gives max error ≈ 6.7·10⁻⁶, h = 0.05 gives ≈ 4.1·10⁻⁷ (ratio 16.07).
• Amplitude loss per step for an oscillator of frequency ω: ≈ (ωh)⁶/144 (half of the |R|² defect). At ωh = 0.1: about 7·10⁻⁹ per step, which is still a secular drift over 10⁸ steps.
• Versus Euler at equal cost (4 Euler steps of h/4): for smooth problems RK4 is many orders of magnitude more accurate.

Efficiency

• Cost: 4 f-evaluations plus about 10n flops of vector updates per step; for cheap f and large n it is memory-bound (5 vectors touched per stage). Measured with the code above on y'' = −y over [0, 10]: h = 0.1 gives max error 6.7·10⁻⁶ and h = 0.05 gives 4.1·10⁻⁷ (ratio 16.07). Bound: cost of f, or memory bandwidth for cheap f.
• Levers:
  • larger h at 4th order — error ÷ 16 per halving, so 4 Euler steps of h/4 cost the same and lose by orders of magnitude (page above) — smooth problems
  • low-storage RK (Williamson 2N) — two registers per variable instead of five — huge PDE state vectors
  • Hermite cubic interpolation for outputs — no extra f-evaluations — several outputs per step
  • fused stage updates — fewer memory passes — large n with cheap f
  • h at the stability limit (|hλ| < 2.78) with margin — fewest steps that stay stable — stiff-ish modes
• Naive → standard → state of the art:
  • explicit Euler at tiny h · error O(h), 1 evaluation per step
  • fixed-step RK4 · 4 evaluations per step, error 4.1·10⁻⁷ at h = 0.05 (measured above)
  • embedded Dormand-Prince 5(4) with step control · 6 evaluations per step, tolerance-driven h (adaptive RK)

Implementation notes

• Allocation-free: stack arrays for small N, preallocated stage buffers for large N. The f(t, y, &mut out) signature avoids per-call vectors.
• For y'' = a(x, v) (mechanics), you can still use RK4 on the first-order system (x, v); Verlet is cheaper (1 force evaluation) and conserves energy better over long runs.
• If you need multiple outputs per step, interpolate with Hermite cubic (above) instead of shrinking the step.
• Choose h from the stability limit: max over eigenvalues of |hλ| below 2.78 (decay) or 2.82 (oscillation), with margin.
• Fixed-step RK4 is deterministic given identical floating-point operations, which lockstep simulations rely on; avoid FMA contraction differences across targets.

Pitfalls

• Using RK4 for orbit or molecular dynamics over long times: energy slowly decays (or grows outside the stability region) instead of oscillating around the true value.
• Updating y in place while computing stages: k2..k4 must all use the original y.
• Accumulating t by t += h for millions of steps: drift; compute t0 + i·h.
• Assuming 4th order with non-smooth f (friction switches, contacts, abs, min): the order drops to 1 at each kink; step to the event or smooth it.

Variants and alternatives

In production

Numerical libraries: Boost.Odeint runge_kutta4, GSL gsl_odeiv2_step_rk4, DifferentialEquations.jl RK4. Robotics and physics: MuJoCo's RK4 integrator option; many game simulators for rigid-body attitude and trajectory prediction. Machine learning: torchdiffeq's fixed-grid rk4 solver. Scientific computing: method-of-lines PDE codes.

Implementations

Sources

• Ernst Hairer, Syvert P. Nørsett, Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993. book
• John C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd ed., Wiley, 2016. book
• Chi-Wang Shu, Stanley Osher, Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes, Journal of Computational Physics 77, 1988. paper
• Emanuel Todorov, Tom Erez, Yuval Tassa, MuJoCo: A physics engine for model-based control, IROS 2012. paper · code