Stiff Systems

Stiff ODEs combine widely separated time scales; explicit solvers require impractically small steps, so implicit methods are essential.

Widely separated time scales

A system is stiff when it contains dynamics on very different time scales, some components decaying far faster than the scale of interest. Even after the fast transients have died, an explicit solver must keep steps small enough to remain stable with respect to those fast modes, wasting enormous effort.

A concrete symptom

On the test equation y' = lambda y with a large negative lambda, forward Euler is stable only for h below 2/|lambda|. If lambda is -10^6, steps must stay below 2 microseconds even to integrate over seconds, though the solution is nearly constant. The stability limit, not accuracy, dictates the step.

Implicit methods to the rescue

A-stable implicit methods, such as backward Euler, the implicit trapezoidal rule, and the BDF family, have stability regions that cover the entire left half-plane. They remain stable for any step size on decaying problems, so the step can be chosen for accuracy alone. The trade-off is solving a (possibly nonlinear) system each step via Newton iteration, which needs the Jacobian.

from scipy.integrate import solve_ivp
import numpy as np
# Robertson problem: classic stiff test
def f(t, y):
    return [-0.04*y[0]+1e4*y[1]*y[2],
             0.04*y[0]-1e4*y[1]*y[2]-3e7*y[1]**2,
             3e7*y[1]**2]
sol = solve_ivp(f, [0, 1e4], [1,0,0], method='BDF')

Recognizing stiffness

Stiffness is not a property of the equation alone but of the equation, the interval, and the accuracy required. A practical sign is an explicit adaptive solver taking many rejected or tiny steps. The remedy is to switch to an implicit solver (BDF or Radau) that supplies or approximates the Jacobian.

Stiffness is pervasive in reacting and magnetized-plasma physics, where fast atomic or gyration time scales coexist with slow transport, making implicit integrators standard in breeder Hyperion modeling.