Stiff differential equations and implicit methods

The formula

The system being integrated, and the test equation everything gets judged against:

dy/dt = f(t, y)          general initial-value problem
dy/dt = λ·y              Dahlquist test equation, Re(λ) < 0

Reading: any solver's stiffness behavior is diagnosed on the single decaying mode e^{λt}. A real system's λ values are the eigenvalues of the Jacobian ∂f/∂y; the fastest one runs the show.

Explicit (forward) Euler versus backward Euler, the simplest implicit method:

y_{n+1} = y_n + h·f(t_n, y_n)              explicit: uses the slope you already have
y_{n+1} = y_n + h·f(t_{n+1}, y_{n+1})      implicit: the unknown appears on both sides

Reading: the implicit step is an equation to be solved, not a formula to be evaluated — each step costs a Newton iteration and a Jacobian, and buys unconditional stability.

Apply each to the test equation and the difference is one fraction:

explicit Euler:   y_{n+1} = (1 + h·λ)·y_n        stable only if |1 + h·λ| ≤ 1  →  h ≤ 2/|λ|
backward Euler:   y_{n+1} = y_n / (1 − h·λ)      |amplification| < 1 for every h > 0 when Re(λ) < 0

Reading: the explicit factor 1 + h·λ is a polynomial in h·λ — it always escapes the unit circle for a big enough step. The implicit factor 1/(1 − h·λ) only shrinks as h·|λ| grows. That property — stability region covering the whole left half-plane — is Dahlquist's A-stability ^[1].

The backward differentiation formulas (BDFs), the workhorse family behind Gear's method and ode15s. BDF1 is backward Euler; BDF2 is:

y_{n+2} − (4/3)·y_{n+1} + (1/3)·y_n = (2/3)·h·f(t_{n+2}, y_{n+2})

Reading: fit a polynomial through the last few solution points plus the unknown new one, differentiate it, and demand the derivative match f at the new point. Orders 1 through 6 exist; order 7 and up are unstable and unusable ^[2].

Where you meet it

The circuit model that brings ode45 to its knees. A power-supply or snubber model carries parasitic time constants in the nanoseconds while the waveform you care about plays out over milliseconds. An explicit solver must hold h ≲ 2·τ_fastest forever — a 1 ns parasitic in a 1 ms simulation is a floor of 500,000 steps, even though the parasitic transient died in the first 5 ns. This is exactly why SPICE has used implicit integration — trapezoidal and Gear — since the 1970s; explicit methods were never on the table for circuits ^[3].

Chemical kinetics and propellant burn models. Reaction networks stack rate constants spanning six or more orders of magnitude: radicals that equilibrate in microseconds alongside bulk species that evolve over seconds. This is the application that produced the field — the 1952 paper that named "stiff equations" came out of Navy Bureau of Ordnance work on kinetics ^[4]^[5].

Thermal networks on the test stand. A lumped model with a thin thermocouple sheath (time constant of milliseconds) bolted to a massive fixture (time constant of minutes) is stiff by construction. Same for structural models mixing very stiff and very soft elements.

The review board asking why the sim takes four hours. The signature of stiffness is unmistakable: the solution looks smooth and boring, and the explicit solver is taking millions of microscopic steps anyway. If ode45 reports step sizes thousands of times smaller than anything visibly changing, the problem is stiff. That is not a tuning problem — tightening or loosening tolerances won't fix it. Switch to ode15s and the same run finishes in seconds ^[6].

How it works

Stiffness is a mismatch of time scales plus a stability trap. After a fast mode decays, it contributes nothing to the solution — but its eigenvalue still sits in the Jacobian, and an explicit method's stability region is a bounded blob around the origin. For λ = −1000, explicit Euler is stable only up to h = 0.002; at h = 0.0021 the amplification factor is −1.1 per step, and by t = 5 the "solution" has grown to around 10⁹⁸ while the true answer decays quietly to zero. Higher-order explicit methods don't escape: classical RK4 buys you h·|λ| ≤ 2.785 on the negative real axis — barely a 40% improvement over Euler for four times the work per step ^[1].

Backward Euler at the same h·λ = −1000 multiplies the fast mode by about 1/1001 per step — it annihilates the transient in one stride and then takes steps sized to the slow physics. On a two-mode system with eigenvalues −1 and −1000 integrated to t = 5, explicit Euler needs at least 2,500 steps to survive; backward Euler at h = 0.01 takes 500 steps and lands the slow mode at 0.00691 against the exact 0.00674. Both numbers above were run, not recited.

