Time stepping: accuracy, stability and stiffness

Before this lesson: eigenvalues and diagonalization (linear algebra 05); the matrix exponential and linear ODEs (linear algebra 16); Taylor expansion (calculus 03) and finite differences (calculus 07). Complex numbers $z = x + iy$ and $\lvert z\rvert$. Time: about 90 minutes, plus the notebook.

---

Why this matters #

• Every weather or climate emulator is trained on the output of a time-stepping solver, and many hybrid models put a learned term inside such a solver.
• An autoregressive ML weather model (a network that maps the state now to the state 6 hours later) is itself an explicit one-step method. Long rollouts stay bounded only if errors are not amplified step after step.
• A neural ODE is trained through an adaptive ODE solver. When the learned dynamics become stiff, the solver takes more and more steps, and training slows down.
• The stability limit often dominates the cost of a simulation. For explicit diffusion, $\Delta t \propto \Delta x^2$: halve the grid spacing and you need 4 times more steps.
• Diffusion-model samplers such as DPM-Solver treat the linear part of their ODE exactly, an idea from the stiff-solver toolbox of section 3.

1. Predict first #

2. The picture #

A time stepper is a map #

Take an ODE $u' = f(t, u)$ with a state $u(t) \in \R^n$ and a start value $u(0) = u_0$. A one-step method replaces the true flow over a time $h$ by a map $u_{k+1} = \Phi_h(u_k)$.

• Forward Euler follows the tangent: $u_{k+1} = u_k + h\,f(t_k, u_k)$.
• Backward Euler uses the slope at the end of the step: $u_{k+1} = u_k + h\,f(t_{k+1}, u_{k+1})$. The unknown $u_{k+1}$ appears on both sides, so each step solves an equation. Such a method is implicit; forward Euler is explicit.
• RK4 (the classical Runge–Kutta method) evaluates $f$ four times per step and combines the slopes:

Order. If the error at a fixed final time behaves like $C h^p$, the method has order $p$. On log–log axes, $\log e = p\log h + \log C$ is a straight line with slope $p$. So the order is a slope you can read from a plot.

Interactive figure: loading…

The test equation and the stability function #

Accuracy is one half of the story. The other half is stability. Start with the simplest ODE,

with a number $\lambda$ that may be complex. The exact solution is $y(t) = e^{\lambda t}$, so one step of size $h$ multiplies it by $e^{h\lambda}$. A one-step method applied to this ODE also multiplies $y$ by a fixed number per step:

$R$ is the stability function of the method. For example, forward Euler gives $y_{k+1} = y_k + h\lambda y_k = (1 + z)\,y_k$, so $R(z) = 1 + z$. Backward Euler gives $y_{k+1} = y_k + z\,y_{k+1}$, so $y_{k+1} = y_k/(1 - z)$ and $R(z) = 1/(1 - z)$. Two questions decide everything:

• Accuracy: is $R(z)$ close to $e^z$ for small $z$? A method of order $p$ has $R(z) = e^z + O(z^{p+1})$.
• Stability: is $\lvert R(z)\rvert \le 1$? If not, $\lvert y_k\rvert = \lvert R(z)\rvert^k$ grows without bound, even when the true solution decays ($\operatorname{Re}\lambda < 0$).

Method	$R(z)$	stable on the negative real axis	$R(-\infty)$
forward Euler (FE)	$1 + z$	$-2 \le z \le 0$	$-\infty$
backward Euler (BE)	$1/(1 - z)$	all $z \le 0$	$0$
trapezoid (Crank–Nicolson, CN)	$(1 + z/2)/(1 - z/2)$	all $z \le 0$	$-1$
RK4	$1 + z + z^2/2 + z^3/6 + z^4/24$	$-2.785 \le z \le 0$	$+\infty$

The set of $z$ with $\lvert R(z)\rvert \le 1$ is the stability region. It is a picture in the complex plane, and every method has its own.

Interactive figure: loading…

Why one scalar equation is enough. Take a system $u' = f(u)$ and a small perturbation $\delta$ of a solution. To first order it obeys $\delta' = J\delta$, where $J = \partial f/\partial u$ is the Jacobian. If $J = V\Lambda V^{-1}$ is diagonalizable, the coordinates $V^{-1}\delta$ decouple into scalar test equations, one per eigenvalue $\lambda_i$ (linear algebra 05). So the rule is: every $h\lambda_i$ must lie in the stability region. This is exact for a linear system with a normal matrix. For a nonlinear or strongly non-normal problem it is a rule of thumb.

