Stiff solvers (implicit, BDF)

Deep dive

How it works

In: y' = f(t, y), y(t₀), tolerance, and a Jacobian (analytic or by finite differences) · Out: y at the requested times · State: past solution values, step size, order, and the cached factorization of I − βhJ · Cost: per step one Newton solve, 1-3 iterations; dense LU 2n³/3 flops amortized over many steps

BDF of order k fits a polynomial through y_{n+1}, y_n, …, y_{n+1−k} and requires its derivative at t_{n+1} to equal f:

BDF1 (backward Euler):  y₁ − y₀                    = h·f(t₁, y₁)
BDF2:                   y₂ − 4/3·y₁ + 1/3·y₀       = 2/3·h·f(t₂, y₂)

Each step solves G(y) = y − β·h·f(t, y) − r = 0 (r from past values, β = 1 or 2/3):

1. Predict by extrapolating past values (y_pred = 2y₁ − y₀ for BDF2).
2. Newton: solve (I − βhJ)·δ = −G(y), y += δ; stop when the scaled ‖δ‖ is a small fraction of the error tolerance.
3. Reuse: the iteration matrix I − βhJ changes slowly; factor it once and keep it across iterations and steps (simplified Newton). Refresh J when the contraction ‖δₖ‖/‖δₖ₋₁‖ is poor, refactor when h changes much.
4. Error control (production): compare with a predictor-based estimate, change h and order k.

Stability on y' = λy: backward Euler gives R(z) = 1/(1 − z), |R| < 1 for every Re z < 0 and R(∞) = 0 (L-stable). Trapezoidal gives R(∞) = −1: stiff modes flip sign each step instead of decaying (ringing).

Code

/// LU with partial pivoting, in place, row-major n×n. Returns None if singular.
fn lu_factor(a: &mut [f64], n: usize, piv: &mut [usize]) -> Option<()> {
    for k in 0..n {
        let p = (k..n).max_by(|&i, &j| a[i * n + k].abs().total_cmp(&a[j * n + k].abs()))?;
        if a[p * n + k] == 0.0 { return None; }
        piv[k] = p;
        if p != k { for j in 0..n { a.swap(k * n + j, p * n + j); } }
        for i in k + 1..n {
            let l = a[i * n + k] / a[k * n + k];
            a[i * n + k] = l;
            for j in k + 1..n { a[i * n + j] -= l * a[k * n + j]; }
        }
    }
    Some(())
}
fn lu_solve(a: &[f64], n: usize, piv: &[usize], b: &mut [f64]) {
    for k in 0..n { b.swap(k, piv[k]); for i in k + 1..n { b[i] -= a[i * n + k] * b[k]; } }
    for k in (0..n).rev() { for j in k + 1..n { b[k] -= a[k * n + j] * b[j]; } b[k] /= a[k * n + k]; }
}

/// Fixed-step BDF2 (first step BDF1) with Newton and a finite-difference Jacobian.
pub struct Bdf2 {
    n: usize, pub h: f64, pub rtol: f64, pub atol: f64, // Newton stops when |dz_i| < atol + rtol |z_i|
    prev: Option<Vec<f64>>, m: Vec<f64>, piv: Vec<usize>, pub newton_iters: usize,
}

impl Bdf2 {
    pub fn new(n: usize, h: f64) -> Self {
        Bdf2 { n, h, rtol: 1e-8, atol: 1e-14, prev: None, m: vec![0.0; n * n], piv: vec![0; n], newton_iters: 0 }
    }

    /// M = I - beta h J(t1, z), J by forward differences, then LU-factor M.
    fn iteration_matrix<F: FnMut(f64, &[f64], &mut [f64])>(
        &mut self, f: &mut F, t1: f64, z: &[f64], f0: &[f64], beta: f64,
    ) -> Result<(), &'static str> {
        let (n, h) = (self.n, self.h);
        let (mut zp, mut f1) = (z.to_vec(), vec![0.0; n]);
        for j in 0..n {
            let d = f64::EPSILON.sqrt() * z[j].abs().max(1.0);
            zp[j] = z[j] + d;
            f(t1, &zp, &mut f1);
            zp[j] = z[j];
            for i in 0..n { self.m[i * n + j] = -beta * h * (f1[i] - f0[i]) / d; }
        }
        for i in 0..n { self.m[i * n + i] += 1.0; }
        lu_factor(&mut self.m, n, &mut self.piv).ok_or("singular iteration matrix")
    }

    pub fn step<F: FnMut(f64, &[f64], &mut [f64])>(&mut self, f: &mut F, t: f64, y: &mut [f64]) -> Result<(), &'static str> {
        let (n, h, t1) = (self.n, self.h, t + self.h);
        let (beta, r, mut z): (f64, Vec<f64>, Vec<f64>) = match &self.prev {
            None => (1.0, y.to_vec(), y.to_vec()),
            Some(p) => (2.0 / 3.0,
                        (0..n).map(|i| 4.0 / 3.0 * y[i] - 1.0 / 3.0 * p[i]).collect(),
                        (0..n).map(|i| 2.0 * y[i] - p[i]).collect()), // linear extrapolation
        };
        let mut f0 = vec![0.0; n];
        f(t1, &z, &mut f0);
        self.iteration_matrix(f, t1, &z, &f0, beta)?;
        let (mut g, mut last) = (vec![0.0; n], f64::INFINITY);
        for _ in 0..20 {
            for i in 0..n { g[i] = -(z[i] - beta * h * f0[i] - r[i]); } // -G(z)
            lu_solve(&self.m, n, &self.piv, &mut g);                    // simplified Newton: M is reused
            self.newton_iters += 1;
            let mut dnorm = 0.0f64;
            for i in 0..n { z[i] += g[i]; dnorm = dnorm.max(g[i].abs() / (self.atol + self.rtol * z[i].abs())); }
            if dnorm < 1.0 {
                self.prev = Some(y.to_vec());
                y.copy_from_slice(&z);
                return Ok(());
            }
            f(t1, &z, &mut f0);
            if dnorm > 0.5 * last { self.iteration_matrix(f, t1, &z, &f0, beta)?; } // slow: refresh J
            last = dnorm;
        }
        Err("Newton did not converge: reduce h")
    }
}

