Skip to main content

Module bdf

numra_ode

Module bdf 

Source
Expand description

BDF / NDF: variable-order, variable-step backward differentiation formulas.

This is a Rust port of SciPy’s scipy.integrate.BDF (the canonical implementation of Shampine–Reichelt’s NDF method from MATLAB’s ode15s). The previous revision of this file was algorithmically incorrect — see “Why a rewrite” below.

§Mathematical formulation

The k-th order BDF method approximates M · ∑ⱼ αⱼ · y_{n+1-j} = h · β · f(t_{n+1}, y_{n+1}) where M is an optional (possibly singular) mass matrix. For ODEs M = I.

Rather than storing past y values and applying fixed-coefficient BDF formulas (which is not order-k accurate when the step size varies — this was the previous file’s fundamental error), we work with modified divided differences D, defined by: D⁰yₙ = yₙ Dʲyₙ = Dʲ⁻¹yₙ − Dʲ⁻¹y_{n-1} (j ≥ 1)

The differences are stored in an array D[0..k+2] whose rows are scaled so that the BDF formula holds as if the step size were constant. When h changes, the array is rescaled in-place via a small transformation matrix (change_d); this preserves order accuracy. This is the “fixed-leading-coefficient” / “quasi-constant step size” trick from Shampine & Reichelt (1997).

§NDF coefficients (Klopfenstein–Shampine)

For each order k we store a KAPPA[k] “boost” parameter; setting κ = 0 recovers pure BDF, while small negative κ gives the NDF method, which reduces the error constant by ~25–30% at the cost of slight stability tightening. Defaults follow Shampine–Reichelt.

§DAE support (index-1)

For singular M, the iteration matrix becomes M − c·J (instead of I − c·J), and the Newton residual carries M·ψ and M·d terms. Algebraic constraint rows (zero rows of M) reduce to f_alg(t, y) = 0, enforced at every step.

§Why a rewrite

The previous bdf.rs had three independent algorithmic bugs that manifested catastrophically (silent wrong answers, runaway step counts):

  1. Newton convergence test on the residual norm — when h·f was small in absolute terms, Newton declared convergence at iteration 0 with y_new = predictor = y_old, leaving the solution unchanged while time advanced. This is the source of the reported “y_final = y_initial at success=true” bug.
  2. Error estimator identically ≈ 0 — the formula referenced history[0], which was always equal to the current y (just pushed), so the estimate’s leading term cancelled.
  3. Fixed BDF coefficients applied to variable-h history — without the modified-divided-differences rescaling, the method silently loses order whenever h changes. The textbook coefficients [1, -4/3, 1/3, 2/3·h] only define BDF2 when y_n, y_{n-1} are spaced equally.

§References

  • L. F. Shampine and M. W. Reichelt, “The MATLAB ODE Suite”, SIAM J. Sci. Comput. 18(1), 1–22 (1997). [NDF coefficients, fixed-leading-coefficient BDF, error estimator.]
  • G. D. Byrne and A. C. Hindmarsh, “A Polyalgorithm for the Numerical Solution of ODEs”, ACM TOMS 1(1), 71–96 (1975).
  • SciPy scipy/integrate/_ivp/bdf.py (BSD-3 / Apache-2.0). Reference implementation that this file ports.

Author: Moussa Leblouba Date: 5 March 2026 Modified: 6 May 2026 (full rewrite based on SciPy’s NDF/MDD algorithm)

Structs§

Bdf
BDF / NDF solver for stiff ODEs and index-1 DAEs (orders 1-5).