🛰️Ordinary Differential Equations: Initial Value Problems

Stiff Equations and Implicit Methods

Stiffness and stability regions, backward Euler and the trapezoidal rule, A- and L-stability, BDF, Newton inside each step

Madison, 1952 — Berkeley, 1973

Charles Curtiss and Joseph Hirschfelder, chemists at the University of Wisconsin, publish Integration of Stiff Equations in the Proceedings of the National Academy of Sciences. Their problem is chemical kinetics: in a flame or a reactor some radicals form and disappear in microseconds, while the overall composition evolves over seconds. Once the fast species settle into balance, nothing in the solution changes quickly. Yet every explicit method they try keeps taking microsecond steps, for ever, as if the fast reactions were still happening. They name the phenomenon — stiffness — and propose an implicit method that escapes it. In the early 1970s Bill Gear at Illinois turns implicit backward-differentiation formulas into robust, variable-order software (1971). In 1973 Laurence Nagel and Donald Pederson at Berkeley release SPICE, a circuit simulator whose equations are stiff for the same reason: picosecond transistor time constants inside microsecond circuits. SPICE uses implicit integration — the trapezoidal rule and Gear's formulas — and every integrated circuit designed since has been simulated by its descendants.

Recall lesson 9-3's test equation

y' = \lambda y

with

\operatorname{Re}\lambda < 0

: one step multiplies

y

by an amplification factor

R(z)

z = h\lambda

. Explicit Euler gives

R = 1 + z

: the numerical solution decays only if

|1 + h\lambda| \le 1

. For

\lambda = -1000

that means

h < 0.002

— whatever accuracy we want. The set of

z

with

|R(z)| \le 1

is the method's region of absolute stability; every explicit method has a bounded one.

💡What stiffness is

A system is stiff when it contains decaying modes much faster than the behaviour we want to follow. For a linear system

\mathbf y' = A\mathbf y

the modes are the eigenvalues of

A

(Unit 6): if

\lambda_{fast} = -10^6

and the solution of interest varies on a time scale of 1, an explicit method needs

h \approx 10^{-6}

for stability although

h \approx 10^{-2}

would be accurate. Stiffness is not about the solution being hard — it is about the method's stability limit being far stricter than its accuracy needs. On

y' = -1000(y - \cos t) - \sin t

, whose solution is just

\cos t

after a 1 ms transient, explicit Euler with

h = 0.0019

is fine and with

h = 0.0021

the error is

10^{33}

Derivation

Why one test equation speaks for every linear system

Take

\mathbf y' = A\mathbf y

with

A = V\Lambda V^{-1}

(eigenvectors in

V

, eigenvalues on the diagonal of

\Lambda

, Unit 6). Change variables to

\mathbf w = V^{-1}\mathbf y

, the components along the eigenvectors:

\mathbf w' = V^{-1}A\mathbf y = V^{-1}AV\,\mathbf w = \Lambda\mathbf w, \qquad\text{that is}\qquad w_i' = \lambda_iw_i \ \text{ for each } i.

The system falls apart into independent copies of the test equation, one per eigenvalue. Euler, RK4, the Adams methods and the implicit methods below are all linear in

\mathbf y

when

f

is linear, so applying them to

\mathbf y

is the same as applying them to each

w_i

separately. The step is stable only if every

z_i = h\lambda_i

lies in the stability region — including the fastest mode, even after it has decayed to nothing. With

\lambda = -10^6

in the system, an explicit method needs

h \lesssim 2 \times 10^{-6}

for ever. For a nonlinear system the same reasoning applies locally, with the eigenvalues of the Jacobian

\partial\mathbf f/\partial\mathbf y

\text{backward Euler: } y_{n+1} = y_n + hf(t_{n+1}, y_{n+1}), \quad R(z) = \frac{1}{1-z};\qquad \text{trapezoidal: } y_{n+1} = y_n + \tfrac h2\big(f_n + f_{n+1}\big),\quad R(z) = \frac{1 + z/2}{1 - z/2}

Implicit methods evaluate the slope at the unknown new point. Each step requires solving an equation — linear for linear problems, otherwise by Newton's method (Unit 2) with the Jacobian ∂f/∂y.

Derivation

The amplification factors, and Newton inside the step

Backward Euler on

y' = \lambda y

y_{n+1} = y_n + h\lambda y_{n+1}

. Solve for the new value:

y_{n+1} = \frac{1}{1 - z}\,y_n.

For

\operatorname{Re}z \le 0

|1 - z| \ge 1

(the point

z

is at least as far from 1 as 0 is), so

|R| \le 1

: stable for every step. As

z \to -\infty

R \to 0

The trapezoidal rule:

y_{n+1} = y_n + \tfrac z2(y_n + y_{n+1})

, so

(1 - \tfrac z2)y_{n+1} = (1 + \tfrac z2)y_n

and

R(z) = \frac{1 + z/2}{1 - z/2}.

|R| \le 1

exactly when

z

is at least as close to

-2

as to

+2

— the whole left half-plane. As

z \to -\infty

R \to -1

Newton inside the step. For a nonlinear

f

, backward Euler's new value is the root of

G(\mathbf y) = \mathbf y - \mathbf y_n - h\,\mathbf f(t_{n+1}, \mathbf y) = \mathbf 0