Complexity and constants

• Per step (dense, as above): n + 1 f-evaluations for the Jacobian, 2n³/3 flops for the LU, 2n² per Newton iteration. Production codes amortize J and LU over many steps (order of 10-50), so the steady cost is a few f-evaluations plus back-solves.
• Large sparse systems (PDEs): sparse LU (KLU, SuperLU) or Newton-Krylov (GMRES with a preconditioner, Jacobian-vector products by finite differences) instead of dense LU.
• Example: Robertson's chemistry problem (rate constants 0.04, 10⁴, 3·10⁷). Its stiff eigenvalue is ≈ −3.4·10³ on [0, 40] (measured) and grows toward order of −10⁴ as y₃ → 1, so an explicit method needs h < 6·10⁻⁴ forever, even at t = 10¹⁰ when nothing changes. Implicit codes grow h with t. The fixed-step BDF2 above at h = 10⁻⁴ gives y₁(40) = 0.715827069 (reference 0.7158270687).
• Newton usually converges in 1-3 iterations per step from a good predictor (1.01 on average in the Robertson run above).

Efficiency

Cost: a dense step costs n + 1 f-evaluations for the Jacobian, 2n³/3 flops for the LU and 2n² flops per Newton iteration; production codes amortize J and LU over order of 10-50 steps, so the steady cost is a few f-evaluations plus back-solves. Bound: LU factorization for dense n (compute), memory traffic for sparse systems. On Robertson's problem (measured above) Newton needed 1.01 iterations per step.

Levers:

• Reuse J and the LU (simplified Newton) — 2n³/3 flops paid once per many steps — Jacobian changes slowly.
• Analytic or colour-compressed sparse Jacobian — evaluations drop from n to the number of colours — known sparsity.
• Variable step and order (CVODE-style, order 1-5) — step grows with t where nothing changes (Robertson at t = 10¹⁰) — long, smooth tails.
• Sparse LU (KLU, SuperLU) or Newton-Krylov — O(nnz)-ish memory instead of n² — large PDE systems.
• Rosenbrock methods — no Newton loop, one LU per step — moderate n.

Naive → standard → state of the art:

• Explicit RK with h < 6·10⁻⁴ for ever (Robertson) · cost scales with the fastest mode
• Fixed-step BDF2 with dense LU · stable at h = 10⁻⁴ (above)
• Variable-order BDF with LU reuse and sparse or Krylov linear solves (CVODE, ode15s) · step follows the solution, not the stiffest mode

Implementation notes

• Scaling: test Newton convergence in the same weighted norm (atol + rtol·|y|) as the error control.
• Analytic or sparse Jacobians: finite differences cost n evaluations; with a known sparsity pattern, column colouring (Curtis-Powell-Reid) needs only as many evaluations as colours.
• Mass matrices / DAEs: M·y' = f(y) changes the iteration matrix to M − βhJ; BDF handles index-1 DAEs directly (DASSL/IDA).
• Variable step BDF: coefficients depend on the step ratios; keep ratios moderate (production solvers limit growth per step) or use the fixed-leading-coefficient form.
• Damping of physical oscillations: BDF2 adds numerical damping to all modes; fine for diffusion, bad for structural vibration (use generalized-α or Newmark there).
• Order of the start-up: BDF2 started with one BDF1 step loses no global order (the single first-order step's local error is O(h²)).

Pitfalls

• Too large a first step on a nonlinear stiff problem: Newton converges to a spurious root. The code above with h = 10⁻³ from Robertson's initial state lands on y₂ < 0. Start with a tiny h and ramp up, as production codes do.
• Treating Newton failure as a fatal error instead of reducing h (production: halve h, refactor, retry).
• Finite-difference increments too small or not scaled by |y_j|: Jacobian is noise, Newton stalls.
• Using high-order BDF on problems with eigenvalues near the imaginary axis (lightly damped oscillations): BDF3+ can be unstable there.
• Trapezoidal rule on stiff problems with discontinuous inputs: ringing that never decays (R(∞) = −1); TR-BDF2 fixes this.
• Recomputing J and LU every Newton iteration when n is large: most of the runtime for no gain.

Variants and alternatives

In production

SUNDIALS CVODE (variable-order BDF 1-5 with Newton or Newton-Krylov) and IDA (DAEs), LSODE/LSODA (ODEPACK; LSODA switches between Adams and BDF), scipy solve_ivp(method='BDF' | 'Radau' | 'LSODA'), MATLAB ode15s (NDF) and ode23s, SPICE transient analysis (trapezoidal or Gear/BDF), DifferentialEquations.jl FBDF, QNDF, Rodas5.

Implementations

Sources

• Ernst Hairer, Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer, 1996. book
• C. William Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, 1971.
• Germund Dahlquist, A Special Stability Problem for Linear Multistep Methods, BIT 3, 1963. paper
• Alan C. Hindmarsh et al., SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers, ACM TOMS 31(3), 2005. paper · code
• Lawrence F. Shampine, Mark W. Reichelt, The MATLAB ODE Suite, SIAM Journal on Scientific Computing 18(1), 1997. paper