PDEs become big ODE systems. Discretize space first and keep time continuous: this is the method of lines. The heat equation $u_t = \nu u_{xx}$ ($\nu$ is the diffusivity) on a grid with spacing $\Delta x$ becomes $u' = \nu A u$, where $A = \frac{1}{\Delta x^2}\operatorname{tridiag}(1, -2, 1)$ is the second-difference matrix. Its eigenvalues lie between about $-4/\Delta x^2$ and $0$. A fine grid creates very fast decaying modes, and their $h\lambda$ must stay inside the stability region. This is where $\Delta t \propto \Delta x^2$ comes from (section 3).

3. Definitions and main results #

Consistency, order and convergence #

Consider $y' = f(t, y)$ and a one-step method $y_{k+1} = y_k + h\,\phi(t_k, y_k, h)$; forward Euler has $\phi = f$. Put the exact solution $y(t)$ into the method. What is left over, divided by $h$, is the local truncation error

• The method is consistent if $\tau \to 0$ as $h \to 0$. For a continuous $\phi$ this means $\phi(t, y, 0) = f(t, y)$.
• It has order $p$ if $\tau = O(h^p)$. Forward Euler has $p = 1$, RK4 has $p = 4$.
• Many texts count the error of one step instead, $h\tau = O(h^{p+1})$, and call it the "local error of order $p + 1$". It is the same fact with different bookkeeping.

Example (forward Euler has order 1). Taylor: $y(t + h) = y(t) + h y'(t) + \tfrac{h^2}{2}y''(\xi)$. With $y' = f$, $\tau = \frac{y(t + h) - y(t)}{h} - f(t, y(t)) = \tfrac h2 y''(\xi) = O(h)$.

Theorem (convergence of one-step methods). Assume an exact start $y_0 = y(t_0)$, $\lvert\tau_k\rvert \le C h^p$ for every $k$, and that $\phi$ is Lipschitz in $y$ with constant $K > 0$: $\lvert\phi(t, y, h) - \phi(t, \tilde y, h)\rvert \le K\lvert y - \tilde y\rvert$. Then the global error $e_k = y(t_k) - y_k$ satisfies

Why this is true, in words. Each step adds a local error of size $h^{p+1}$. There are $T/h$ steps, so together they give $h^p$. The factor $e^{KT}$ bounds how much the earlier errors are amplified on the way. The amplification step is the stability part of the proof.

Proof of the convergence theorem

1. By the definition of $\tau$: $y(t_{k+1}) = y(t_k) + h\phi(t_k, y(t_k), h) + h\tau_{k+1}$.
2. Subtract the method $y_{k+1} = y_k + h\phi(t_k, y_k, h)$:

$$e_{k+1} = e_k + h\,[\phi(t_k, y(t_k), h) - \phi(t_k, y_k, h)] + h\tau_{k+1}.$$

3. Lipschitz: $\lvert e_{k+1}\rvert \le (1 + hK)\,\lvert e_k\rvert + C h^{p+1}$.
4. Unroll from $e_0 = 0$ and sum the geometric series:

$$\lvert e_k\rvert \le C h^{p+1} \sum_{j=0}^{k-1}(1 + hK)^j = \frac{C h^p}{K}\big((1 + hK)^k - 1\big).$$

5. Use $1 + hK \le e^{hK}$ and $kh = t_k - t_0$. $\blacksquare$

Consistency + stability = convergence. This is the central idea of numerical analysis. For linear multistep methods it is the Dahlquist equivalence theorem: a method converges if and only if it is consistent and zero-stable. (Zero-stable means that small perturbations stay bounded on a fixed interval as $h \to 0$; one-step methods with a Lipschitz $\phi$ are always zero-stable.) For linear PDE schemes it is the Lax equivalence theorem: a consistent scheme for a well-posed linear initial-value problem converges if and only if it is stable.

Measuring the order. If $e(h) \approx C h^p$, then with two step sizes $p \approx \log_2\big(e(h)/e(h/2)\big)$. With many, fit a line to $\log e$ against $\log h$. The fit is meaningful only in the asymptotic range: $h$ small enough that the leading term dominates, but large enough that the error is far above rounding (about $10^{-15}$ relative).

Absolute stability, A-stability, L-stability #

Convergence is a statement about $h \to 0$. In practice $h$ is fixed, and we need absolute stability: the method is absolutely stable at $z = h\lambda$ if $\lvert R(z)\rvert \le 1$.

• Forward Euler: the disk $\lvert 1 + z\rvert \le 1$ (centre $-1$, radius 1).
• Backward Euler: everything outside the open disk $\lvert 1 - z\rvert < 1$; it contains the whole left half-plane.
• Trapezoid: exactly the closed left half-plane, because $\lvert 1 + z/2\rvert \le \lvert 1 - z/2\rvert$ means that $z$ is at least as close to $-2$ as to $+2$.
• RK4: a bounded region; on the real axis $[-2.785, 0]$, on the imaginary axis $\lvert\operatorname{Im} z\rvert \le 2\sqrt{2} \approx 2.83$.

