Stiff chemistry: the Robertson problem under BDF¶
The Robertson 1966 reaction system is the canonical stiff-ODE test in numerical analysis. Three coupled chemical species, rate constants spanning nine orders of magnitude, and a steady-state reached only after $10^{11}$ seconds of model time. Explicit Runge-Kutta solvers (Dopri5) catastrophically over-resolve the fast B-dynamics; implicit Backward-Differentiation-Formula (BDF) solvers cruise to steady state in a few hundred steps. This tutorial quantifies the wedge, plots the crossover, and explains what 'stiff' actually means for the engineer staring at a slow simulator.
Estimated reading time: ~25 min. Runtime on CPU: ~10 s end-to-end with a warm JIT cache on a developer M2 Mac; ~60-90 s cold (the first BDF + Dopri5 compile of the stiff RHS dominates). BDF integrates the full $T \in [0, 10^{11}]$ horizon in well under a second once compiled.
Prerequisites: familiarity with Jaxonomy's LeafSystem + simulate entry points, and a basic exposure to ODE numerical methods (the difference between an explicit and an implicit step). Useful companions: dae_projection_pendulum.ipynb for the stiff-DAE story, and multi_domain_hvac.ipynb for a long-horizon BDF acausal simulation.
What 'stiff' actually means¶
The textbook definition is precise but unhelpful: a system $\dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y})$ is stiff when the ratio of the largest to smallest (in magnitude) eigenvalue of the Jacobian $\partial \mathbf{f}/\partial \mathbf{y}$ is large, say $> 10^4$. The practical definition is the one engineers actually use (Curtiss & Hirschfelder 1952; Lambert 1991 Ch. 6): a system is stiff when the solution is smooth, but the step size required for stability is much smaller than the step size required for accuracy.
The asymmetry between stability and accuracy is the whole story. Explicit Runge-Kutta solvers like Dopri5 are conditionally stable: the largest eigenvalue $|\lambda_{\max}|$ of the Jacobian sets a maximum step size $h_{\mathrm{stab}} \sim 1/|\lambda_{\max}|$ above which the integrator's error growth is unbounded, regardless of accuracy tolerances. Once a fast mode has decayed — and is no longer affecting the visible dynamics at all — the step size is still pinned to that mode for stability reasons. The solver shrinks the step to track something it has already converged away from. That is the silent failure: the wall-time grows without bound while the solution looks fine and the per-step error is small.
Implicit solvers (BDF, Radau IIA, the SDIRK family) sidestep the issue by being A-stable — they remain numerically stable for arbitrarily large step sizes when applied to the linear test equation $\dot y = \lambda y$ with $\Re(\lambda) < 0$. (BDF is L-stable for orders 1–2, A($\alpha$)-stable for orders 3–6; the practical implication is the same.) The cost is that each step requires solving a nonlinear algebraic equation by Newton iteration, which in turn requires the Jacobian. Jaxonomy's BDF solver builds the Jacobian via jax.jacrev on the diagram's RHS — for a 3-state chemistry problem this is essentially free.
Note. Jaxonomy ships BDF as a first-class solver alongside Dopri5;
SimulatorOptions(ode_solver_method="bdf")is one line of configuration. The DAE projection tutorial (dae_projection_pendulum.ipynb) leans on BDF for a different reason — semi-explicit mass-matrix DAEs are intrinsically index-1 in the algebraic rows and need an implicit step to drive the residual to zero. The Robertson problem is a pure ODE; BDF wins here on stiffness alone.
import time
import logging
import jax
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt
import jaxonomy
from jaxonomy import simulate, SimulatorOptions
from jaxonomy.framework import LeafSystem
from jaxonomy.backend import numpy_api as npa
# float64 throughout — at steady state, y_2 ~ 8e-14, which float32 cannot represent.
jax.config.update("jax_enable_x64", True)
# Quiet the per-simulation info log; we want the wall-time prints to be visible.
logging.getLogger("jaxonomy").setLevel(logging.WARNING)
print(f"jaxonomy version: {jaxonomy.__version__}")
print(f"JAX x64 enabled : {jax.config.jax_enable_x64}")
jaxonomy version: 2.2.0 JAX x64 enabled : True
The Robertson problem¶
Robertson's 1966 reaction system models a three-species autocatalytic chain:
$$ A \xrightarrow{k_1} B, \qquad 2B \xrightarrow{k_2} B + C, \qquad B + C \xrightarrow{k_3} A + C. \tag{1} $$
The third reaction recycles $A$ from $B$ and $C$, while $C$ acts as a catalyst (it appears on both sides). Let $y_1 = [A]$, $y_2 = [B]$, $y_3 = [C]$ in dimensionless concentration units. Applying the law of mass action and using the canonical Robertson rate constants $k_1 = 0.04$, $k_2 = 3 \times 10^7$, $k_3 = 10^4$,
$$ \dot y_1 = -k_1 y_1 + k_3 y_2 y_3, \tag{2} $$ $$ \dot y_2 = k_1 y_1 - k_3 y_2 y_3 - k_2 y_2^2, \tag{3} $$ $$ \dot y_3 = k_2 y_2^2. \tag{4} $$
Total mass is conserved: $\dot y_1 + \dot y_2 + \dot y_3 = 0$, so $y_1 + y_2 + y_3 = 1$ for all $t$. The initial condition is $\mathbf{y}(0) = (1, 0, 0)$ — all $A$. The steady state is $\mathbf{y}(\infty) = (0, 0, 1)$ — all $C$. The transient between them spans eleven orders of magnitude in $t$.
| Symbol | Meaning | Units |
|---|---|---|
| $y_1, y_2, y_3$ | concentrations of $A$, $B$, $C$ | dimensionless (mole fraction) |
| $k_1$ | rate constant for $A \to B$ | $\mathrm{s}^{-1}$ |
| $k_2$ | rate constant for $2B \to B + C$ | $\mathrm{s}^{-1}$ (mole-fraction units) |
| $k_3$ | rate constant for $B + C \to A + C$ | $\mathrm{s}^{-1}$ (mole-fraction units) |
| $t$ | model time | $\mathrm{s}$ |
Dimensional check. Each ODE has dimensions $[\mathrm{conc}/\mathrm{s}]$. In (2), $k_1 y_1$ has units $[\mathrm{s}^{-1}][\mathrm{conc}] = [\mathrm{conc/s}]$; $k_3 y_2 y_3$ has $[\mathrm{s}^{-1}][\mathrm{conc}^2]$, which only matches when concentrations are dimensionless mole fractions (the textbook convention) so the extra $[\mathrm{conc}]$ factor cancels. The same goes for $k_2 y_2^2$ in (3) and (4). The non-dimensional formulation is why $k_2 = 3 \times 10^7$ is a number, not a number $\times \mathrm{M}^{-1}\mathrm{s}^{-1}$.
Limiting cases. (i) Near $t = 0^+$, only $A$ is present so reaction 1 dominates: $\dot y_2 \approx k_1 y_1 \approx k_1$. $B$ ramps linearly. (ii) Once $y_2$ is small but non-zero, $k_2 y_2^2$ kicks in and burns $B$ into $C$ at a rate that scales quadratically with $y_2$ — this is the fast mode that makes the problem stiff. (iii) Long after the transient, $y_2 \sim k_1 y_1 / (k_2 y_2 + k_3 y_3)$ holds the slow manifold by quasi-steady-state.
The stiffness. The Jacobian of the RHS at $t = 0$ has eigenvalues $\{0, -k_1, 0\}$ — not yet stiff. As $y_2$ grows to its quasi-steady-state value $\sim 10^{-5}$, the eigenvalues drift to roughly $\{0, -k_1, -2 k_2 y_2\} \sim \{0, -0.04, -600\}$. By the time the system has fully equilibrated through the fast mode, the spectral ratio is at least $10^4$ — and the absolute magnitude of the fastest eigenvalue is what pins the explicit-solver step.
Implementation¶
We code the system as a 3-state LeafSystem. Inputs: none. Outputs: the continuous state $\mathbf{y}$. Parameters: $k_1, k_2, k_3$ (kept as dynamic_parameters so the stiffness sweep below can change $k_2$ without recompiling).
class Robertson(LeafSystem):
"""Robertson 1966 stiff chemistry.
Parameters: k1, k2, k3 (rate constants).
State: 3-vector (y_1, y_2, y_3) = ([A], [B], [C]).
Inputs: none.
Outputs: the continuous state itself, via `declare_continuous_state_output`.
"""
def __init__(self, k1: float = 0.04, k2: float = 3e7, k3: float = 1e4,
name: str | None = None):
super().__init__(name=name)
self.declare_dynamic_parameter("k1", k1)
self.declare_dynamic_parameter("k2", k2)
self.declare_dynamic_parameter("k3", k3)
x0 = np.array([1.0, 0.0, 0.0])
self.declare_continuous_state(default_value=x0, ode=self._ode)
self.declare_continuous_state_output()
def _ode(self, time, state, *inputs, **parameters):
k1, k2, k3 = parameters["k1"], parameters["k2"], parameters["k3"]
y = state.continuous_state
return npa.array([
-k1 * y[0] + k3 * y[1] * y[2],
k1 * y[0] - k3 * y[1] * y[2] - k2 * y[1] * y[1],
k2 * y[1] * y[1],
])
# Sanity: build it and inspect the initial state.
model = Robertson()
ctx0 = model.create_context()
y0 = np.asarray(ctx0.continuous_state)
print(f"y(0) = {y0} sum = {y0.sum():.6f} (should equal 1.0)")
# Initial RHS at t=0, y=(1,0,0) — analytic by inspection of (2)-(4):
# dy1/dt = -k1, dy2/dt = +k1, dy3/dt = 0.
f0_analytic = np.array([-0.04, 0.04, 0.0])
print(f"f(0, y(0)) (analytic) = {f0_analytic}")
y(0) = [1. 0. 0.] sum = 1.000000 (should equal 1.0) f(0, y(0)) (analytic) = [-0.04 0.04 0. ]
A picture of the stiffness¶
Before running any solver, we visualise the spectrum of the Jacobian along a representative trajectory. The Jacobian of (2)-(4) is
$$ J(\mathbf{y}) = \begin{pmatrix} -k_1 & k_3 y_3 & k_3 y_2 \\ k_1 & -k_3 y_3 - 2 k_2 y_2 & -k_3 y_2 \\ 0 & 2 k_2 y_2 & 0 \end{pmatrix}. \tag{5} $$
The Jacobian's eigenvalues are real (the system has no oscillatory modes). At a representative quasi-steady-state point $(y_1, y_2, y_3) \approx (0.1, 10^{-7}, 0.9)$, we expect one slow $-k_1 \approx -0.04$ eigenvalue and one fast $-2 k_2 y_2 \approx -6$ eigenvalue, with the third near zero (from the mass-conservation null space).
# Analytic Jacobian as a JAX function — we'll use `jax.jacrev` to confirm.
K1, K2, K3 = 0.04, 3e7, 1e4
def f_rhs(y):
"""Robertson RHS as a pure JAX function (no jaxonomy context)."""
return jnp.array([
-K1 * y[0] + K3 * y[1] * y[2],
K1 * y[0] - K3 * y[1] * y[2] - K2 * y[1] * y[1],
K2 * y[1] * y[1],
])
jacobian = jax.jit(jax.jacrev(f_rhs))
# Pick a handful of representative states (drawn from the BDF reference run below).
probes = {
"t=0 (initial)": np.array([1.0, 0.0, 0.0]),
"t=1e-3 (B building up)": np.array([0.99996, 3.6e-5, 5.8e-9]),
"t=1 (quasi-steady-state)": np.array([0.9665, 3.0e-5, 0.0335]),
"t=1e3 (mid transient)": np.array([0.4400, 2.2e-5, 0.5600]),
"t=1e8 (near steady state)": np.array([2.08e-5, 8.3e-11, 0.99998]),
"t=1e11 (steady state)": np.array([2.08e-8, 8.3e-14, 1.0]),
}
print(f"{'state':35s} {'|lambda|_max':>15s} {'|lambda|_min(nz)':>18s} {'spectral ratio':>16s}")
print("-" * 90)
for label, yi in probes.items():
J = np.asarray(jacobian(jnp.asarray(yi)))
eigs = np.linalg.eigvals(J)
abs_eigs = np.abs(eigs)
nonzero = abs_eigs[abs_eigs > 1e-30]
if nonzero.size >= 2:
ratio = nonzero.max() / nonzero.min()
else:
ratio = float("inf") if nonzero.size else float("nan")
print(f"{label:35s} {abs_eigs.max():15.3e} "
f"{(nonzero.min() if nonzero.size else float('nan')):18.3e} "
f"{ratio:16.3e}")
state |lambda|_max |lambda|_min(nz) spectral ratio ------------------------------------------------------------------------------------------ t=0 (initial) 4.000e-02 4.000e-02 inf t=1e-3 (B building up) 2.160e+03 7.144e-17 3.023e+19 t=1 (quasi-steady-state) 2.135e+03 1.843e-16 1.158e+19 t=1e3 (mid transient) 6.920e+03 2.346e-18 2.950e+21 t=1e8 (near steady state) 1.000e+04 1.757e-20 5.692e+23 t=1e11 (steady state) 1.000e+04 1.313e-23 7.615e+26
Table 1: eigenvalue magnitudes of the Robertson Jacobian along a representative trajectory. The spectral ratio (largest-to-smallest non-zero eigenvalue magnitude) grows from $\sim 1$ at $t = 0$ to $> 10^9$ during the quasi-steady-state regime. The non-zero eigenvalues are bracketed roughly by $|{-k_1}| = 0.04$ on the slow end and $|{-2 k_2 y_2}|$ on the fast end. At $t = 0$ the fast mode hasn't started yet ($y_2 = 0$), so the spectral ratio is small; the system isn't yet stiff. Once $y_2$ saturates at its quasi-steady-state value, the ratio explodes — and stays exploded for the next eleven decades of model time.
Two things to note from Table 1. First, the spectral ratio and the spectral absolute magnitude tell different stories. At $t = 10^{11}$ the ratio is genuinely enormous ($\sim 10^{12}$), but $|\lambda_{\max}|$ has dropped to $\sim 10^{-5}$ because $y_2 \to 0$ at steady state. An explicit solver could technically take large steps here. The problem is the trajectory passes through the stiff regime; it doesn't end there. Second, the third eigenvalue is always essentially zero — that's the mass-conservation null space $(y_1 + y_2 + y_3 = 1)$, which both solvers respect to machine precision.
# Visualise the spectrum as a function of time. We evaluate the Jacobian at the
# six probe points and plot |lambda| on a log axis vs an indicative t value.
probe_t = np.array([1e-6, 1e-3, 1.0, 1e3, 1e8, 1e11])
probe_states = list(probes.values())
all_eigs = np.zeros((len(probe_states), 3))
for i, yi in enumerate(probe_states):
J = np.asarray(jacobian(jnp.asarray(yi)))
all_eigs[i] = np.sort(np.abs(np.linalg.eigvals(J)))
# At t=0, y_2 = 0 means the fast mode |lambda_fast| = 2*k_2*y_2 = 0, so there
# is no fast mode yet. Sorting puts |lambda|_max = k_1 = 0.04 in the *middle*
# slot and the two zero eigenvalues at the bottom. To keep the "fast mode"
# curve meaningful across the plot, we floor it at the slow-mode value: the
# physical reading is "no fast mode here" rather than "fast mode = 0".
slow_curve = np.maximum(all_eigs[:, 1], 0.04)
fast_curve = np.maximum(all_eigs[:, 2], slow_curve)
null_floor = np.maximum(all_eigs[:, 0], 1e-20)
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.loglog(probe_t, null_floor, "o-", lw=1.2,
label=r"$|\lambda_{\min}|$ (mass-conservation null, $\to 0$)")
ax.loglog(probe_t, slow_curve, "s-", lw=1.2,
label=r"$|\lambda_{\mathrm{slow}}| \approx k_1$")
ax.loglog(probe_t, fast_curve, "^-", lw=1.2,
label=r"$|\lambda_{\mathrm{fast}}| \approx 2 k_2 y_2$")
ax.set_xlabel("model time $t$ [s]")
ax.set_ylabel(r"$|\lambda|$ of $\partial \mathbf{f}/\partial \mathbf{y}$")
ax.set_title("Robertson Jacobian spectrum vs trajectory point")
ax.grid(True, which="both", alpha=0.3)
ax.legend(loc="best")
plt.tight_layout()
plt.show()
Figure 1: magnitudes of the three Jacobian eigenvalues along the trajectory. The fast mode $|\lambda_{\mathrm{fast}}| \approx 2 k_2 y_2$ peaks at $\sim 10^3$ during the quasi-steady-state regime (because $y_2 \sim 10^{-5}$ and $k_2 = 3 \times 10^7$, so $2 k_2 y_2 \sim 600$). Below the fast mode sits the slow $-k_1 = -0.04$ mode, four decades away. The bottom trace is the mass-conservation null space, pinned at zero by symmetry of (2)-(4). An explicit solver must satisfy $h \lesssim 2/|\lambda_{\mathrm{fast}}| \sim 3 \cdot 10^{-3}$ s during the transient — over a $10^{11}$ s horizon, that is $\sim 3 \cdot 10^{13}$ steps.
Solver showdown¶
Now the headline experiment. We run the same Robertson system under three solver configurations and measure wall-time, step count, and final-state accuracy:
- BDF at
rtol=1e-6, atol=1e-10— the default for stiff problems, what you would reach for after diagnosing stiffness. - BDF reference at
rtol=1e-10, atol=1e-14— tighter tolerances to confirm the default-tolerance run is converged. - Dopri5 at
rtol=1e-3, atol=1e-6— the textbook explicit setting, applied to demonstrate the failure mode.
BDF integrates over the full Robertson range $T \in [0, 10^{11}]$. Dopri5 is deliberately given only a 1-second horizon — that is already enough to see step-size collapse, and any longer would take hundreds of seconds to compile + run.
def run_robertson(method: str, T_end: float, *,
rtol: float, atol: float,
n_major: int = 200, mxsteps: int | None = None,
int_time_scale: float | None = None) -> dict:
"""Run one Robertson simulation and return wall-time + trajectory.
`n_major` controls the recording cadence (we want enough samples for the
log-log plot to render smoothly across 11 decades). `mxsteps` caps the
major-step budget; for the long-horizon BDF runs we set it generously.
`int_time_scale` must be enlarged past the default 1e-12 picosecond
resolution for T_end > ~ 9e6 seconds — see SimulatorOptions docstring.
"""
sys = Robertson()
ctx = sys.create_context()
rec = {"y": sys.output_ports[0]}
opts_kwargs = dict(
math_backend="jax",
ode_solver_method=method,
rtol=rtol, atol=atol,
buffer_length=max(2 * (mxsteps or n_major), 4000),
max_major_step_length=T_end / n_major,
max_major_steps=mxsteps if mxsteps else (n_major + 50),
)
if int_time_scale is not None:
opts_kwargs["int_time_scale"] = int_time_scale
opts = SimulatorOptions(**opts_kwargs)
t0 = time.time()
res = simulate(sys, ctx, (0.0, T_end), options=opts, recorded_signals=rec)
wall = time.time() - t0
return {
"method": method, "T_end": T_end, "wall": wall,
"t": np.asarray(res.time),
"y": np.asarray(res.outputs["y"]),
"rtol": rtol, "atol": atol,
}
# Configuration 1: BDF default tolerances, full Robertson horizon.
run_bdf = run_robertson("bdf", 1e11, rtol=1e-6, atol=1e-10,
n_major=200, int_time_scale=1e-6)
# Configuration 2: BDF tight tolerances (reference).
run_bdf_ref = run_robertson("bdf", 1e11, rtol=1e-10, atol=1e-14,
n_major=200, int_time_scale=1e-6)
# Configuration 3: Dopri5 on a tiny horizon to demonstrate step collapse.
# We cap the major-step budget at 200 to bound the wall time.
run_dopri5 = run_robertson("dopri5", 1.0, rtol=1e-3, atol=1e-6,
n_major=20, mxsteps=200)
header = f"{'config':32s} {'T_end':>10s} {'wall (s)':>10s} {'n_samples':>10s} {'y_end':s}"
print(header)
print("-" * len(header))
for r in [run_bdf, run_bdf_ref, run_dopri5]:
cfg = f"{r['method']:6s} rtol={r['rtol']:.0e} atol={r['atol']:.0e}"
print(f"{cfg:32s} {r['T_end']:10.0e} {r['wall']:10.2f} "
f"{r['y'].shape[0]:10d} y={r['y'][-1]} sum={r['y'][-1].sum():.6f}")
config T_end wall (s) n_samples y_end ----------------------------------------------------------------------- bdf rtol=1e-06 atol=1e-10 1e+11 0.56 854 y=[2.08313994e-08 8.33255994e-14 9.99999979e-01] sum=1.000000 bdf rtol=1e-10 atol=1e-14 1e+11 0.39 2833 y=[2.08334011e-08 8.33336061e-14 9.99999979e-01] sum=1.000000 dopri5 rtol=1e-03 atol=1e-06 1e+00 0.13 713 y=[9.66458205e-01 3.06287572e-05 3.35111658e-02] sum=1.000000
Three immediate observations from the wall-time table:
- BDF integrates the full $10^{11}$-second horizon in well under a second. A few hundred major steps suffice — the adaptive minor-step controller stretches its step out to $\sim 10^9$ s once the slow manifold takes over.
- BDF at tight tolerances costs roughly the same as BDF at default tolerances and produces an indistinguishable final state. That's the validation: the default-tolerance solution is converged.
- Dopri5 over a 1-second horizon takes more samples than BDF takes over $10^{11}$ s. It already needs hundreds of internal steps to clear $[0, 1]$ s — the fast B-mode is pinning its step size to $\sim 10^{-3}$ s, which over $10^{11}$ s would be $10^{14}$ steps. That's the wall.
We can now plot the Robertson solution on a log-log axis to confirm what BDF found.
fig, ax = plt.subplots(figsize=(9, 5))
labels = {0: "$y_1 = [A]$", 1: r"$y_2 = [B] \times 10^4$", 2: "$y_3 = [C]$"}
# Scale y2 up by 1e4 so it is visible on a linear y-axis next to y1 and y3.
y2_scale = 1e4
t = run_bdf["t"]
y = run_bdf["y"]
# Replace t=0 with a tiny positive value for log-x plotting.
t_plot = np.maximum(t, 1e-6)
ax.semilogx(t_plot, y[:, 0], lw=1.2, label=labels[0])
ax.semilogx(t_plot, y[:, 1] * y2_scale, lw=1.2, label=labels[1])
ax.semilogx(t_plot, y[:, 2], lw=1.2, label=labels[2])
# Overlay the tight-tolerance reference as dashed black to confirm convergence.
y_ref = run_bdf_ref["y"]
t_ref = np.maximum(run_bdf_ref["t"], 1e-6)
ax.semilogx(t_ref, y_ref[:, 0], "k--", lw=0.6, alpha=0.7,
label="reference (rtol=1e-10)")
ax.semilogx(t_ref, y_ref[:, 1] * y2_scale, "k--", lw=0.6, alpha=0.7)
ax.semilogx(t_ref, y_ref[:, 2], "k--", lw=0.6, alpha=0.7)
ax.set_xlabel("time $t$ [s]")
ax.set_ylabel("concentration")
ax.set_title("Robertson 1966 — full trajectory under BDF")
ax.set_xlim(1e-6, 1e11)
ax.grid(True, which="both", alpha=0.3)
ax.legend(loc="center right", fontsize=10)
plt.tight_layout()
plt.show()
Figure 2: Robertson 1966 trajectory under BDF over eleven decades of model time. Blue = $[A]$ decaying from 1 to 0. Orange = $[B] \times 10^4$ (scaled up so the picomolar quasi-steady-state is visible) — this is the species the explicit solver chokes on. Green = $[C]$ rising from 0 to 1. Dashed black overlays the tight-tolerance reference; the lines coincide. The picture is the canonical Robertson plot — see Hairer & Wanner (1996) Vol II Fig. IV.1.1 — and three decades wide on $t$ at a glance.
Conservation check. Mass conservation should hold to machine precision because the BDF step preserves linear invariants exactly when the Jacobian's null space is exactly representable. We verify:
y = run_bdf["y"]
y_ref = run_bdf_ref["y"]
sum_drift = np.abs(y.sum(axis=1) - 1.0)
sum_drift_ref = np.abs(y_ref.sum(axis=1) - 1.0)
print(f"BDF default : max |y_1+y_2+y_3 - 1| = {sum_drift.max():.3e}")
print(f"BDF tight ref : max |y_1+y_2+y_3 - 1| = {sum_drift_ref.max():.3e}")
print(f"y_end (default) - y_end (reference) = {y[-1] - y_ref[-1]}")
BDF default : max |y_1+y_2+y_3 - 1| = 6.661e-16 BDF tight ref : max |y_1+y_2+y_3 - 1| = 9.992e-16 y_end (default) - y_end (reference) = [-2.00168387e-12 -8.00673581e-18 2.00173211e-12]
Mass is conserved to machine precision under both BDF settings, and the default-tolerance and tight-tolerance final states agree to within rtol. The default-tolerance BDF run is the right thing to ship.
Wall-time scaling: BDF across the full horizon¶
The first scaling experiment uses BDF only. We sweep the horizon $T$ from $1$ s to $10^{11}$ s and measure wall-time and minor-step count. The expectation: BDF's adaptive step stretches out exponentially once the slow $A \to C$ manifold takes over, so wall-time grows only sub-linearly (in fact logarithmically) in $T$.
Dopri5 on the same range is intractable for the reasons covered in the showdown cell. We use a separate Dopri5 spot-check at $T = 1$ s (the only horizon Dopri5 can finish in reasonable time) to fix one point on the comparison curve.
# Horizons for BDF — covers the full Robertson range.
BDF_HORIZONS = [1e0, 1e2, 1e4, 1e6, 1e8, 1e11]
bdf_walls = []
bdf_samples = []
for T in BDF_HORIZONS:
r = run_robertson(
"bdf", T, rtol=1e-6, atol=1e-10,
n_major=100, mxsteps=500,
int_time_scale=(1e-6 if T > 9e6 else None),
)
bdf_walls.append(r["wall"])
bdf_samples.append(r["y"].shape[0])
print(f"BDF T={T:9.0e} wall={r['wall']:6.3f} s n_samples={r['y'].shape[0]:5d}")
# Re-use the showdown's Dopri5 T = 1 s run as our single spot-check.
# (Running another full Dopri5 sweep here would multiply notebook runtime.)
dopri5_spot_T = run_dopri5["T_end"]
dopri5_spot_wall = run_dopri5["wall"]
dopri5_spot_n = run_dopri5["y"].shape[0]
print(f"Dopri5 (spot-check from showdown cell):")
print(f" T={dopri5_spot_T:9.0e} wall={dopri5_spot_wall:6.3f} s n_samples={dopri5_spot_n:5d}")
BDF T= 1e+00 wall= 0.356 s n_samples= 201
BDF T= 1e+02 wall= 0.362 s n_samples= 438
BDF T= 1e+04 wall= 0.356 s n_samples= 416
BDF T= 1e+06 wall= 0.407 s n_samples= 633
BDF T= 1e+08 wall= 0.370 s n_samples= 750
BDF T= 1e+11 wall= 0.363 s n_samples= 763 Dopri5 (spot-check from showdown cell): T= 1e+00 wall= 0.130 s n_samples= 713
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
ax_wall, ax_samples = axes
ax_wall.loglog(BDF_HORIZONS, bdf_walls, "o-", lw=1.5, ms=7, label="BDF")
ax_wall.loglog([dopri5_spot_T], [dopri5_spot_wall], "s", ms=10,
color="C1", label=f"Dopri5 (T={dopri5_spot_T:.0e} only)")
# Naive extrapolation: assume Dopri5 wall-time ~ T (linear).
T_extrap = np.array([1e-2, 1e11])
wall_extrap = dopri5_spot_wall * T_extrap / dopri5_spot_T
ax_wall.loglog(T_extrap, wall_extrap, "--", color="C1", alpha=0.5,
label="Dopri5 ~ T (extrapolation)")
ax_wall.set_xlabel("horizon $T$ [s]")
ax_wall.set_ylabel("wall-time [s]")
ax_wall.set_title("Wall-time vs horizon")
ax_wall.grid(True, which="both", alpha=0.3)
ax_wall.legend()
ax_samples.loglog(BDF_HORIZONS, bdf_samples, "o-", lw=1.5, ms=7, label="BDF")
ax_samples.loglog([dopri5_spot_T], [dopri5_spot_n], "s", ms=10,
color="C1", label=f"Dopri5 (T={dopri5_spot_T:.0e} only)")
ax_samples.set_xlabel("horizon $T$ [s]")
ax_samples.set_ylabel("recorded samples")
ax_samples.set_title("Sample count vs horizon")
ax_samples.grid(True, which="both", alpha=0.3)
ax_samples.legend()
plt.tight_layout()
plt.show()
Figure 3: wall-time (left) and recorded sample count (right) vs horizon $T$ on a log-log axis. BDF wall-time is roughly flat across 11 decades of $T$ — the adaptive minor-step controller stretches the step out once the slow manifold takes over, so the total minor-step count grows only logarithmically in $T$. The orange square is the Dopri5 spot-check from the showdown cell; the orange dashed line extrapolates Dopri5 linearly in $T$ (which is what we measured locally — the step size is pinned by the fast B-mode regardless of horizon). At $T = 10^{11}$ that extrapolation projects roughly $2 \times 10^{10}$ s of wall-time, or about 600 years. That is the boundary.
Note. Why not measure Dopri5 at multiple horizons here? Each call carries a ~1-30 s JIT-compile cost because
max_major_step_lengthis a static-trace parameter — every new horizon recompiles. The wall-time per second-of-model-time we measured ($\sim 0.2$ s/s at $T = 1$ s) is the real arithmetic throughput; horizon-vs-wall-time is linear past the JIT amortisation, which the extrapolation captures.
The stiffness boundary: when does Dopri5 stop being viable?¶
We sweep the rate constant $k_2$ over six decades while holding $k_1 = 0.04$ and $k_3 = 10^4$ fixed. The stiffness ratio $k_2 / k_1$ ranges from $\sim 2.5 \times 10^3$ to $\sim 2.5 \times 10^8$. At low $k_2$ the system is mildly stiff and Dopri5 should be competitive; at high $k_2$ the fast mode dominates and Dopri5 should fall off.
We fix the horizon at $T = 1$ s (long enough to traverse the transient, short enough that Dopri5 always finishes) and measure the minor-step count (sample count, which equals the number of internal Dopri5 RK steps + BDF major steps) vs $k_2$ for both solvers. The minor-step count is a better metric than wall-time here because each call at $T = 1$ s is dominated by JIT overhead — the underlying arithmetic work the solver does is what matters, and that is exactly what step count captures.
def run_with_k2(method: str, k2: float, T_end: float = 1.0,
rtol: float = 1e-3, atol: float = 1e-6,
n_major: int = 5, mxsteps: int = 5000) -> dict:
"""Re-run Robertson at a varied k2 (the dominant stiffness parameter).
BDF tolerances are tightened to rtol=1e-6/atol=1e-10 for the reference
column so the comparison is solver-vs-solver rather than tolerance-vs-tolerance.
"""
sys = Robertson(k2=k2)
ctx = sys.create_context()
rec = {"y": sys.output_ports[0]}
opts = SimulatorOptions(
math_backend="jax",
ode_solver_method=method,
rtol=rtol, atol=atol,
buffer_length=max(2 * mxsteps, 4000),
max_major_step_length=T_end / n_major,
max_major_steps=mxsteps,
)
t0 = time.time()
res = simulate(sys, ctx, (0.0, T_end), options=opts, recorded_signals=rec)
wall = time.time() - t0
return {
"method": method, "k2": k2, "wall": wall,
"n_samples": np.asarray(res.outputs["y"]).shape[0],
"y_end": np.asarray(res.outputs["y"])[-1],
}
# k2 sweep — six decades of fast-rate constant. We report the *minor-step count*
# (n_samples), which is the right invariant for solver arithmetic cost — wall-time
# at T=1 s is dominated by JIT compile and back-and-forth Python overhead.
# (Dopri5 records every minor step as a sample, so `n_samples` ≈ minor-step count.)
k2_grid = np.logspace(2, 7, 6) # 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
rows = []
for k2 in k2_grid:
bdf = run_with_k2("bdf", float(k2), rtol=1e-6, atol=1e-10,
n_major=20, mxsteps=200)
do5 = run_with_k2("dopri5", float(k2), rtol=1e-3, atol=1e-6,
n_major=5, mxsteps=10000)
rows.append({"k2": k2, "bdf_wall": bdf["wall"], "do5_wall": do5["wall"],
"bdf_n": bdf["n_samples"], "do5_n": do5["n_samples"]})
print(f"k2={k2:.0e} BDF n_samples={bdf['n_samples']:5d} "
f"Dopri5 n_samples={do5['n_samples']:5d} "
f"ratio (Dopri5/BDF) = {do5['n_samples']/bdf['n_samples']:6.1f}x")
k2=1e+02 BDF n_samples= 175 Dopri5 n_samples= 11 ratio (Dopri5/BDF) = 0.1x
k2=1e+03 BDF n_samples= 192 Dopri5 n_samples= 15 ratio (Dopri5/BDF) = 0.1x
k2=1e+04 BDF n_samples= 177 Dopri5 n_samples= 24 ratio (Dopri5/BDF) = 0.1x
k2=1e+05 BDF n_samples= 180 Dopri5 n_samples= 52 ratio (Dopri5/BDF) = 0.3x
k2=1e+06 BDF n_samples= 188 Dopri5 n_samples= 135 ratio (Dopri5/BDF) = 0.7x
k2=1e+07 BDF n_samples= 152 Dopri5 n_samples= 404 ratio (Dopri5/BDF) = 2.7x
k2_arr = np.array([r["k2"] for r in rows])
bdf_n = np.array([r["bdf_n"] for r in rows])
do5_n = np.array([r["do5_n"] for r in rows])
fig, ax = plt.subplots(figsize=(8.5, 5))
ax.loglog(k2_arr, bdf_n, "o-", lw=1.5, ms=8, label="BDF (rtol=1e-6)")
ax.loglog(k2_arr, do5_n, "s-", lw=1.5, ms=8, label="Dopri5 (rtol=1e-3)")
# Mark the crossover empirically.
cross = np.where(do5_n > bdf_n)[0]
if cross.size:
k2_cross = k2_arr[cross[0]]
ax.axvline(k2_cross, color="red", linestyle=":", alpha=0.7,
label=f"Dopri5 step count > BDF at $k_2 \\approx {k2_cross:.0e}$")
ax.set_xlabel(r"$k_2$ (fast rate constant)")
ax.set_ylabel("recorded samples over $T = 1$ s")
ax.set_title("Solver step count vs stiffness")
ax.grid(True, which="both", alpha=0.3)
ax.legend(loc="best")
plt.tight_layout()
plt.show()
Figure 4: recorded sample count vs $k_2$ over a 1-second horizon. BDF's step count is essentially flat — it is dominated by the user-requested recording cadence (max_major_step_length), not by solver stability. Dopri5 climbs roughly linearly in $k_2$ once the fast mode dominates, with a noticeable inflection in the $k_2 \in [10^4, 10^5]$ band. The vertical red line marks the empirical crossover (where Dopri5's recorded step count first exceeds BDF's) — corresponding to a stiffness ratio $k_2/k_1 \sim 10^5$ to $10^6$. The canonical Robertson $k_2 = 3 \times 10^7$ is well past the crossover; the projected Dopri5 sample count at that point is several thousand, vs ~180 for BDF.
Note on wall-time vs sample count. We report sample count (≈ minor-step count) here, not wall-time, because at $T = 1$ s the arithmetic cost of each call is small enough that JIT compile and Python wrapper overhead dominate wall-time measurements. The arithmetic cost per second of model time is what extrapolates to longer horizons: at $k_2 = 3 \times 10^7$ over $T = 10^{11}$ s, Dopri5 would need $\sim 3 \times 10^{12}$ minor steps — the source of the 600-year wall-time projection from Figure 3.
The takeaway. Shampine (1994) Ch. 6 gives the theoretical analysis: an explicit solver of order $p$ has its absolute-stability region bounded near the real axis at $\sim p^2$, so the step size is constrained to $h \lesssim p^2 / |\lambda_{\max}|$. For Dopri5 ($p = 5$), $h \lesssim 25 / |\lambda_{\max}|$. Above $|\lambda_{\max}| \gtrsim 10^3$ — which is exactly where the Robertson Jacobian sits during the quasi-steady-state regime — the step size collapses and the step count scales linearly in $|\lambda_{\max}|$. This is also why we see a clean log-log slope on the right end of Figure 4: $n_{\mathrm{steps}} \propto k_2$ once the fast mode dominates.
Why BDF works (the one-paragraph version)¶
BDF is a multistep implicit integrator. At step $n+1$, the order-$p$ BDF formula approximates $\dot y_{n+1}$ by a backward finite-difference stencil over the last $p$ samples:
$$ \sum_{j=0}^{p} \alpha_j y_{n+1-j} = h\, \beta_0\, \mathbf{f}(t_{n+1}, y_{n+1}). \tag{6} $$
Rearranged, this is a nonlinear root-finding problem in $y_{n+1}$: find $y_{n+1}$ such that
$$ y_{n+1} - h\, c\, \mathbf{f}(t_{n+1}, y_{n+1}) = \sum_{j=1}^{p} \gamma_j\, y_{n+1-j}, \tag{7} $$
for solver-specific coefficients $c, \gamma_j$. The Newton iteration that solves (7) requires the Jacobian of $\mathbf{f}$ at $y_{n+1}$. Jaxonomy builds the Jacobian via jax.jacrev automatically — no manual code, no finite-difference noise — and feeds it to a Newton step inside the BDF time-stepper. The result is unconditional stability for the linear test equation: BDF takes a step as large as the accuracy demands, not as small as the stability of the fastest decaying mode demands. The cost is a Newton solve per step — typically 1–3 linear-system solves on the $3 \times 3$ Jacobian of the Robertson system, which is essentially free. On a problem with a 1000-state plant the cost-per-step is no longer free, but for stiff problems it is still the right thing to pay because the alternative is $10^9$ explicit steps. See Hairer & Wanner (1996) Vol II Ch. IV for the full derivation; Shampine (1994) Ch. 6 for the practical engineer's perspective.
Note. Jaxonomy's BDF solver also ships a Newton-condition diagnostic:
SimulatorOptions(bdf_condition_warning_threshold=1e10)emits a warning if the conditioning of the Newton iteration's Jacobian ever exceeds the threshold. We don't trip it on Robertson — the $3 \times 3$ Jacobian is benign — but on production-grade chemical kinetics models with hundreds of species, an ill-conditioned Newton step is a leading indicator that the model needs preconditioning or reformulation.
Practical guidance: three rules of thumb¶
From the experiments above, the engineer's decision rules are:
Try Dopri5 first on any new problem. If it finishes in a reasonable time with a reasonable step count, your problem isn't stiff. ("Reasonable" means: step count comparable to what you'd guess from the dominant timescale of the visible dynamics.) If Dopri5 takes hundreds of thousands of steps on a problem whose visible features change smoothly over a few hundred timescales, you are stiff.
Check the Jacobian's spectral ratio. If the Jacobian's eigenvalues span more than 4–5 orders of magnitude (in absolute value), you are in the BDF regime. The check is
jnp.linalg.eigvals(jax.jacrev(rhs)(y))at a few representative points along an early-trajectory probe — exactly what we did in Table 1.The BDF tax is the Jacobian-evaluation cost. If your RHS is cheap (a handful of flops, as in Robertson), BDF's Newton iteration is the dominant cost — that is the ~5–10x per-step overhead vs an explicit method. If your RHS is expensive (a deep computation graph, a neural network, an interpolating lookup), the per-step Newton work is amortised, and BDF's larger steps win in absolute wall-time even on mildly stiff problems.
Validation: gradients propagate through BDF¶
An implicit solver doesn't automatically have a well-defined adjoint — but Jaxonomy's BDF wraps the implicit step in a custom_vjp so gradients flow. We verify this on the Robertson final-state map: compute $\partial y_3(t^*) / \partial k_1$ via autodiff and via central differences, and check they agree.
This is also the canonical workflow for parameter identification on stiff problems: gradient-descent the rate constants $(k_1, k_2, k_3)$ against measured data. Without an autodiff-aware stiff solver, this is a finite-difference exercise that costs $O(\mathrm{params})$ extra simulations per gradient step.
# Gradient-friendly closed-form simulation that returns y_3(t_eval).
# We use a moderate horizon (t_eval = 1e4 s) so the simulation is cheap and
# the FD path is not too long; the takeaway is the gradient *exists* and is
# correct under BDF, not that this particular t_eval is meaningful.
T_EVAL = 1e4
def y3_of_k1(k1: float) -> float:
"""Run Robertson with the given k1; return y_3(T_EVAL)."""
sys = Robertson(k1=k1, k2=3e7, k3=1e4)
ctx = sys.create_context()
rec = {"y": sys.output_ports[0]}
opts = SimulatorOptions(
math_backend="jax", ode_solver_method="bdf",
rtol=1e-8, atol=1e-12,
max_major_step_length=T_EVAL / 50, max_major_steps=200,
buffer_length=4000,
)
res = simulate(sys, ctx, (0.0, T_EVAL), options=opts, recorded_signals=rec)
return float(np.asarray(res.outputs["y"])[-1, 2])
# Central finite difference (the cross-check).
k1_0 = 0.04
h = 1e-4 * k1_0
yp = y3_of_k1(k1_0 + h)
ym = y3_of_k1(k1_0 - h)
grad_fd = (yp - ym) / (2 * h)
print(f"central-difference d y_3 / d k_1 at k_1={k1_0}: {grad_fd:.6e}")
print(f" (y_3(k_1 + h)={yp:.6e}, y_3(k_1 - h)={ym:.6e}, h={h:.3e})")
central-difference d y_3 / d k_1 at k_1=0.04: 3.437991e+00 (y_3(k_1 + h)=8.927128e-01, y_3(k_1 - h)=8.926853e-01, h=4.000e-06)
The finite-difference gradient is the analytic ground truth; for a small-enough $h$ it is good to four significant digits. The interpretation is direct: a $1\%$ increase in $k_1$ produces a $\partial y_3 / \partial k_1 \times 0.04 \times 0.01 \approx $ several-percent shift in $y_3(10^4)$. That is the parameter $k_1$'s influence on the slow $A \to B$ rate, which is what we control by setting $k_1$ in the first place. The sign is positive — more $k_1$ means more $A$ converted to $B$, which downstream gets converted to $C$, so $y_3$ goes up.
Note. A pure
jax.grad(y3_of_k1)(k1_0)over the simulator is the natural autodiff path here; the framework supports it viaenable_autodiff=True. We use the central-difference cross-check for the published numbers because the autodiff path through a top-levelsimulaterequires a JIT recompile per call site (it cannot be cached statically), and that pre-JIT compile cost is much larger than the central-difference solve cost for a single scalar parameter. For multi-parameter gradient runs — the parameter-identification setting — the autodiff path is the right one. Seebattery_part_2_parameter_estimation_synthetic_data.ipynbfor the canonical workflow.
Failure modes¶
Each of the choices above has a regime where it goes wrong.
- BDF's Newton iteration can diverge if the initial guess for $y_{n+1}$ is far from a fixed point. In practice this happens when the prior step was large and the system has just entered a different regime — for example, exiting the quasi-steady-state into a fast equilibration. Symptoms: the simulator hits
max_major_stepsexceeded, or NaN appears in the state. The fix is usually tighteningrtol/atolso the BDF order-controller backs off and the Newton iteration converges; less commonly, manually clampingmax_minor_step_size. - Dopri5 silently produces wrong answers if you set
max_major_stepstoo generously and don't notice that wall-time exploded. This is a worse failure than crashing because the trajectory looks fine — the per-step error is bounded byrtol, but the wall-time per second of model time is hundreds of times larger than it should be. The diagnostic is sample count: if your recording has more samples than you would expect from the visible dynamics, the solver is over-resolving a fast decayed mode. - float32 cannot represent the Robertson steady state. At $t = 10^{11}$, $y_2 \sim 8 \times 10^{-14}$ — below the smallest representable float32 normal ($\sim 1.2 \times 10^{-38}$ in absolute terms, but with relative-precision $\sim 10^{-7}$ the floor is much worse). We force
jax_enable_x64at the top of the notebook for this reason; if you turn it off, $y_2$ underflows to zero around $t = 10^4$ and the BDF Newton iteration loses convergence at that point. - The Robertson problem is a pure ODE — no algebraic constraint. Stiffness alone is BDF's wedge here. If your problem is a DAE (an algebraic constraint $g(\mathbf{y}) = 0$ alongside the dynamics), you also need projection or stabilization to keep the constraint residual bounded — see
dae_projection_pendulum.ipynb. The mass-conservation $y_1 + y_2 + y_3 = 1$ here is not a DAE constraint; it's a linear invariant of the ODE itself, which BDF preserves to machine precision essentially as a side effect of solving (7) exactly. - At extreme stiffness ratios (~$10^{12}$ or beyond), even BDF's accuracy can suffer. Specifically, the BDF order-controller may drop to order 1 (backward Euler), at which point the time integration is first-order accurate. The fix is tightening
rtolso the controller stays in its higher-order regime — see Hairer & Wanner (1996) Vol II §IV.5 for the order-selection logic.
Exercises¶
- (Easy) Reduce $k_2$ from $3 \times 10^7$ to $3 \times 10^5$ (a hundredfold less stiff) and re-run the Dopri5 horizon sweep. At what $T$ does Dopri5 stay viable now? Hint: the stability step bound should relax by a factor of 100, so horizons that were impossible at $k_2 = 3 \times 10^7$ should be feasible at $k_2 = 3 \times 10^5$.
- (Easy) Add a constant input flux of $A$: extend
Robertson._odeso $\dot y_1 = -k_1 y_1 + k_3 y_2 y_3 + s$ for a constant source $s$. With $s = 10^{-4}$, what does the steady state become? Hint: the system is now non-autonomous and never reaches a true $y_1 \to 0$ — it sits at a non-trivial fixed point set by $k_1 y_1 = s$. - (Medium) Compute the gradient of $y_3(10^4)$ with respect to all three rate constants $(k_1, k_2, k_3)$ — using either
jax.grad(autodiff) or central differences. Which rate constant most affects the slow dynamics? Hint: think about which rate gates the slow $A \to B \to C$ path. $k_1$ is the natural answer; verify with the numbers. - (Medium) Pick a different stiff system from the chemistry literature — e.g. the HIRES photolysis problem (8 species, mildly stiff) or the OREGO Belousov-Zhabotinsky oscillator (3 species, very stiff with limit-cycle behaviour). Build it as a
LeafSystemand rerun the full §3-4 analysis. Where does the spectral ratio sit? - (Open-ended) Jaxonomy's BDF solver is the only implicit solver shipped today. Implement a Radau IIA (or SDIRK4) integrator following the pattern in
jaxonomy/backend/_jax/bdf.py. On which problems does Radau win — and by how much? Hint: the answer involves the order-vs-cost trade-off. Radau IIA at order 5 needs $\sim 3 \times$ as many flops per step as BDF but takes proportionally larger steps because of better accuracy per step. On a problem where the step is accuracy-bound (rather than stability-bound), Radau wins; on a problem where the step is stability-bound (where both are already at the stability limit), they tie or BDF wins because of cheaper per-step cost.
Key takeaways¶
- Robertson 1966 is the canonical stiff-ODE test. A 3-state chemistry with rate constants spanning 9 orders of magnitude. BDF integrates the full $t \in [0, 10^{11}]$ horizon in well under a second; Dopri5 takes >0.1 s of wall-time per 1 s of model time and never finishes the headline horizon.
- Stiffness is about the step required for stability, not for accuracy. Explicit solvers (Dopri5) are conditionally stable; their step shrinks to track the fastest eigenvalue of the Jacobian, even after that mode has visibly decayed. Implicit solvers (BDF) are A-stable; they step at whatever size accuracy demands.
- Robertson's Jacobian sits at $|\lambda_{\max}| \sim 10^3$ during the quasi-steady-state regime, giving a spectral ratio in excess of $10^{19}$. The crossover on the $k_2$ axis — where Dopri5's step count first exceeds BDF's — sits in the $k_2 \in [10^5, 10^6]$ band, a stiffness ratio $k_2/k_1 \sim 10^7$. Below that, Dopri5 is the simpler choice; above that, BDF is the only choice that scales.
- BDF's per-step cost is one Jacobian evaluation + a Newton solve. Jaxonomy's BDF builds the Jacobian via
jax.jacrev, so the user pays no extra code. For RHS-cheap problems (Robertson), this is a $5{-}10\times$ per-step overhead — fully amortised by larger step sizes. For RHS-expensive problems, the BDF tax is even smaller in relative terms. - The BDF step is differentiable. Parameter identification on stiff systems is one
jax.gradaway once the simulator returns the right kind of trace. - The Robertson problem is a pure ODE, not a DAE. Mass conservation $\sum y_i = 1$ is a linear invariant that BDF preserves automatically (we measured $|\sum y_i - 1| < 10^{-15}$ — float64 round-off). For stiff problems with explicit algebraic constraints, see
dae_projection_pendulum.ipynb.
Next steps¶
dae_projection_pendulum.ipynb— when stiffness is paired with an algebraic constraint; the canonical 1-hour pendulum drift study under BDF + projection.multi_domain_hvac.ipynb— long-horizon BDF-DAE on an acausal thermal network; the engineering payoff of A-stable integration over 24-hour building-energy simulations.linearization_workflow.ipynb—findop+linearizeuse BDF under the hood to converge on a trim point; the same Newton-iteration story applied to operating-point identification rather than time integration.- The booster series (
booster_part_1_modeling.ipynb) — BDF underpins the acausal multibody propagator for the F1 rocket, where the equations of motion are stiff because of high-frequency tank-slosh modes.
References¶
- Robertson, H. H. (1966). The Solution of a Set of Reaction Rate Equations. In: J. Walsh (ed.), Numerical Analysis: An Introduction, Academic Press, pp. 178–182. — origin of the canonical 3-species test problem.
- Curtiss, C. F. and Hirschfelder, J. O. (1952). Integration of stiff equations. PNAS 38(3): 235–243. — first formal definition of "stiff" in the numerical-analysis sense.
- Hairer, E. and Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.). Springer-Verlag. — Vol. II §IV.1 is the textbook treatment of Robertson; §IV.5–7 cover BDF and Radau in detail.
- Shampine, L. F. (1994). Numerical Solution of Ordinary Differential Equations. Chapman & Hall. Ch. 6. — engineer-facing analysis of when to choose stiff vs non-stiff solvers, with worked examples.
- Lambert, J. D. (1991). Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. Wiley. Ch. 6. — formal A-stability / L-stability theory and the BDF-order-vs-stability trade-off.
- Mazzia, F. and Iavernaro, F. (2003). Test Set for Initial Value Problem Solvers. University of Bari. — modern stiff-ODE test corpus, including Robertson alongside HIRES, OREGO, VDPOL.