. Newton's method (lesson 2-4) needs its Jacobian,

I - h\,\partial\mathbf f/\partial\mathbf y

, and solves

(I - hJ)\boldsymbol\Delta = -G

at each iteration, starting from the previous step's value — which is usually so close that one to three iterations suffice. The matrix

I - hJ

changes slowly, so solvers factorise it once (Unit 3) and reuse the factors for many steps.

💡A-stability and L-stability

A method is A-stable if its stability region contains the entire left half-plane: no decaying mode can make it blow up, whatever

h

. Backward Euler and the trapezoidal rule are both A-stable. Dahlquist's second barrier (1963): no explicit method, and no linear multistep method of order above 2, is A-stable. For very stiff modes a stronger property matters. As

z \to -\infty

, backward Euler's

R(z) \to 0

— the stiff mode is damped out in one step (L-stability). The trapezoidal rule's

R(z) \to -1

: a stiff mode is not amplified, but it rings, flipping sign every step and decaying only slowly. With

z = -100

, backward Euler reduces a stiff component by a factor 100 per step; the trapezoidal rule only to 0.96 of its size, with alternating sign. Circuit simulators know this ringing well.

BDF methods. The backward differentiation formulas differentiate an interpolating polynomial through

y_{n+1}, y_n, \dots

and set its derivative at

t_{n+1}

equal to

f(t_{n+1}, y_{n+1})

. BDF1 is backward Euler; BDF2 is

\tfrac32y_{n+1} - 2y_n + \tfrac12y_{n-1} = hf_{n+1}

. Orders 1 and 2 are A-stable; orders 3–6 are stable in almost the whole left half-plane — enough for most stiff problems. They are the core of Gear's codes, of solve_ivp(method='BDF') and of LSODA. Radau IIA implicit Runge–Kutta methods (method='Radau') are A- and L-stable at order 5, at the price of solving a larger system each step.

Derivation

BDF2 from a backward difference

Approximate

y'(t_{n+1})

using the new value and the two before it. Expand the old values about

t_{n+1}

y_n = y_{n+1} - hy' + \tfrac{h^2}{2}y'' - \dots, \qquad y_{n-1} = y_{n+1} - 2hy' + 2h^2y'' - \dots

Look for the combination that keeps

y'

and cancels

y''

3y_{n+1} - 4y_n + y_{n-1}

. The

y

terms give

3 - 4 + 1 = 0

, the

y'

terms

4h - 2h = 2h

, the

y''

terms

-2h^2 + 2h^2 = 0

. So

y'(t_{n+1}) = \frac{3y_{n+1} - 4y_n + y_{n-1}}{2h} + O(h^2).

Setting this equal to

f(t_{n+1}, y_{n+1})

and multiplying by

h

gives

\tfrac32y_{n+1} - 2y_n + \tfrac12y_{n-1} = hf_{n+1}

— BDF2. (This one-sided formula is also the derivative of the quadratic through the three points, which is how the higher BDF formulas are built.) Like backward Euler it evaluates

f

only at the new point, which is what gives it stability far into the left half-plane.

My methods are really methods of working and thinking; this is why they have crept in everywhere anonymously.

– Emmy Noether, letter to Helmut Hasse (1931)

📐Robertson's kinetics

H. H. Robertson's 1966 test problem: three species, rate constants

0.04

10^4

and

3 \times 10^7

— nine orders of magnitude apart. solve_ivp with RK45 needs 7,088 steps and 49,718 evaluations just to reach $t = 10$ , and about 1.7 million steps to reach

t = 1000

. Radau reaches

t = 1000

in 115 steps, BDF in 211 steps with 569 evaluations, LSODA (which switches automatically between Adams and BDF) in 327 steps — all agreeing that

y_1(1000) = 0.33687

, and conserving total mass to

10^{-16}

. Supplying the Jacobian jac= lets the implicit methods' Newton iterations converge quickly.

💡Where you will use this

Rocket Propulsion: finite-rate chemistry in combustion chambers and nozzles is the classic stiff problem. Power Electronics and Circuits: switching converters mix nanosecond switching with millisecond control loops. Structures: very stiff elements beside flexible ones make the semi-discretised equations of motion stiff; implicit Newmark methods are the norm. Heat transfer and PDEs (Unit 10): the method of lines turns diffusion into a stiff ODE system whose stiffness grows as the mesh is refined.

📚Further reading

Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed., Springer, 1996) — stiffness, A- and L-stability, BDF and Radau methods.
C. F. Curtiss and J. O. Hirschfelder, 'Integration of stiff equations', Proceedings of the National Academy of Sciences 38(3), 235–243 (1952) — the paper of the cold open.
C. William Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, 1971) — the book that turned BDF methods into software.

↓But there is a problem

Implicit methods make stiff problems tractable by being stable: errors decay instead of grow. But stability is a weak promise. A method can be perfectly stable and still slowly, systematically get the physics wrong — losing energy from an orbit that should conserve it, or gaining it. In 1907 the Norwegian mathematician Carl Størmer wanted to know why the aurora appears in rings around the magnetic poles. He computed, by hand, the paths of charged particles spiralling in the Earth's magnetic field — hundreds of trajectories over thousands of hours of calculation. The method he chose was so simple that physicists rediscovered it sixty years later, and found that it preserved something every other method lost.