The stability function of an explicit method is a polynomial, so $\lvert R(z)\rvert \to \infty$ as $\lvert z\rvert \to \infty$. Every explicit method has a bounded stability region, so a very negative eigenvalue always forces a small step.

Definition. A method is A-stable if its stability region contains the whole left half-plane $\operatorname{Re} z \le 0$: every decaying mode decays numerically, for every $h$. A method is L-stable if it is A-stable and $R(z) \to 0$ as $z \to -\infty$: very fast modes are not only stable, they are damped almost to zero in one step, as in reality.

Backward Euler is A-stable and L-stable. The trapezoid rule is A-stable but not L-stable: $R(-\infty) = -1$, so stiff modes keep their size and flip sign every step.

Dahlquist's second barrier and stiff solvers of higher order

An A-stable linear multistep method is implicit and has order at most 2 (Dahlquist's second barrier). Among them, the trapezoid rule has the smallest error constant. Higher-order stiff solvers therefore use implicit Runge–Kutta methods (Radau IIA, which is L-stable and of order $2s - 1$ with $s$ stages) or BDF methods. BDF1 (backward Euler) and BDF2 are A-stable; BDF3 to BDF6 are only $A(\alpha)$-stable: stable in a wedge around the negative real axis, not in the whole left half-plane.

The explicit heat equation: $\Delta t \le \Delta x^2/(2\nu)$ #

Theorem. The scheme $u_j^{m+1} = u_j^m + r\,(u_{j+1}^m - 2u_j^m + u_{j-1}^m)$ with $r = \nu\Delta t/\Delta x^2$ (grid index $j$, time level $m$) is stable for every Fourier mode if and only if $r \le 1/2$.

Proof (von Neumann analysis).

1. The scheme is linear with constant coefficients, so Fourier modes do not mix. Insert one mode $u_j^m = g^m e^{i\kappa x_j}$ (amplification factor $g$, wavenumber $\kappa$, $i = \sqrt{-1}$, $x_j = j\Delta x$).
2. Divide by $g^m e^{i\kappa x_j}$, and use $e^{i\theta} + e^{-i\theta} = 2\cos\theta$ and $1 - \cos\theta = 2\sin^2(\theta/2)$:

3. $g$ is real and $g \le 1$ always. Stability needs $g \ge -1$ for every $\kappa$.
4. The worst mode is the grid-scale zigzag (sawtooth), with $\sin^2 = 1$: $1 - 4r \ge -1 \iff r \le 1/2$, that is $\Delta t \le \Delta x^2/(2\nu)$. $\blacksquare$

The same number in the method-of-lines picture: $g = 1 + \Delta t\,\nu\lambda_k$, where $\nu\lambda_k = -(4\nu/\Delta x^2)\sin^2(\cdot)$ are the eigenvalues of $\nu A$. So von Neumann analysis is the stability region of forward Euler, applied to the spectrum of the space operator.

Interactive figure: loading…

Three remarks.

• On a finite grid with fixed ends, the largest $\lvert\lambda_k\rvert$ is slightly below $4/\Delta x^2$, so the exact limit is slightly above $\Delta x^2/(2\nu)$. The bound is sharp as the grid is refined.
• In $d$ dimensions with equal spacing the limit is $\Delta x^2/(2d\nu)$.
• With $r$ fixed, $\Delta t \propto \Delta x^2$, so the $O(\Delta t)$ error of forward Euler is $O(\Delta x^2)$, the same size as the spatial error. The price is the number of steps.

The CFL condition #

Take the advection equation $u_t + c\,u_x = 0$ with speed $c > 0$. Its solution $u(x, t) = u_0(x - ct)$ moves to the right with speed $c$. The value at $(x, t)$ depends only on the initial value at the point $x - ct$, the foot of the characteristic line through $(x, t)$. This point is the domain of dependence.

The upwind scheme $u_j^{m+1} = u_j^m - \mathrm{Co}\,(u_j^m - u_{j-1}^m)$, with the Courant number $\mathrm{Co} = c\,\Delta t/\Delta x$, uses only the grid point itself and its left neighbour (the side the wind comes from). After $m$ steps, $u_j^m$ depends on the initial data in $[x_j - m\Delta x,\ x_j]$: this interval is the numerical domain of dependence.

CFL condition (Courant–Friedrichs–Lewy). A scheme can converge only if its numerical domain of dependence contains the exact one as $\Delta x, \Delta t \to 0$. Here the point $x_j - c\,m\Delta t$ must lie in $[x_j - m\Delta x,\ x_j]$, so $\mathrm{Co} \le 1$. In words: in one step, information may not travel further than the stencil reaches.

Interactive figure: loading…

Von Neumann analysis agrees: $g = 1 - \mathrm{Co} + \mathrm{Co}\,e^{-i\kappa\Delta x}$ gives $\lvert g\rvert^2 = 1 - 2\,\mathrm{Co}(1 - \mathrm{Co})(1 - \cos\kappa\Delta x)$, which is $\le 1$ for every $\kappa$ exactly when $0 \le \mathrm{Co} \le 1$.

Stiffness #

Definition. A problem is stiff when the step size is set by stability, not by accuracy. The typical cause: the Jacobian has eigenvalues with large negative real parts (fast modes that decay very quickly) next to the slow dynamics you care about. The stiffness ratio $\max\lvert\operatorname{Re}\lambda\rvert / \min\lvert\operatorname{Re}\lambda\rvert$ measures the gap. Stiffness depends on the question: if you need to resolve the fast transient itself, small steps are needed anyway, and an explicit method is fine.

Interactive figure: loading…

Typical stiff problems: diffusion on fine grids (ratio $\propto 1/\Delta x^2$), chemical kinetics with rates over many orders of magnitude, and vertical mixing in an atmospheric column. Fast gravity and sound waves in a weather model are the oscillatory kind: large imaginary eigenvalues that carry little energy but set the explicit step.

A practical symptom: an adaptive explicit solver takes thousands of tiny steps on a solution that looks smooth. It has hit the stability boundary, not an accuracy problem.

Choosing a time stepper #

Many PDEs, after discretization in space, have the semilinear form

with a stiff linear part $\mathcal{L}$ (diffusion, hyperdiffusion, dispersion, fast linear waves) and a non-stiff nonlinear part $\mathcal{N}$ (advection, reactions). There are four families of methods.

Family	Examples	Cost per step	Use when
explicit	RK4, adaptive Dormand–Prince ("RK45")	a few evaluations of $f$	accuracy limits the step anyway: non-stiff ODEs, advection at $\mathrm{Co}$ near 1
implicit	backward Euler, BDF, Radau IIA	a linear or Newton solve	stiff nonlinear problems (chemistry); choose L-stable when stiff modes must be damped
IMEX	IMEX Euler, CNAB2, IMEX Runge–Kutta	a linear solve with $\mathcal L$	stiff linear part, non-stiff nonlinear part
exponential (ETD)	exponential Euler, ETDRK4	applying $e^{h\mathcal L}$ and $\varphi$-functions	$e^{h\mathcal L}$ is cheap, e.g. $\mathcal L$ diagonal in Fourier space

• IMEX (implicit–explicit) treats $\mathcal{L}$ implicitly and $\mathcal{N}$ explicitly. The simplest is IMEX Euler, $(I - h\mathcal{L})\,u_{k+1} = u_k + h\,\mathcal{N}(u_k)$: one linear solve per step (tridiagonal for 1-D diffusion, diagonal in Fourier space). The diffusion limit disappears; the explicit part keeps its own limit (for example the advective CFL condition). Semi-implicit weather models treat fast gravity waves this way.
• Exponential integrators solve the linear part exactly with the variation-of-constants formula

and approximate only the integral. With $\mathcal{N}$ frozen at $u_k$ you get exponential Euler (ETD1): $u_{k+1} = e^{h\mathcal{L}}u_k + h\,\varphi_1(h\mathcal{L})\,\mathcal{N}(u_k)$, with $\varphi_1(z) = (e^z - 1)/z$. The stiff linear part is integrated exactly, so it imposes no step limit.

A short rule: if accuracy sets the step, go explicit; if stability sets the step and the fast modes do not matter, go implicit, IMEX or exponential. An implicit step costs a linear or Newton solve. It is worth it when the explicit step would be smaller by more than that cost ratio.

ETDRK4 and the cancellation problem

ETDRK4 (Cox and Matthews, 2002) uses four RK4-like stages. For an autonomous $\mathcal{N}$ (Kassam and Trefethen, 2005, §2):

$$\begin{aligned} a_k &= e^{\mathcal{L}h/2}u_k + \mathcal{L}^{-1}(e^{\mathcal{L}h/2} - I)\,\mathcal{N}(u_k),\\ b_k &= e^{\mathcal{L}h/2}u_k + \mathcal{L}^{-1}(e^{\mathcal{L}h/2} - I)\,\mathcal{N}(a_k),\\ c_k &= e^{\mathcal{L}h/2}a_k + \mathcal{L}^{-1}(e^{\mathcal{L}h/2} - I)\,\big(2\mathcal{N}(b_k) - \mathcal{N}(u_k)\big),\\ u_{k+1} &= e^{\mathcal{L}h}u_k + f_1\,\mathcal{N}(u_k) + 2f_2\,\big(\mathcal{N}(a_k) + \mathcal{N}(b_k)\big) + f_3\,\mathcal{N}(c_k), \end{aligned}$$

where $a_k, b_k, c_k$ are stage values and $f_1, f_2, f_3$ are fixed functions of $h\mathcal{L}$, for example $f_2 = h^{-2}\mathcal{L}^{-3}[2 + \mathcal{L}h + e^{\mathcal{L}h}(-2 + \mathcal{L}h)]$. As $\mathcal{L} \to 0$ each $f_i \to h/6$, and ETDRK4 becomes classical RK4. The order is 4.

The trap. $(e^z - 1)/z$ and the $f_i$ lose all their digits to cancellation for small $\lvert z\rvert$: at $z = 10^{-10}$ the direct formula gives $1.000000082740371$ instead of $1.00000000005$ (only 7 correct digits). Kassam and Trefethen evaluate these functions as the mean over equally spaced points on a circle around $z$ in the complex plane (32 points), computed once before time stepping. By the Cauchy integral formula, this mean equals the value at the centre, and no subtraction of nearly equal numbers occurs. ETDRK4 fits best when $e^{h\mathcal{L}}$ is cheap, as for a diagonal $\mathcal{L}$ in Fourier space on periodic domains (lesson 08).

Answers to section 1.

1. A factor of 8. Halving $\Delta x$ doubles the number of grid points. Since $\Delta t \le \Delta x^2/2$, it also needs 4 times more steps. In $d$ dimensions the factor is $2^{d+2}$.
2. It explodes. The Jacobian is $\lambda = -1000$, so $z = h\lambda = -10$ and $\lvert R(z)\rvert = \lvert 1 - 10\rvert = 9$. Any difference between $y_k$ and the slow solution (even a rounding error) is multiplied by $-9$ every step, although the true solution is smooth. Backward Euler gives $R = 1/(1 + 10) \approx 0.09$ and works. This is stiffness.

4. Worked example #

(a) Measuring an order. Solve $y' = -y$, $y(0) = 1$ up to $T = 1$. The exact value is $e^{-1} = 0.367879$.

1. Forward Euler multiplies by $R(-h) = 1 - h$ per step. With $h = 0.1$: $y_{10} = 0.9^{10} = 0.348678$, error $0.019201$.
2. With $h = 0.05$: $y_{20} = 0.95^{20} = 0.358486$, error $0.009394$.
3. Ratio $0.019201/0.009394 = 2.044$, so $p \approx \log_2 2.044 = 1.03$. Order 1. ✓
4. RK4 with $h = 0.1$ multiplies by

So $y_{10} = 0.36787977$, error $3.33\cdot10^{-7}$.

5. RK4 with $h = 0.05$: error $2.00\cdot10^{-8}$. Ratio $16.7$, so $p \approx 4.06$. Order 4. ✓

(b) The heat equation on three points. $u_t = u_{xx}$ ($\nu = 1$) on $[0, 1]$ with $\Delta x = 0.25$: three interior points and $u = 0$ at both ends.

1. The matrix is $A = 16\,\operatorname{tridiag}(1, -2, 1)$, a $3\times 3$ matrix.
2. Its eigenvalues are $\lambda_k = -\tfrac{4}{\Delta x^2}\sin^2(k\pi\Delta x/2)$, here $\lambda_k = -64\sin^2(k\pi/8)$:

3. Stiffness ratio: $54.63/9.37 = 5.83$. On 99 interior points it is about 4,050; it grows like $1/\Delta x^2$.
4. Forward Euler is stable when $\Delta t\,\lvert\lambda_3\rvert \le 2$, that is $\Delta t \le 2/54.63 = 0.0366$. The von Neumann bound $\Delta x^2/2 = 0.03125$ is a little stricter, as explained in section 3.
5. Take $\Delta t = 0.04$. The fast mode has $g_3 = 1 - 0.04 \cdot 54.63 = -1.185$. After 20 steps it has grown by $1.185^{20} \approx 30$, while the true mode has decayed by $e^{-54.63 \cdot 0.8} \approx 10^{-19}$.
6. Same $\Delta t = 0.04$, implicit: $z = -2.185$. Backward Euler gives $1/(1 + 2.185) = 0.314$ (damped); Crank–Nicolson gives $(1 - 1.093)/(1 + 1.093) = -0.044$ (damped, with a sign flip).
7. Now a step 100 times above the explicit limit on a fine grid: $z \approx -200$. Backward Euler gives $0.005$, but Crank–Nicolson gives $-0.980$: the fast mode survives and oscillates. A-stability alone is not enough.

5. In code #

NumPy and SciPy check the claims of this lesson with closed forms and library solvers.

import numpy as np
from scipy.optimize import brentq

# Stability functions R(z) of four one-step methods on y' = lam*y, z = h*lam
R = {
    "forward Euler":  lambda z: 1 + z,
    "backward Euler": lambda z: 1 / (1 - z),
    "trapezoid (CN)": lambda z: (1 + z / 2) / (1 - z / 2),
    "RK4":            lambda z: 1 + z + z**2 / 2 + z**3 / 6 + z**4 / 24,
}
# Order: |R(z) - e^z| = O(z^{p+1}); halving z divides the error by 2^{p+1}
for name, Rz in R.items():
    e1, e2 = abs(Rz(0.01) - np.exp(0.01)), abs(Rz(0.005) - np.exp(0.005))
    print(f"{name:15s} p = {np.log2(e1 / e2) - 1:.2f}   |R(-1e6)| = {abs(Rz(-1e6)):.1e}")
# forward Euler    p = 1.00   |R(-1e6)| = 1.0e+06
# backward Euler   p = 1.01   |R(-1e6)| = 1.0e-06
# trapezoid (CN)   p = 2.01   |R(-1e6)| = 1.0e+00
# RK4              p = 3.99   |R(-1e6)| = 4.2e+22

# Where RK4 leaves its stability region on the real and the imaginary axis
x_real = brentq(lambda x: R["RK4"](x) - 1, -3, -2)
y_imag = brentq(lambda y: abs(R["RK4"](1j * y)) - 1 - 1e-15, 2.5, 3.0)
print(f"{x_real:.4f} {y_imag:.4f}")   # -2.7853 2.8284

Only backward Euler damps a very stiff mode. The heat equation on 99 interior points:

nu, n = 1.0, 99
dx = 1 / (n + 1)
k = np.arange(1, n + 1)
lam = -(4 * nu / dx**2) * np.sin(k * np.pi / (2 * (n + 1)))**2   # eigenvalues of nu*A
dt_max = dx**2 / (2 * nu)
print(f"{lam.min():.0f} {lam.max():.2f} {lam.min() / lam.max():.0f}")   # -39990 -9.87 4052

for factor in (0.95, 1.05):
    g = np.abs(1 + factor * dt_max * lam).max()        # forward Euler, worst mode
    print(f"{factor}: max|g| = {g:.4f}, after 400 steps x{g**400:.1e}")
# 0.95: max|g| = 0.9995, after 400 steps x8.3e-01
# 1.05: max|g| = 1.0995, after 400 steps x3.0e+16

z = 100 * dt_max * lam.min()                          # 100 times above the explicit limit
print(f"{z:.0f} {1 / (1 - z):.4f} {(1 + z / 2) / (1 - z / 2):.4f}")   # -200 0.0050 -0.9802

The stiff equation of section 1 with SciPy's adaptive solvers:

from scipy.integrate import solve_ivp

a = 1000.0
f = lambda t, y: -a * (y - np.cos(t))
exact = lambda t: (a**2 * np.cos(t) + a * np.sin(t) - a**2 * np.exp(-a * t)) / (a**2 + 1)
for method in ("RK45", "Radau", "BDF"):
    sol = solve_ivp(f, (0.0, 10.0), [0.0], method=method, rtol=1e-6, atol=1e-9)
    print(method, sol.t.size - 1, f"{abs(sol.y[0, -1] - exact(10.0)):.0e}")
# RK45 3587 5e-07
# Radau 81 2e-08
# BDF 209 2e-09

The explicit RK45 needs 3,587 steps (the exact counts can change a little between SciPy versions); Radau needs 81 and BDF 209, and both are more accurate. The average RK45 step is $2.8\cdot10^{-3}$, so $h\lambda \approx -2.8$: the adaptive controller sits at the edge of the stability interval of its scheme (about $[-3.3, 0]$). That is the stiffness symptom of section 3.

6. Where it shows up #

• Weather and climate models. The fastest waves set the explicit limit $\Delta t \lesssim \Delta x / c$: sound at about 340 m/s on a 9 km grid allows only about 26 s. Many global models therefore treat fast waves semi-implicitly (often with semi-Lagrangian advection) and take steps of several minutes. Vertical diffusion in a column is treated implicitly, and atmospheric chemistry uses linearly implicit (Rosenbrock) solvers.
• Hybrid ML–physics models. A learned tendency $f_\theta$ inside an explicit scheme is stable only if $h\lambda$ lies in the stability region for the eigenvalues $\lambda$ of its Jacobian $\partial f_\theta/\partial u$. Good offline skill does not guarantee this. Get the eigenvalues from an autodiff Jacobian (lesson 01) at states from a long run.
• Autoregressive emulators. A network that maps the state to the state 6 hours later is an explicit one-step map. Long rollouts stay bounded only if errors are not amplified step after step (lesson 16).
• Neural ODEs. When the learned dynamics become stiff, an adaptive explicit solver shrinks its step, and the number of function evaluations grows during training. Kinetic-energy and Jacobian penalties (Finlay et al., 2020) keep the dynamics cheap to solve (lesson 03).
• Diffusion samplers. The probability-flow ODE of a diffusion model is semilinear. DPM-Solver and DEIS are exponential integrators: they integrate the linear part exactly and approximate only the network term, so 10–20 steps can be enough (lesson 13).
• Adjoints. Integrating $u$ backward in time flips the sign of every eigenvalue: decaying modes become growing ones. This is why re-integrating a dissipative system backward fails over long horizons (lesson 03).
• Seismology. Explicit finite-difference and spectral-element codes for seismic waves run at a fixed Courant number: the time step is set by the smallest element and the fastest (P-wave) speed. Refining the mesh near the surface makes the whole simulation slower, which is why local time stepping exists.

7. Common mistakes #

• Mixing up local and global order. One step of an order-$p$ method has error $O(h^{p+1})$; after $T/h$ steps the error is $O(h^p)$.
• A wrong slope. Too large $h$, too small $h$ (rounding floor), a non-smooth solution, or an inaccurate reference solution all spoil the measured order.
• "Implicit is more accurate." Implicit buys stability, not accuracy. With a large step, backward Euler damps everything fast, and Crank–Nicolson leaves stiff modes ringing.
• "A-stable is enough." For strongly damped fast modes you want L-stability (or a few backward Euler steps at the start, as in Rannacher smoothing).
• "The CFL condition guarantees stability." It is only necessary. Centred differences with forward Euler for advection are unstable at every Courant number.
• Forgetting the dimension. In $d$ dimensions the explicit diffusion limit is $\Delta x^2/(2d\nu)$.
• Evaluating $(e^z - 1)/z$ naively in exponential integrators.

8. Check yourself #

1. The explicit midpoint method (RK2) has $R(z) = 1 + z + z^2/2$. Find its stability interval on the negative real axis. Is it stable for pure advection with centred differences?
   Answer

   $\lvert 1 + x + x^2/2\rvert \le 1$ for $-2 \le x \le 0$. On the imaginary axis $\lvert R(iy)\rvert^2 = 1 + y^4/4 > 1$ for $y \ne 0$, so it is unstable for every step: centred advection has purely imaginary eigenvalues. (RK4 has $\lvert R(iy)\rvert^2 = 1 - y^6/72 + y^8/576$, which is $\le 1$ up to $\lvert y\rvert = 2\sqrt2$, so RK4 works.)
2. A convergence study gives errors $3.2\cdot10^{-3}$, $2.0\cdot10^{-4}$, $1.25\cdot10^{-5}$ for $h$, $h/2$, $h/4$. What is the order?
   Answer

   Each ratio is 16, so $p = \log_2 16 = 4$ (for example RK4).
3. Derive the forward Euler limit for the 2-D heat equation $u_t = \nu(u_{xx} + u_{yy})$ with the 5-point stencil and $\Delta x = \Delta y$.
   Answer

   With $u_{jl}^m = g^m e^{i(\kappa x_j + \eta y_l)}$: $g = 1 - 4r[\sin^2(\kappa\Delta x/2) + \sin^2(\eta\Delta x/2)]$, $r = \nu\Delta t/\Delta x^2$. The worst mode has both $\sin^2 = 1$: $1 - 8r \ge -1 \iff r \le 1/4$, so $\Delta t \le \Delta x^2/(4\nu)$.
4. Why is Crank–Nicolson a poor choice for a diffusion problem stepped far above the explicit limit, although it is A-stable and second order?
   Answer

   $R(z) \to -1$ as $z \to -\infty$: the stiff grid-scale modes are hardly damped and flip sign every step. This gives spurious oscillations, especially after a rough initial condition. An L-stable method (backward Euler, BDF2, Radau IIA), or a few backward Euler start-up steps, damps them.
5. A wind of 50 m/s blows over a grid with $\Delta x = 25$ km. What is the largest stable step for explicit upwind advection, and why is it a limit of the CFL type?
   Answer

   $\mathrm{Co} = c\Delta t/\Delta x \le 1$ gives $\Delta t \le 25{,}000/50 = 500$ s. In one step the air moves at most one grid cell, so the upwind stencil still contains the departure point (the domain of dependence).
6. For $u' = \mathcal{L}u + \mathcal{N}(u)$ with $\mathcal{L}$ diagonal in Fourier space and very stiff, what do ETDRK4 and IMEX Euler each remove, and what limit remains?
   Answer

   Both remove the step limit from $\mathcal{L}$: ETDRK4 integrates the linear part exactly (order 4), IMEX Euler treats it implicitly (order 1). The explicit treatment of $\mathcal{N}$ keeps its own limit, for example an advective CFL condition. ETDRK4 also needs a stable evaluation of its $\varphi$-type coefficients (the contour mean).

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Define the order $p$ of a one-step method.

   Show answer

2. Give the stability functions $R(z)$ of forward Euler, backward Euler and the trapezoid rule.

   Show answer

3. Forward Euler on $y' = \lambda y$ with real $\lambda < 0$: what is the largest stable step?

   Show answer

4. What is the difference between A-stable and L-stable?

   Show answer

5. What is the explicit stability limit for the heat equation $u_t = \nu u_{xx}$ with centred differences in 1-D?

   Show answer

6. State the CFL condition for upwind advection, and say why it holds.

   Show answer

7. Define stiffness.

   Show answer

8. How do local and global error relate for an order-$p$ method?

   Show answer

9. What does an IMEX method split, and why?

   Show answer

10. Why does an explicit method always have a bounded stability region?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What is stiffness, and how do you detect it?" — Stability, not accuracy, sets the step. Look at the eigenvalues of the Jacobian, or at an adaptive explicit solver that takes many tiny steps on a smooth solution.
• "Explicit or implicit for the heat equation on a fine grid?" — Implicit (or IMEX, or ETD): the explicit step $\propto \Delta x^2$ is far below what accuracy needs; an implicit step in 1-D is a tridiagonal solve, $O(n)$.
• "What does the CFL condition say, and is it sufficient?" — In one step, information may not travel further than the stencil reaches; necessary, not sufficient.
• "How do you verify a new solver?" — A convergence study (log–log slope = theoretical order), exact solutions, conservation checks, stability tests at the limit.
• "Why can a learned or hybrid weather model blow up after many steps?" — Error amplification above 1 per step, or eigenvalues of $h\,\partial f_\theta/\partial u$ outside the stability region. Separate integrator problems from model problems.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It has five exercises, each with a check cell and a folded solution:

1. Write an RK4 step, solve $y' = -y$, and measure the order from the errors at several step sizes.
2. Draw the stability region of a method from its $R(z)$ (a contour plot of $\lvert R(z)\rvert = 1$) and find where it crosses the negative real axis.
3. Run the explicit heat equation, find the blow-up threshold numerically, and compare it with $\Delta x^2/(2\nu)$.
4. Stiff problem: compare forward and backward Euler by the number of steps they need for the same accuracy.
5. Evaluate $\varphi_1(z) = (e^z - 1)/z$ near $z = 0$ by the direct formula and by the contour mean, and count the correct digits.

11. Watch and read #

Stability of Forward Euler and Backward Euler Integration Schemes for Differential EquationsSteve Brunton · 33 min

The test equation and the stability regions of forward and backward Euler, as in section 2.

Runge-Kutta Integrator Overview: All Purpose Numerical Integration of Differential EquationsSteve Brunton · 30 min

How RK2 and RK4 combine several slopes per step, and why that raises the order.

Stability function of a Runge-Kutta methodMelvin Leok · 8 min

How to compute R(z) for any Runge–Kutta method.

Stiff ODEsMelvin Leok · 20 min

Stiffness, A-stability and why implicit methods are needed.

Von Neumann analysis: derivation of the CFL criterionLMU Seismology · 20 min

Von Neumann analysis of a wave-equation scheme, from a seismology course.

• Driscoll and Braun, Fundamentals of Numerical Computation (free online): §6.2 (local truncation error, consistency, convergence), §6.4 (Runge–Kutta methods), §6.8 (zero-stability and the Dahlquist equivalence theorem), §11.2–11.3 (method of lines, absolute stability, the heat equation), §11.4 (stiffness), §12.2 (upwinding and the CFL condition).
• Kassam and Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput. 26 (2005): §2 (ETDRK4 and the contour integral), §3 (comparison with IMEX and other methods).
• exponax: ETDRK solvers of orders 1–4 for periodic semilinear PDEs in JAX, used to generate training data for neural PDE emulators.