What you pay: each implicit step is a nonlinear solve. The solver forms the Jacobian ∂f/∂y, factors a matrix, and Newton-iterates. On a non-stiff problem that overhead is pure loss — ode45 beats ode15s on clean, well-scaled dynamics, which is why the stiff solver isn't simply the default.

Two theoretical fences shape every production code. First, the second Dahlquist barrier: no A-stable linear multistep method can exceed order 2 ^[1]. Production BDF codes get to order 5 or 6 by accepting slightly less than full A-stability — stability regions that cover most, but not all, of the left half-plane — which works because real eigenvalues rarely hug the imaginary axis. Second, stability is not accuracy. Backward Euler is unconditionally stable and only first-order accurate; it will happily hand you a smooth, plausible, wrong answer at too coarse a step. Stable garbage is still garbage — converge the step size or the tolerance before believing the plot.

One more trap, familiar to anyone who has simulated a switching converter: the trapezoidal rule is A-stable and second-order, but its amplification factor tends to −1 (not 0) as h·λ → −∞ — at h·λ = −1000 it multiplies by −0.996 per step. Very stiff modes don't die; they flip sign every step and linger as a sawtooth riding the waveform. That is the "trapezoidal ringing" artifact in SPICE transient runs. Backward Euler and the BDFs drive the factor to zero instead — the property is called L-stability, and it's why stiff codes are built on BDFs rather than the trapezoidal rule alone.

In MATLAB the practical decision tree is short: try ode45; if it crawls on a smooth solution, use ode15s, which runs Klopfenstein–Shampine numerical differentiation formulas (NDFs) — a lightly modified BDF family — at variable order 1 through 5, with classic BDFs available as an option ^[6]^[9]. In Python, scipy.integrate.solve_ivp(method='BDF') or the LSODA option, which detects stiffness and switches automatically. The stiffness ratio max|Re λ| / min|Re λ| you'll see quoted is a rough screen, not a definition — it's neither necessary nor sufficient ^[1]. The operational test remains the honest one: if the step size is being set by stability rather than by anything happening in the solution, the problem is stiff.

History

The problem announced itself as a defense problem. In March 1952, Charles Curtiss and Joseph Hirschfelder — two theoretical chemists at the University of Wisconsin's Naval Research Laboratory (Curtiss had earned his chemistry Ph.D. under Hirschfelder just four years earlier), working under a Navy Bureau of Ordnance contract — published "Integration of Stiff Equations" in the Proceedings of the National Academy of Sciences, naming the phenomenon and pointing at where it lived: chemical kinetics, electrical circuit theory, and missile guidance ^[4]^[5]. Their proposed cure, computing the derivative at the forward point instead of the current one, is the backward differentiation formula in embryo ^[2]^[5].

The theory arrived a decade later from Sweden. Germund Dahlquist, who had spent the 1950s at the Swedish Board of Computer Machinery before joining Stockholm's Royal Institute of Technology, published "A special stability problem for linear multistep methods" in BIT in 1963 — the paper that defined A-stability, proved that no A-stable multistep method can beat second order, and became one of the most cited papers in numerical analysis ^[1]^[7]. His barriers told the field exactly what it could and couldn't have, and every stiff solver since has been engineered right up against them.

The software came from Illinois. C. William Gear, a London-born computer scientist at the University of Illinois, turned the Curtiss–Hirschfelder idea into robust variable-order BDF codes in papers from the late 1960s and his 1971 book Numerical Initial Value Problems in Ordinary Differential Equations — so thoroughly that BDF integration is still called "Gear's method," and it was built into the Berkeley SPICE2 circuit simulator that every modern SPICE descends from ^[2]^[3]^[8]. The lineage lands on your desktop in 1997, when Lawrence Shampine and Mark Reichelt published the MATLAB ODE suite, including ode15s and its NDF variant of the BDFs ^[6]^[9] — which is why the fix for a crawling simulation is, forty-plus years after Curtiss and Hirschfelder, a five-character edit.

Related tools

/tools/rc-charge-time
/tools/rc-snubber
/tools/rc-filter
/tools/half-life-decay
/tools/heat-conduction

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →