Numerical solution of differential equations

Integrators for ODEs and discretizations of PDEs.

In this chapter 9

Tick a unit to add it to your PDF basket.

Ordinary differential equations
Partial differential equations

Ordinary differential equations

Method

Integrators for ODEs

Adapted

Add to PDF

Builds the one-step and multistep integrators for $y' = f(x,y)$ from the Taylor expansion they imitate, defines the local truncation error and the order, proves that a method of order $p$ whose increment is Lipschitz has global error at most $(C/L)(e^{L(b-a)}-1)h^{p}$ , and verifies the orders $1$ , $2$ , $2$ and $4$ of Euler, Heun, Adams-Bashforth and classical Runge-Kutta on a nonlinear problem with a known solution.

Prerequisites: Existence and uniqueness of solutions, Taylor's theorem, Dependence on initial data and Gronwall's inequality, Linear systems and the matrix exponential

Definition(The problem, the grid, and a one-step method)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6, Eq. (27.60), p. 1020; Hirsch, Smale & Devaney (2013), §7.5, p. 154

Throughout, the object to be solved is the initial-value problem

y' = f(x, y), \qquad y(a) = y_a ,

on an interval $[a,b]$ , where $y$ may be a vector in $\mathbb{R}^{m}$ and $f$ is defined on a region $R \subseteq [a,b]\times\mathbb{R}^{m}$ . By Existence and uniqueness of solutions a unique solution exists as long as $f$ is continuous and Lipschitz in $y$ , and that is the function the methods below try to produce. Nothing is lost by restricting attention to first-order systems: an equation $y^{(R)} = F(x, y, y', \dots, y^{(R-1)})$ becomes the system $y' = p_1,\ p_1' = p_2,\ \dots,\ p_{R-1}' = F$ in the $R$ unknowns $y, p_1, \dots, p_{R-1}$ , and every method here is applied componentwise, all components being advanced through a step before any is advanced through the next. Riley, Hobson and Bence set this reduction out in full; it is the first thing done to a high-order equation before it is handed to software, and it is the reason the units on linear systems (Linear systems and the matrix exponential) and on integrators speak the same language.

Fix a step size $h > 0$ and the grid $x_n = a + nh$ . The methods produce numbers $y_n$ intended to approximate $y(x_n)$ ; $y_n$ is a vector when $y$ is.

A one-step method is a rule

y_{n+1} \;=\; y_n + h\,\Phi(x_n, y_n; h) ,

in which the increment function $\Phi$ is built from evaluations of $f$ alone. The method is explicit when $\Phi$ can be evaluated from $(x_n, y_n, h)$ by a finite list of substitutions, as all the methods in this unit can; methods for which $y_{n+1}$ appears on both sides are implicit and are the subject of Stability and stiff equations.

A linear multistep method with $k$ steps is a rule

y_{n+1} \;=\; \sum_{j=0}^{k-1}\alpha_j\,y_{n-j} \;+\; h\sum_{j=-1}^{k-1}\beta_j\,f(x_{n-j},\,y_{n-j}) ,

with constant coefficients; it is explicit when $\beta_{-1} = 0$ . Such a method needs $k$ starting values, and only one of them is given by the problem.

Definition(Local truncation error, order, convergence)Original

Let $y$ denote the exact solution. The local truncation error of a one-step method at $(x, y(x))$ is the amount by which the exact solution fails to satisfy the method's own recurrence over one step:

\tau(x, h) \;=\; y(x+h) \;-\; y(x) \;-\; h\,\Phi\bigl(x, y(x); h\bigr) .

It is the error a single step would commit if it started from exact data. The method has order $p$ on $[a,b]$ if there is a constant $C$ , independent of $x$ and $h$ , with

\lVert\tau(x,h)\rVert \;\le\; C\,h^{\,p+1} \qquad\text{for all } x \in [a,b] \text{ and all small } h ,

and if $p$ is the largest exponent for which this holds. A method of order $p \ge 1$ is consistent: dividing by $h$ , the recurrence reproduces the differential equation in the limit.

The global error is $e_n = y(x_n) - y_n$ , the quantity the user actually suffers. The method converges if $\max_{n}\lVert e_n\rVert \to 0$ as $h \to 0$ with $b - a$ fixed. The one power of $h$ that separates $\tau$ from the global error is the bookkeeping of the theorem below: roughly $(b-a)/h$ local errors of size $h^{p+1}$ accumulate into a global error of size $h^{p}$ .

For a multistep method the local truncation error is defined by the same substitution of the exact solution into the recurrence, and the order is again the largest $p$ with $\lVert\tau\rVert \le Ch^{p+1}$ .

Definition(The five methods used here)AdaptedAdapted from Strogatz (1994), §2.8, pp. 32-34; Hirsch, Smale & Devaney (2013), §7.5, pp. 154-155; Riley, Hobson & Bence (2006), §27.6, Eqs. (27.63), (27.73), (27.78), (27.80), pp. 1021, 1025, 1027, 1028

Write $f_n = f(x_n, y_n)$ . Each rule below says which slope is used to cross the interval $[x_n, x_{n+1}]$ ; Hirsch, Smale and Devaney's way of putting it is that a method is nothing but a choice of the straight line along which the step is taken.

Euler's method uses the slope at the left endpoint:

y_{n+1} = y_n + h f_n .

Strogatz reads it off the vector field directly: travel for time $h$ at the velocity found where you are.

The improved Euler (Heun) method takes a trial step with Euler, evaluates the slope there, and crosses with the average of the two slopes:

\tilde y = y_n + h f_n , \qquad y_{n+1} = y_n + \tfrac{h}{2}\Bigl( f_n + f(x_{n+1}, \tilde y) \Bigr) .

The midpoint method instead uses a half-step probe and crosses with the slope found there:

y_{n+1} = y_n + h\,f\Bigl(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2} f_n\Bigr) .

The classical fourth-order Runge-Kutta method (RK4) uses four slopes:

\begin{aligned} k_1 &= f(x_n,\, y_n), &\qquad k_2 &= f\bigl(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2}k_1\bigr), \\ k_3 &= f\bigl(x_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2}k_2\bigr), &\qquad k_4 &= f\bigl(x_n + h,\; y_n + h k_3\bigr), \end{aligned}

y_{n+1} = y_n + \tfrac{h}{6}\bigl(k_1 + 2k_2 + 2k_3 + k_4\bigr) .

The two-step Adams-Bashforth method (AB2) is not a one-step method: it reuses the slope already computed at the previous grid point,

y_{n+1} = y_n + \tfrac{h}{2}\bigl(3 f_n - f_{n-1}\bigr) .

Riley, Hobson and Bence derive it by fitting a straight line to $f$ through the two points $(x_{n-1}, f_{n-1})$ and $(x_n, f_n)$ and integrating that line across $[x_n, x_{n+1}]$ ; the extra term $\tfrac{h}{2}(f_n - f_{n-1})$ is a correction to the Euler step. It costs one evaluation of $f$ per step, against two for Heun and the midpoint rule and four for RK4, because $f_{n-1}$ is kept from the previous step.

Method(Convergence of a one-step method, and the orders of the five methods)Original

1 (Convergence). Let $R = [a,b]\times U$ with $U \subseteq \mathbb{R}^{m}$ , let the exact solution $y$ of $y' = f(x,y)$ , $y(a) = y_a$ stay in $R$ on $[a,b]$ , and consider a one-step method started exactly, $y_0 = y_a$ , whose iterates $y_n$ also stay in $R$ . Suppose

(L) the increment function satisfies $\lVert \Phi(x, u; h) - \Phi(x, v; h)\rVert \le L\,\lVert u - v\rVert$ for all $(x,u), (x,v) \in R$ and all $h \in (0, h_0]$ , with $L > 0$ ; and
(T) the local truncation error obeys $\lVert\tau(x,h)\rVert \le C h^{p+1}$ on $[a,b]$ for $h \in (0, h_0]$ .

Then for every $h \in (0,h_0]$ and every grid point $x_n \in [a,b]$ ,

\lVert y(x_n) - y_n \rVert \;\le\; \frac{C}{L}\Bigl(e^{L(x_n - a)} - 1\Bigr)\,h^{p} .

In particular the method converges, with global error $O(h^{p})$ : one power of $h$ is lost relative to the local error.

2 (Orders). Let $f$ be as smooth as each statement requires on a neighbourhood of the solution. Then

Euler's method has order $1$ , with local truncation error $\tfrac{1}{2}h^{2}y''(x) + O(h^{3})$ ;
the Heun and midpoint methods have order $2$ ;
AB2 has order $2$ , with local truncation error $\tfrac{5}{12}h^{3}y'''(x) + O(h^{4})$ ;
RK4 has order $4$ .

3 (The second-order family). Every method of the two-stage form

y_{n+1} = y_n + h\Bigl[\alpha_1 f_n + \alpha_2 f\bigl(x_n + \beta h,\, y_n + \beta h f_n\bigr)\Bigr]

has order at least $2$ if and only if

\alpha_1 + \alpha_2 = 1 \qquad\text{and}\qquad \alpha_2\beta = \tfrac12 ,

so the family has one free parameter. Heun is $\alpha_1 = \alpha_2 = \tfrac12$ , $\beta = 1$ ; the midpoint method is $\alpha_1 = 0$ , $\alpha_2 = 1$ , $\beta = \tfrac12$ . No member of the family has order $3$ for all $f$ .

Idea of the proofOriginal

Every method here is an attempt to imitate the Taylor expansion of the exact solution without differentiating $f$ . The expansion $y(x+h) = y + hy' + \tfrac12 h^{2}y'' + \cdots$ needs $y'' = f_x + f_y f$ , $y'''$ , and so on; computing those by hand is exactly what Runge-Kutta methods avoid, by evaluating $f$ itself at cleverly placed points and combining the values so that the combination reproduces the Taylor coefficients. Matching through $h^{p}$ is what "order $p$ " means, and part 3 shows the matching conditions in the smallest nontrivial case.

The convergence statement is an accumulation argument with one twist. Over a step, the previous error is carried forward and distorted by the method, and a fresh local error $\tau$ is added. Hypothesis (L) says the distortion is at worst a factor $1 + hL$ , so the errors obey $\varepsilon_{n+1} \le (1+hL)\varepsilon_n + Ch^{p+1}$ . The twist is that summing this geometric recursion produces $\bigl((1+hL)^{n} - 1\bigr)/(hL)$ , and the factor $hL$ in the denominator cancels one power of $h$ from the numerator $Ch^{p+1}$ ; that is where the lost power goes. Finally $(1+hL)^{n} \le e^{nhL}$ converts the discrete growth factor into the exponential that Gronwall's inequality produces for the continuous problem. The exponential is not an artefact: two exact solutions of the same equation really can separate that fast, so no method can do better.

Proof(Part 1)Original

Step 1 (one step of the recursion). Write $e_n = y(x_n) - y_n$ and $\varepsilon_n = \lVert e_n\rVert$ . By the definition of the local truncation error applied at $x_n$ ,

y(x_{n+1}) = y(x_n) + h\,\Phi\bigl(x_n, y(x_n); h\bigr) + \tau(x_n, h) ,

while the method gives $y_{n+1} = y_n + h\Phi(x_n, y_n; h)$ . Subtracting,

e_{n+1} = e_n + h\Bigl[\Phi\bigl(x_n, y(x_n); h\bigr) - \Phi(x_n, y_n; h)\Bigr] + \tau(x_n, h) .

Step 2 (estimate it). Both $(x_n, y(x_n))$ and $(x_n, y_n)$ lie in $R$ , by the hypothesis that the exact solution and the iterates stay there. So hypothesis (L) bounds the bracket by $L\varepsilon_n$ , and hypothesis (T) bounds $\tau$ by $Ch^{p+1}$ . The triangle inequality gives

\varepsilon_{n+1} \;\le\; (1 + hL)\,\varepsilon_n \;+\; C h^{p+1} . \tag{$*$}

Step 3 (an accumulation lemma). If $\varepsilon_0 = 0$ and $(*)$ holds with $A = 1+hL > 1$ and $B = Ch^{p+1} \ge 0$ , then $\varepsilon_n \le B\,(A^{n}-1)/(A-1)$ . Induct on $n$ . The case $n = 0$ reads $\varepsilon_0 \le 0$ , which is the hypothesis. If the bound holds at $n$ then

\varepsilon_{n+1} \le A\cdot B\frac{A^{n}-1}{A-1} + B = B\,\frac{A^{n+1} - A + A - 1}{A-1} = B\,\frac{A^{n+1}-1}{A-1} ,

which is the bound at $n+1$ .

Step 4 (substitute and simplify). Here $A - 1 = hL$ , so

\varepsilon_n \;\le\; \frac{C h^{p+1}}{hL}\Bigl((1+hL)^{n} - 1\Bigr) = \frac{C}{L}\,h^{p}\Bigl((1+hL)^{n} - 1\Bigr) .

The single power of $h$ cancelled against $hL$ is the one lost between local and global order.

Step 5 (exponential form). For $t \ge 0$ the convexity of the exponential gives $1 + t \le e^{t}$ , so $(1+hL)^{n} \le e^{nhL} = e^{L(x_n - a)}$ because $nh = x_n - a$ . Substituting into Step 4 gives the claim. The starting value is exact, $\varepsilon_0 = 0$ , as assumed in Step 3; if instead $\varepsilon_0 > 0$ the same induction yields the extra term $\varepsilon_0 e^{L(x_n-a)}$ . $\square$

Proof(Part 2, except RK4)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6, Eqs. (27.68)-(27.70) and (27.73), pp. 1023-1025

Throughout, $y$ is the exact solution through $(x, y(x))$ , so that $y' = f(x,y(x))$ , and all derivatives are evaluated at $x$ unless marked otherwise. We use the chain rule in the form Riley, Hobson and Bence record, $y'' = f_x + f_y\,f$ , and Taylor's theorem (Taylor's theorem).

Euler. With $\Phi = f(x, y(x))= y'$ ,

\tau = y(x+h) - y(x) - h y' = \tfrac12 h^{2} y''(x) + O(h^{3}) ,

which is $O(h^{2})$ and no better unless $y'' \equiv 0$ . So $p = 1$ .

The two-stage family, and Heun and midpoint with it. Expand the second stage of the family in part 3 of the statement by Taylor's theorem in two variables, about $(x, y(x))$ :

f\bigl(x + \beta h,\; y + \beta h f\bigr) = f + \beta h\bigl(f_x + f_y f\bigr) + \tfrac12\beta^{2}h^{2}\bigl(f_{xx} + 2f f_{xy} + f^{2}f_{yy}\bigr) + O(h^{3}) .

Hence the increment is

h\Phi = (\alpha_1+\alpha_2)h f + \alpha_2\beta h^{2}\bigl(f_x + f_y f\bigr) + O(h^{3}) = (\alpha_1+\alpha_2)h y' + \alpha_2\beta h^{2} y'' + O(h^{3}) .

The exact increment is $y(x+h) - y(x) = h y' + \tfrac12 h^{2}y'' + \tfrac16 h^{3}y''' + O(h^{4})$ . Matching the $h$ and $h^{2}$ terms gives exactly the two conditions $\alpha_1 + \alpha_2 = 1$ and $\alpha_2\beta = \tfrac12$ , and these are necessary as well as sufficient because $f$ and the combination $f_x + f_y f$ can be prescribed independently at a point. Heun and the midpoint method satisfy them, so both have order at least $2$ .

For the last sentence of part 3, carry the expansion one term further. The $h^{3}$ coefficient produced by the family is

\tfrac12\alpha_2\beta^{2}\bigl(f_{xx} + 2ff_{xy} + f^{2}f_{yy}\bigr) ,

while the exact solution has $\tfrac16 y''' = \tfrac16\bigl(f_{xx} + 2ff_{xy} + f^{2}f_{yy}\bigr) + \tfrac16 f_y\bigl(f_x + f_y f\bigr)$ . The second group, $\tfrac16 f_y(f_x+f_yf)$ , cannot be produced at all: no choice of $\alpha_2, \beta$ puts a term containing $f_y$ into the $h^{3}$ coefficient. So the family never reaches order $3$ , and two stages are genuinely not enough. (This is where the extra stages of RK4 earn their keep.)

AB2. Substituting the exact solution, and writing $y'(x-h) = y' - hy'' + \tfrac12 h^{2}y''' + O(h^{3})$ ,

\tfrac{h}{2}\bigl(3y'(x) - y'(x-h)\bigr) = h y' + \tfrac12 h^{2}y'' - \tfrac14 h^{3}y''' + O(h^{4}) .

Subtracting this from $y(x+h) - y(x) = hy' + \tfrac12 h^{2}y'' + \tfrac16 h^{3}y''' + O(h^{4})$ , the first two terms cancel and

\tau = \Bigl(\tfrac16 + \tfrac14\Bigr)h^{3}y''' + O(h^{4}) = \tfrac{5}{12}h^{3}y''' + O(h^{4}) ,

so AB2 has order $2$ with the error constant $5/12$ . $\square$

Remark(RK4, and why its order is asserted rather than derived here)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6, Eq. (27.80), p. 1028; Strogatz (1994), §2.8, pp. 33-34

Matching the Taylor expansion through $h^{4}$ imposes the eight order conditions through order $4$ on the coefficients of a four-stage Runge-Kutta method, and verifying them by hand is several pages of bookkeeping with no idea in it. Riley, Hobson and Bence print the scheme "without derivation" for the same reason, and Strogatz calls the four numbers cunningly chosen and leaves their choice to an exercise. We follow them in stating the order rather than deriving it, and verify it in two independent ways instead.

Symbolically. Expanding $y_{n+1}$ and the exact solution in powers of $h$ with computer algebra, for the non-autonomous field

f(x,y) = \tfrac23 + \tfrac57 x + \tfrac{3}{11}y + \tfrac15 x^{2} + \tfrac{7}{13}xy + \tfrac29 y^{2} + \tfrac{1}{17}x^{3} + \tfrac{4}{19}xy^{2} + \tfrac{1}{23}y^{3} + \tfrac{3}{29}x^{2}y + \tfrac{1}{31}y^{4} + \tfrac{1}{37}x^{4}

at $(x,y) = (\tfrac12, \tfrac13)$ , the two series agree through $h^{4}$ and first differ at $h^{5}$ , with coefficient $0.0250846$ . The same computation returns leading powers $h^{2}$ for Euler, $h^{3}$ for Heun and for the midpoint rule, and $h^{4}$ for the third-order scheme Riley, Hobson and Bence give as their Eq. (27.79) — that is, orders $1$ , $2$ , $2$ and $3$ , confirming the hand derivations above and showing that this particular $f$ is generic enough to detect the true order.

Numerically. The error table in the first example below halves the step four times and shows the RK4 error falling by factors tending to $16$ , which is the signature of order $4$ .

Remark(Where the hypotheses are used)Original

(L), the Lipschitz condition on $\Phi$ , is what makes errors merely grow exponentially instead of arbitrarily. It is inherited from $f$ : for Euler, $\Phi = f$ , so (L) is the Lipschitz condition of Existence and uniqueness of solutions itself. For a multi-stage method each stage is a Lipschitz function of the previous one, so if $f$ has Lipschitz constant $K$ in $y$ on $R$ then $\Phi$ has a Lipschitz constant of the form $L = K(1 + c_1 hK + c_2 h^{2}K^{2} + \cdots)$ , bounded by a constant multiple of $K$ for $h \le h_0$ . Without (L) the recursion in Step 2 has nothing to bound the propagated error by: the first stage of the proof fails immediately, not just the constant.

(T) must be uniform in $x$ . The constant $C$ involves derivatives of the solution, through $y''$ , $y'''$ and so on. On an interval where the solution approaches a singularity these blow up, and the bound is true but empty. The example below is deliberately set on $[0, \tfrac12]$ for this reason: its exact solution $1/(1-\sin x)$ has a pole at $x = \pi/2$ , and the error constants degrade badly as that point is approached.

The iterates must stay in $R$ . This is a hypothesis, not a conclusion, and it is the one that is usually verified after the fact. It can genuinely fail: for $y' = y^{2}$ , $y(0) = 1$ , whose solution blows up at $x = 1$ , an Euler step from near the pole lands far outside any region on which a Lipschitz constant was computed. In practice one integrates, checks that the computed track stayed inside the region used to compute $L$ , and reruns with a smaller step if not.

Order is a statement about smooth $f$ . Each order calculation differentiates $f$ as many times as the expansion requires: order $p$ needs $f \in C^{p}$ near the solution. Applied to an $f$ with a kink — a switch in a control law, a table lookup, a friction law that changes form at zero velocity — RK4 is not fourth order across the kink, and refining the step buys much less than the table of ratios would suggest. The standard remedy is to stop the integration at the discontinuity and restart.

The exponential factor is sharp. Taking $f(x,y) = Ly$ makes $\Phi$ exactly $L$ -Lipschitz, and two exact solutions separate by precisely $e^{L(x-a)}$ times their initial gap. So no sharper general bound than part 1 exists, and the $e^{L(b-a)}$ is what makes long integrations of systems with positive Lyapunov exponents (Lyapunov exponents and sensitive dependence) untrustworthy trajectory by trajectory.

Round-off sets a floor that the theorem does not see. Part 1 says the error decreases like $h^{p}$ forever, which is false in floating-point arithmetic: each of the $(b-a)/h$ steps commits a rounding error of relative size about the unit roundoff $u$ (Floating-point arithmetic and error), and these accumulate roughly as $u(b-a)/h$ , growing as $h$ shrinks. The total error is therefore U-shaped in $h$ , with a minimum near $h \sim u^{1/(p+1)}$ . Strogatz makes exactly this point when he warns against simply choosing a tiny step and leaving it there.

Method(Running an integrator on a real problem)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6.3, p. 1025 and §27.7, Eqs. (27.82)-(27.84), pp. 1028-1029

Reduce to a first-order system. Introduce $p_r = y^{(r)}$ as in the first definition, so the problem becomes $Y' = F(x, Y)$ with $Y \in \mathbb{R}^{m}$ . Advance all components through each step before starting the next.
Non-dimensionalise. Scale $x$ and each component of $Y$ so that the quantities of interest are of order one. Step-size choices, error tolerances and the convergence table all become meaningful only after this.
Estimate the two constants. Compute $K = \sup\lVert \partial F/\partial Y\rVert$ over the region you expect to visit; this gives the Lipschitz constant and the factor $e^{K(b-a)}$ that no method can avoid. If $e^{K(b-a)}$ is astronomically large, no accuracy target is achievable and the question being asked of the computation should change.
Choose the method by the cost of $F$ . Count evaluations per step: $1$ for Euler and AB2, $2$ for Heun and midpoint, $4$ for RK4. When $F$ is expensive and the solution is smooth, the higher order wins easily, because reaching a tolerance $\eta$ needs $h \sim (\eta)^{1/p}$ and so $\sim \eta^{-1/p}$ steps.
Start a multistep method with a one-step method of at least the same order. AB2 needs $y_1$ as well as $y_0$ ; take it from one RK4 step, whose $O(h^{5})$ error cannot pollute an $O(h^{2})$ computation. Riley, Hobson and Bence note the same need for the central-difference and Adams schemes and suggest a Taylor expansion; a Runge-Kutta step is the modern answer and needs no derivatives of $f$ .
Measure the order, do not assume it. Run at $h$ , $h/2$ , $h/4$ and tabulate the differences of the answers. Successive differences falling by $2^{p}$ confirm the order and simultaneously confirm that the solution is smooth enough for it.
Control the step by step doubling. Take one step of size $h$ and two of size $h/2$ from the same point. For a method of order $p$ the difference $\Delta$ of the two answers is about $2^{p}-1$ times the error of the finer one — equivalently $(2^{p}-1)/2^{p}$ times the error of the coarser one — so the finer answer carries an error of about $\Delta/(2^{p}-1)$ , which is $\Delta/15$ for RK4; that is an error estimate for free; accept the step if it is below tolerance, and set the next step to $h\,(\eta/\text{estimate})^{1/(p+1)}$ with a safety factor of about $0.9$ . Richardson extrapolation can also be used to improve the answer, at the price of losing the error estimate.
Predict and correct, if the method is implicit. Riley, Hobson and Bence's pattern — predict $y_{n+1}$ with an explicit formula, evaluate $f$ there, then recompute $y_{n+1}$ with the better formula — is the cheap way to use an implicit formula without solving an equation. Iterating the correction converges only when $h$ times the Lipschitz constant is small, which is precisely what fails for the stiff problems of Stability and stiff equations; there the corrector must be solved properly, by Newton's method (Root finding).
Check a conserved quantity if the problem has one. A drift in the energy of a mechanical system is the clearest diagnostic there is, and for long integrations it is a reason to change the class of method entirely (Symplectic integrators).

Example(Four methods on a nonlinear problem with a known solution)Original

Take

y' = y^{2}\cos x , \qquad y(0) = 1 ,

whose solution is found by separation: $-1/y = \sin x + c$ with $c = -1$ , so

y(x) = \frac{1}{1 - \sin x} ,

valid up to the pole at $x = \pi/2$ . We integrate to $b = \tfrac12$ , where $y(\tfrac12) = 1/(1 - \sin\tfrac12) = 1.920954780068805$ . All arithmetic below was done in $40$ -digit precision so that nothing reported is round-off.

One step by hand, $h = 0.1$ . Here $f(0,1) = 1$ , so Euler gives $y_1 = 1.1$ . For RK4,

k_1 = 1, \quad k_2 = f(0.05,\, 1.05) = 1.101122162, \quad k_3 = f(0.05,\, 1.0550561) = 1.111752252, \quad k_4 = f(0.1,\, 1.1111752) = 1.228541972 ,

and $y_1 = 1 + \tfrac{0.1}{6}(k_1 + 2k_2 + 2k_3 + k_4) = 1.110904846666$ . The exact value is $y(0.1) = 1.110905490709$ . The three errors after one step are

\text{Euler } -1.0905\times10^{-2}, \qquad \text{Heun } -7.0774\times10^{-4}, \qquad \text{RK4 } -6.4404\times10^{-7} ,

in the ratio $1 : 6.49\times10^{-2} : 5.91\times10^{-5}$ — already the $h^{2}$ , $h^{3}$ , $h^{5}$ hierarchy of local errors at a single step size.

Global errors at $x = \tfrac12$ . Each column gives the error $y_N - y(\tfrac12)$ and, beneath, the factor by which it shrinks when $h$ is halved. AB2 was started with one RK4 step, as step 5 of the method recommends; the column depends on that choice, and a cruder start changes its leading digits.

$h$	Euler	Heun	midpoint	AB2	RK4
$0.1$	$-1.56455\times10^{-1}$	$-1.43004\times10^{-2}$	$-1.51474\times10^{-2}$	$-3.18433\times10^{-2}$	$-2.57489\times10^{-5}$
$0.05$	$-8.88446\times10^{-2}$	$-3.88394\times10^{-3}$	$-4.24395\times10^{-3}$	$-1.05151\times10^{-2}$	$-1.70535\times10^{-6}$
$0.025$	$-4.77691\times10^{-2}$	$-1.00876\times10^{-3}$	$-1.12144\times10^{-3}$	$-2.97218\times10^{-3}$	$-1.09364\times10^{-7}$
$0.0125$	$-2.48381\times10^{-2}$	$-2.56773\times10^{-4}$	$-2.88046\times10^{-4}$	$-7.85310\times10^{-4}$	$-6.91714\times10^{-9}$
$0.00625$	$-1.26747\times10^{-2}$	$-6.47541\times10^{-5}$	$-7.29769\times10^{-5}$	$-2.01464\times10^{-4}$	$-4.34786\times10^{-10}$
$0.003125$	$-6.40362\times10^{-3}$	$-1.62578\times10^{-5}$	$-1.83651\times10^{-5}$	$-5.09951\times10^{-5}$	$-2.72496\times10^{-11}$
ratios	$1.761,\,1.860,\,1.923,\,1.960,\,1.979$	$3.682,\,3.850,\,3.929,\,3.965,\,3.983$	$3.569,\,3.784,\,3.893,\,3.947,\,3.974$	$3.028,\,3.538,\,3.785,\,3.898,\,3.951$	$15.10,\,15.59,\,15.81,\,15.91,\,15.96$

The five columns converge to the factors $2$ , $4$ , $4$ , $4$ and $16$ , confirming orders $1$ , $2$ , $2$ , $2$ and $4$ . Heun and the midpoint rule, the two members of the family of part 3 that we singled out, cost the same and differ only in their error constants: $1.665$ against $1.881$ , so on this problem Heun's error is smaller by about $11\%$ at a fixed step. Since the error goes as $h^{2}$ , that buys only a $6\%$ larger step at a given tolerance, which is no reason to prefer one over the other. AB2 approaches its asymptotic factor more slowly than Heun, as a two-step method with a larger error constant should: at $h = 0.1$ its error is $2.2$ times Heun's, although both are second order.

Cost at a fixed accuracy. Read an error constant off the finest row, $|e| \approx C h^{p}$ : it is $C = 2.0492$ for Euler, $1.6648$ for Heun and $0.28573$ for RK4. Solving $Ch^{p} = 10^{-5}$ and dividing $\tfrac12$ by the result gives

method	$h$ for $10^{-5}$	steps	evaluations of $f$
Euler	$4.88\times10^{-6}$	$102\,458$	$102\,458$
Heun	$2.451\times10^{-3}$	$205$	$410$
RK4	$7.691\times10^{-2}$	$7$	$28$

So RK4 costs about $3660$ times fewer evaluations than Euler here and about $15$ times fewer than Heun, and both advantages grow without limit as the tolerance tightens: the evaluation count scales as $\eta^{-1}$ for Euler, $\eta^{-1/2}$ for Heun and $\eta^{-1/4}$ for RK4.

The convergence bound, checked. Take the region $R = [0,\tfrac12]\times[0,2]$ , which contains both the exact solution and every Euler iterate computed above (the largest over all six runs is $1.9146$ , at the finest step; at $h=0.1$ it is $1.7645$ ). On $R$ , $|\partial f/\partial y| = |2y\cos x| \le 4$ , so $L = 4$ serves for Euler, whose $\Phi$ is $f$ . The local error constant is $C = \tfrac12\max|y''|$ , and $y'' = 2y^{3}\cos^{2}x - y^{2}\sin x$ increases on $[0,\tfrac12]$ from $2$ to $9.149247$ , so $C = 4.574624$ . Part 1 then predicts

|y(\tfrac12) - y_N| \;\le\; \frac{C}{L}\bigl(e^{L/2} - 1\bigr)h \;=\; \frac{4.574624}{4}\bigl(e^{2}-1\bigr)h \;=\; 7.30688\,h .

At $h = 0.1$ this is $0.730688$ against the measured $0.156455$ , at $h = 0.05$ it is $0.365344$ against $0.088845$ , and at $h = 0.025$ it is $0.182672$ against $0.047769$ . The bound holds with a factor of between $3.8$ and $4.7$ to spare, and the slack has the expected cause: it assumes the worst case $|y''| = 9.149$ at every step, whereas $y''$ starts at $2$ and only reaches $3$ at $x = 0.1535$ , climbing through $3.424$ at $x = 0.2$ and $3.968$ at $x = 0.25$ before its final $9.149$ ; and it assumes the error is amplified at the maximal rate throughout.

ApplicationOriginal

Every simulation of a continuous system is an instance of this. The equations of a spacecraft trajectory, a reaction network, a power-system transient, a population model or a climate model's physics column are first-order systems after the reduction above, and all of them are advanced by one of the methods here or a descendant.

Why RK4 and not something better. A production solver is usually an embedded Runge-Kutta pair — two methods of adjacent order sharing their stage evaluations, so that the difference of the two answers is an error estimate at almost no extra cost, which drives the step-size controller of step 7. The classical RK4 survives because it is the cheapest scheme whose order is high enough that smooth problems are limited by round-off rather than truncation at a comfortable step size, and because its four stages map cleanly onto a force routine that is called once per stage.

Where order is the wrong figure of merit. Two situations in this chapter defeat the analysis above, and each has its own unit. A system with widely separated time scales forces the step far below what accuracy needs, because an explicit method becomes unstable long before it becomes inaccurate; the diagnosis and the implicit cure are in Stability and stiff equations. A Hamiltonian system integrated over millions of periods suffers a slow drift in its energy that no increase of order removes, because the drift is a structural defect rather than a truncation error; the remedy, a method that is exactly symplectic at every step, is in Symplectic integrators.

Method of lines. Discretising a partial differential equation in space alone (Finite-difference approximations, Spectral methods) leaves a large system of ordinary differential equations in time, to which the methods here apply directly. This is how most time-dependent PDE codes are organised, and it is why the stability restrictions of Stability, the CFL condition and von Neumann analysis and the stiffness of Stability and stiff equations are two views of a single phenomenon: the eigenvalues of the spatial discretisation grow like $(\Delta x)^{-2}$ for a diffusion operator, and the explicit integrator's step must shrink to match.

Theorem

Stability and stiff equations

Adapted

Add to PDF

Feeds the test equation $y' = \lambda y$ through a numerical method to get its stability function $R$ , defines the region of absolute stability as $|R| \le 1$ and A-stability as that region containing the whole left half-plane, proves that no explicit Runge-Kutta method is A-stable while backward Euler and the trapezoidal rule are, computes the real stability intervals $[-2,0]$ for Euler and Heun and $[-2.78529356,\,0]$ for classical Runge-Kutta, and shows on a system with eigenvalues $-1$ and $-1000$ that the explicit step is set by stability rather than accuracy.

Prerequisites: Integrators for ODEs, Eigenvalues and eigenvectors, Linear systems and the matrix exponential, Functions of matrices and the matrix exponential

Definition(The test equation, the stability function, and absolute stability)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6.1, Eqs. (27.61)-(27.63), p. 1021

Riley, Hobson and Bence integrate $y' = -y$ , $y(0) = 1$ by Euler's method and find the recurrence $y_{i+1} = (1-h)y_i$ . They observe that it has the form $y_{i+1} = \lambda y_i$ and that $|\lambda| < 1$ is needed if the computed solution is not to diverge, which here means $0 < h < 2$ ; and they add that this condition does nothing to make the answer accurate. Everything in this unit is the systematic version of that observation.

Take the test equation

y' = \lambda y , \qquad \lambda \in \mathbb{C} ,

whose solution $y(t) = y(0)e^{\lambda t}$ decays exactly when $\operatorname{Re}\lambda < 0$ . Apply a one-step method (Integrators for ODEs) to it with step $h$ and write

z = h\lambda .

For every method of this unit the step then reduces to multiplication by a number depending on $z$ alone,

y_{n+1} = R(z)\,y_n , \qquad\text{so}\qquad y_n = R(z)^{n}y_0 .

The function $R$ is the method's stability function, and $R(z)$ is its amplification factor at step $z$ . The region of absolute stability is

S \;=\; \bigl\{\,z \in \mathbb{C} \;:\; |R(z)| \le 1 \,\bigr\} ,

the set of $z$ for which the computed solution of the test equation does not grow. The method is A-stable when

\mathbb{C}^{-} = \{z : \operatorname{Re} z \le 0\} \subseteq S ,

that is, when every step size reproduces the qualitative decay of every decaying test solution.

Consistency of order $p \ge 1$ forces $R(z) = e^{z} + O(z^{p+1})$ as $z \to 0$ , since the exact one-step map of the test equation is multiplication by $e^{z}$ ; in particular $R(0) = 1$ and $R'(0) = 1$ . So $S$ always touches the origin, and the question is how far into the left half-plane it reaches.

Definition(Stiffness)Original

Let $y' = Ay + g(x)$ with $A$ an $m\times m$ matrix whose eigenvalues $\lambda_1,\dots,\lambda_m$ all have $\operatorname{Re}\lambda_i < 0$ . Its stiffness ratio is

\sigma \;=\; \frac{\max_i |\operatorname{Re}\lambda_i|}{\min_i |\operatorname{Re}\lambda_i|} \;\ge\; 1 .

The problem is called stiff on an interval when $\sigma$ is large and the solution being computed is dominated there by the slow components — that is, when the fast modes have already decayed and are no longer part of the answer.

The operational definition, which is the one that matters, is comparative: a problem is stiff for a given explicit method when the step size permitted by stability is far smaller than the step size permitted by accuracy. Nothing about the exact solution is difficult in such a problem; what is difficult is that an explicit method keeps paying attention to modes that contribute nothing.

For a nonlinear system $y' = F(y)$ the same language is used with $A = \partial F/\partial y$ evaluated along the computed track, a frozen Jacobian. That is a heuristic, not a theorem: the statements below are about constant $A$ , and the remark on hypotheses says what survives.

Definition(Three implicit methods)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.6.3, p. 1025, step (iii)

Backward Euler replaces the slope at the left end of the step by the slope at the right end:

y_{n+1} = y_n + h\,f(x_{n+1}, y_{n+1}) .

The trapezoidal rule averages the two:

y_{n+1} = y_n + \tfrac{h}{2}\bigl(f(x_n,y_n) + f(x_{n+1},y_{n+1})\bigr) .

Riley, Hobson and Bence write exactly this formula as the corrector of their predict-and-correct scheme, on the ground that the average better represents the slope across the step; used on its own, with $y_{n+1}$ solved for rather than guessed, it is a method in its own right. (Applied to a semi-discretised diffusion equation the same formula is the Crank-Nicolson scheme of Stability, the CFL condition and von Neumann analysis.)

The implicit midpoint rule evaluates at the midpoint of the step in both arguments:

y_{n+1} = y_n + h\,f\Bigl(x_n + \tfrac{h}{2},\ \tfrac{1}{2}(y_n + y_{n+1})\Bigr) .

In each case $y_{n+1}$ appears on both sides, so each step requires the solution of an $m$ -dimensional system of equations. For a linear problem $f = Ay$ that is one linear solve; for a nonlinear one it is a Newton iteration (Root finding) on

G(w) = w - y_n - h\,\theta\,f(x_{n+1}, w) - (\text{explicit terms}) = 0 ,

whose Jacobian is $I - h\theta\,\partial f/\partial y$ and whose cost per step is the price of A-stability.

Theorem(Stability functions, stability regions, and A-stability)Original

1 (Diagonalisable linear systems reduce to the test equation). Let $A$ be diagonalisable, $A = P\Lambda P^{-1}$ with $\Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_m)$ , and apply to $y' = Ay$ a method whose test-equation step is multiplication by $R$ . Then

y_n = P\,\operatorname{diag}\bigl(R(h\lambda_1)^{n},\dots,R(h\lambda_m)^{n}\bigr)\,P^{-1}y_0 ,

so $y_n \to 0$ for every $y_0$ if and only if $|R(h\lambda_i)| < 1$ for every $i$ , and $\{y_n\}$ stays bounded for every $y_0$ if and only if $h\lambda_i \in S$ for every $i$ .

2 (The explicit methods). For a method whose stages are substituted one into the next,

method	$R(z)$	$S \cap \mathbb{R}$	$S \cap i\mathbb{R}$
Euler	$1 + z$	$[-2, 0]$	$\{0\}$
Heun, explicit midpoint	$1 + z + \tfrac{z^{2}}{2}$	$[-2, 0]$	$\{0\}$
RK4	$1 + z + \tfrac{z^{2}}{2} + \tfrac{z^{3}}{6} + \tfrac{z^{4}}{24}$	$[z_{*}, 0]$	$[-2\sqrt{2}\,i,\ 2\sqrt{2}\,i]$

where $z_{*}$ is the unique real root of $z^{3} + 4z^{2} + 12z + 24 = 0$ ,

z_{*} = -2.785293563405281623529759\ldots

3 (No explicit Runge-Kutta method is A-stable). If a one-step method is explicit with finitely many stages, each stage being an affine combination of $y_n$ and of earlier stages, then $R$ is a polynomial; if in addition the method is consistent, $R$ is non-constant, so $S$ is bounded and cannot contain $\mathbb{C}^{-}$ .

4 (Backward Euler). $R(z) = (1-z)^{-1}$ . The method has order $1$ with local truncation error $-\tfrac12 h^{2}y''(x) + O(h^{3})$ . It is A-stable; moreover $S$ is the exterior of the disc of radius $1$ about $z = 1$ , so $S$ contains far more than the left half-plane, and $|R(z)| \to 0$ as $\operatorname{Re}z \to -\infty$ .

5 (Trapezoidal rule). $R(z) = (2+z)/(2-z)$ . The method has order $2$ with local truncation error $-\tfrac{1}{12}h^{3}y'''(x) + O(h^{4})$ . It is A-stable, with

S = \mathbb{C}^{-} \quad\text{exactly},

and $|R(z)| = 1$ for every $z$ on the imaginary axis. The implicit midpoint rule has the same stability function and the same order.

Idea of the proofOriginal

Everything follows from one observation: on the test equation, a method built from additions, multiplications by $h\lambda$ and one division produces a rational function $R$ of $z = h\lambda$ , and iterating the method just raises that number to a power. Boundedness of $R(z)^{n}$ over all $n$ is therefore a question about a single complex number, $|R(z)| \le 1$ , and the whole subject is the geometry of the sublevel set $S$ .

That geometry splits the methods in two. An explicit method never divides, so $R$ is a polynomial; polynomials blow up in every direction, so $S$ is a bounded blob and some step size is always too large, however stable the exact solution. An implicit method divides by something like $1 - z$ , which is large exactly where the trouble was, and the same factor that costs a linear solve per step is the factor that makes $S$ unbounded. For backward Euler and the trapezoidal rule the A-stability test reduces to comparing $|1-z|$ with $1$ , or $|2+z|$ with $|2-z|$ , and both comparisons are settled by expanding the squares and cancelling: the imaginary parts drop out and only the sign of $\operatorname{Re}z$ remains.

The real stability intervals in part 2 come from the same squaring. For real $z$ the Euler and Heun amplification factors are real, so one only has to find where they cross $\pm1$ ; for RK4 the truncated exponential stays positive on the negative axis, so the binding condition is $R(z) = 1$ , which factors as $z$ times a cubic. On the imaginary axis the cancellation is more delicate and is the reason RK4 can integrate an oscillation at all while Euler and Heun cannot.

Proof(Parts 1 and 2)Original

Part 1. Substitute $y_n = Pw_n$ . Every method here applies to $y' = Ay$ by replacing each evaluation $f(\cdot) = A(\cdot)$ , so the step map is a matrix polynomial (or rational function) in $hA$ ; call it $R(hA)$ , defined by the same formula as $R(z)$ with $hA$ in place of $z$ and $I$ in place of $1$ . Since $A = P\Lambda P^{-1}$ , every power of $hA$ is $P(h\Lambda)^{k}P^{-1}$ , so $R(hA) = P\,R(h\Lambda)\,P^{-1}$ and $R(h\Lambda) = \operatorname{diag}(R(h\lambda_i))$ because $h\Lambda$ is diagonal. (For an implicit method the inverse $(I - hA)^{-1}$ exists as long as no $h\lambda_i$ equals $1$ , and is likewise diagonalised by $P$ .) Then $y_n = R(hA)^{n}y_0 = P\operatorname{diag}(R(h\lambda_i)^{n})P^{-1}y_0$ . Since $P$ is invertible, $y_n \to 0$ for all $y_0$ exactly when every diagonal entry tends to $0$ , which is $|R(h\lambda_i)| < 1$ ; boundedness for all $y_0$ is $|R(h\lambda_i)| \le 1$ .

Part 2, the stability functions. On $y' = \lambda y$ each stage of a Runge-Kutta method is $\lambda$ times an affine combination of earlier quantities. For Euler, $y_{n+1} = y_n + h\lambda y_n = (1+z)y_n$ . For Heun, $k_1 = \lambda y_n$ and $k_2 = \lambda(y_n + h\lambda y_n)$ , so

y_{n+1} = y_n + \tfrac{h}{2}\lambda\bigl(y_n + y_n + z y_n\bigr) = \Bigl(1 + z + \tfrac{z^{2}}{2}\Bigr)y_n ;

the explicit midpoint method — $y_{n+1} = y_n + h f(x_n + \tfrac h2,\, y_n + \tfrac h2 f_n)$ , not the implicit rule of the definitions above — gives $y_{n+1} = y_n + h\lambda(y_n + \tfrac{z}{2}y_n)$ , the same function. For RK4, writing $k_j = \lambda c_j y_n$ ,

c_1 = 1,\quad c_2 = 1 + \tfrac{z}{2},\quad c_3 = 1 + \tfrac{z}{2}c_2 = 1 + \tfrac{z}{2} + \tfrac{z^{2}}{4},\quad c_4 = 1 + z c_3 = 1 + z + \tfrac{z^{2}}{2} + \tfrac{z^{3}}{4} ,

and

R(z) = 1 + \tfrac{z}{6}\bigl(c_1 + 2c_2 + 2c_3 + c_4\bigr) = 1 + \tfrac{z}{6}\Bigl(6 + 3z + z^{2} + \tfrac{z^{3}}{4}\Bigr) = 1 + z + \tfrac{z^{2}}{2} + \tfrac{z^{3}}{6} + \tfrac{z^{4}}{24} ,

the Taylor polynomial of $e^{z}$ of degree $4$ — as it must be, since the method has order $4$ and the test equation has $y^{(k)} = \lambda^{k}y$ .

Part 2, the real intervals. For Euler, $|1+z| \le 1$ with $z$ real means $-2 \le z \le 0$ . For Heun, $g(z) = 1 + z + z^{2}/2 = \tfrac12\bigl((1+z)^{2}+1\bigr) > 0$ for all real $z$ , so $|g| \le 1$ reduces to $g(z) \le 1$ , that is $z(1 + z/2) \le 0$ , that is $-2 \le z \le 0$ .

For RK4 put $g(z) = 1 + z + z^{2}/2 + z^{3}/6 + z^{4}/24$ . In fact $g > 0$ on all of $\mathbb{R}$ , so $|g| \le 1$ is the single inequality $g(z) \le 1$ . To see the positivity, locate the minimum: $g'(z) = 1 + z + z^{2}/2 + z^{3}/6$ is a cubic whose derivative $g''(z) = 1 + z + z^{2}/2 = \tfrac12\bigl((1+z)^{2}+1\bigr)$ is everywhere positive, so $g'$ is strictly increasing and has one real root, $z = -1.5960716$ . There $g = 0.2703948 > 0$ , and since $g \to +\infty$ at both ends that value is the global minimum. (In the variable $t = -z$ the same polynomial $1 - t + t^{2}/2 - t^{3}/6 + t^{4}/24$ is the degree- $4$ Taylor polynomial of $e^{-t}$ , truncated just before the negative term $-t^{5}/120$ that would eventually pull it down.) So the stability condition is

z + \frac{z^{2}}{2} + \frac{z^{3}}{6} + \frac{z^{4}}{24} \le 0 \qquad\Longleftrightarrow\qquad \frac{z}{24}\bigl(z^{3} + 4z^{2} + 12z + 24\bigr) \le 0 .

The cubic $c(z) = z^{3}+4z^{2}+12z+24$ has $c'(z) = 3z^{2}+8z+12$ with discriminant $64 - 144 < 0$ , so $c' > 0$ everywhere and $c$ is strictly increasing: it has exactly one real root $z_{*}$ , and the two values $c(-3) = -27 + 36 - 36 + 24 = -3$ and $c(-2) = -8 + 16 - 24 + 24 = 8$ bracket it, $-3 < 0 < 8$ , so $z_{*} \in (-3,-2)$ . Solving numerically, $z_{*} = -2.785293563405281623529759$ (to the digits shown; substituting it back leaves the residual $c(z_{*}) = 2.5\times10^{-24}$ ). The product $z\,c(z)$ is $\le 0$ exactly on $[z_{*}, 0]$ , which is the interval claimed.

Part 2, the imaginary axes. For Euler, $|1+iy|^{2} = 1 + y^{2} > 1$ unless $y = 0$ . For Heun,

|R(iy)|^{2} = \Bigl(1 - \tfrac{y^{2}}{2}\Bigr)^{2} + y^{2} = 1 - y^{2} + \tfrac{y^{4}}{4} + y^{2} = 1 + \tfrac{y^{4}}{4} > 1

unless $y = 0$ . For RK4, separating $R(iy)$ into real and imaginary parts,

\operatorname{Re}R(iy) = 1 - \tfrac{y^{2}}{2} + \tfrac{y^{4}}{24}, \qquad \operatorname{Im}R(iy) = y - \tfrac{y^{3}}{6} ,

and squaring and adding,

|R(iy)|^{2} = 1 - \frac{y^{6}}{72} + \frac{y^{8}}{576} = 1 + \frac{y^{6}}{576}\bigl(y^{2} - 8\bigr) .

This is $\le 1$ exactly when $y^{2} \le 8$ , that is $|y| \le 2\sqrt{2} = 2.828427124746$ . $\square$

Proof(Parts 3, 4 and 5)Original

Part 3. Let the method have stages $Y_1,\dots,Y_s$ , each an affine combination of $y_n$ and of $f$ evaluated at earlier stages, and let $y_{n+1}$ be an affine combination of the same quantities. On the test equation every evaluation of $f$ is $\lambda$ times a stage, so by induction each $Y_j$ is $p_j(z)y_n$ with $p_j$ a polynomial, and $y_{n+1} = R(z)y_n$ with $R$ a polynomial. Consistency gives $R(z) = 1 + z + O(z^{2})$ , so $\deg R \ge 1$ .

A non-constant polynomial satisfies $|R(z)| \to \infty$ as $|z| \to \infty$ : writing $R(z) = a_dz^{d} + \cdots$ with $a_d \ne 0$ and $d \ge 1$ , $|R(z)| \ge |a_d||z|^{d} - \sum_{k<d}|a_k||z|^{k} \to \infty$ . Hence $S = \{|R| \le 1\}$ is bounded. The left half-plane is unbounded, so it is not contained in $S$ , and the method is not A-stable. Concretely, along the negative real axis a method of this kind always has a finite stability limit, as the three computed intervals in part 2 show.

Part 4. On the test equation, backward Euler reads $y_{n+1} = y_n + zy_{n+1}$ , so $(1-z)y_{n+1} = y_n$ and $R(z) = (1-z)^{-1}$ , defined for $z \ne 1$ . Now $|R(z)| \le 1$ is $|1 - z| \ge 1$ , that is $|z - 1| \ge 1$ : the complement of the open unit disc centred at $1$ . If $\operatorname{Re}z = x \le 0$ then

|z-1|^{2} = (x-1)^{2} + y^{2} = 1 - 2x + x^{2} + y^{2} \ge 1 ,

because $-2x \ge 0$ and $x^{2}+y^{2}\ge 0$ . So $\mathbb{C}^{-} \subseteq S$ , and the containment is strict since, for instance, $z = 3$ lies in $S$ . As $\operatorname{Re}z \to -\infty$ , $|1-z| \to \infty$ and $|R(z)| \to 0$ .

For the order, substitute the exact solution. Using $y'(x+h) = y' + hy'' + \tfrac12 h^{2}y''' + O(h^{3})$ with everything on the right at $x$ ,

\tau = y(x+h) - y(x) - h y'(x+h) = \Bigl(hy' + \tfrac{h^{2}}{2}y'' + \tfrac{h^{3}}{6}y'''\Bigr) - \Bigl(hy' + h^{2}y'' + \tfrac{h^{3}}{2}y'''\Bigr) + O(h^{4}) = -\tfrac12 h^{2}y'' - \tfrac13 h^{3}y''' + O(h^{4}) ,

so the order is $1$ with the stated leading term.

Part 5. On the test equation the trapezoidal rule reads $y_{n+1} = y_n + \tfrac{z}{2}(y_n + y_{n+1})$ , so $\bigl(1 - \tfrac{z}{2}\bigr)y_{n+1} = \bigl(1 + \tfrac{z}{2}\bigr)y_n$ and

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

With $z = x+iy$ ,

|2+z|^{2} - |2-z|^{2} = \bigl((2+x)^{2}+y^{2}\bigr) - \bigl((2-x)^{2}+y^{2}\bigr) = 8x .

So $|R(z)| \le 1$ if and only if $x \le 0$ : the region of absolute stability is exactly the closed left half-plane, and $|R| = 1$ precisely on the imaginary axis. The method is A-stable, and no more than A-stable.

For the order,

\tau = y(x+h) - y(x) - \tfrac{h}{2}\bigl(y'(x) + y'(x+h)\bigr) .

Expanding both sides to fourth order, $y(x+h)-y(x) = hy' + \tfrac{h^2}{2}y'' + \tfrac{h^3}{6}y''' + \tfrac{h^4}{24}y''''$ , while $\tfrac{h}{2}(y' + y'(x+h)) = hy' + \tfrac{h^{2}}{2}y'' + \tfrac{h^{3}}{4}y''' + \tfrac{h^{4}}{12}y''''$ , all derivatives at $x$ . Subtracting, the first two terms cancel and

\tau = \Bigl(\tfrac16 - \tfrac14\Bigr)h^{3}y''' + \Bigl(\tfrac{1}{24}-\tfrac{1}{12}\Bigr)h^{4}y'''' + O(h^{5}) = -\tfrac{1}{12}h^{3}y''' - \tfrac{1}{24}h^{4}y'''' + O(h^{5}) ,

the order is $2$ , and the constant is $-1/12$ . The implicit midpoint rule gives $\bigl(1 - \tfrac z2\bigr)y_{n+1} = \bigl(1+\tfrac z2\bigr)y_n$ on the test equation as well, hence the same $R$ ; its order is also $2$ . $\square$

Remark(Where the hypotheses are used, and what A-stability does not buy)Original

Diagonalisability, in part 1. If $A$ is defective the diagonal form is replaced by a Jordan form (Jordan normal form) and $R(hA)^{n}$ acquires polynomial-in- $n$ prefactors from the nilpotent part. The strict inequality $|R(h\lambda_i)| < 1$ still forces $y_n \to 0$ , because $n^{k}\rho^{n}\to0$ for $\rho<1$ ; but the boundary case $|R| = 1$ no longer gives boundedness, so the closed region $S$ is the right object only in the diagonalisable case. The same caveat attaches to $P$ : the bound on $\lVert y_n\rVert$ carries the condition number of $P$ (Matrix norms and condition numbers), so a stable-but-highly-non-normal $A$ can produce large transient growth even when every $|R(h\lambda_i)| < 1$ .

Constant coefficients. Parts 1 to 5 are theorems about $y' = \lambda y$ and $y' = Ay$ . Applying them to a nonlinear problem through the frozen Jacobian is an extrapolation that usually works and sometimes does not; it has no proof here and should be treated as a design rule, confirmed by the step-doubling error estimate of Integrators for ODEs, not as a guarantee.

Being inside $S$ is not enough; one needs to be well inside. This is the hypothesis most often mishandled, and it is visible in the example below. At $h\lambda = -2$ exactly, Euler has $|R| = 1$ : the corresponding mode neither grows nor decays, so a fast transient that should have vanished within a thousandth of the interval instead survives to the end, alternating in sign. Riley, Hobson and Bence's $h > 1$ oscillation is the same effect one step earlier: $R(z) < 0$ for $z < -1$ means that every decaying mode is computed as an alternating sequence.

A-stability says nothing about accuracy. It is a statement about the sign of the computed decay, not its rate. Backward Euler is A-stable and first order; at $z = -50$ it returns $R = 0.0196$ where the exact factor is $e^{-50} = 2\times10^{-22}$ , which is a wild overestimate that happens to be harmless because both are negligible.

A-stability does not damp. The trapezoidal rule has $|R| = 1$ on the imaginary axis and $R(z) \to -1$ as $\operatorname{Re}z \to -\infty$ . So a stiff transient is not removed, it is turned into an alternating sequence that decays only as $|(2+z)/(2-z)|^{n}$ , which is close to $1$ for very negative $z$ . The example shows this costing more than half the answer. The property backward Euler has and the trapezoidal rule lacks — $|R(z)| \to 0$ as $\operatorname{Re}z \to -\infty$ , so that infinitely stiff modes are annihilated in one step — is called L-stability, and it is why production stiff solvers are built from backward differentiation formulas and from implicit Runge-Kutta methods with this extra property rather than from the trapezoidal rule.

The cost is not symmetric. Part 3 forbids A-stability only for explicit methods. The price of the implicit alternative is one Newton solve per step, which for a system of dimension $m$ costs $O(m^{3})$ once per Jacobian factorisation and $O(m^{2})$ per iteration (Solving linear systems numerically). Stiffness is therefore worth exploiting only when $\sigma$ is large enough that the saved steps outnumber the added linear algebra; the crossover is problem-dependent and is why every serious library ships both kinds of solver.

Method(Diagnosing and treating stiffness)Original

Get the eigenvalues. For $y' = Ay + g$ compute the spectrum of $A$ ; for a nonlinear system, the spectrum of $\partial F/\partial y$ at a few representative states along the track (Computing eigenvalues). Keep the most negative real part $\mu_{\max}$ and the smallest $\mu_{\min}$ .
Form the two candidate step sizes. The stability step is $h_{\mathrm{stab}} = |z_{\mathrm{edge}}|/\mu_{\max}$ , where $z_{\mathrm{edge}}$ is the left endpoint of $S\cap\mathbb{R}$ for the method in hand ( $2$ for Euler and Heun, $2.785$ for RK4). The accuracy step is the one the error estimator asks for, which is governed by the slowest active mode, so $h_{\mathrm{acc}} \sim \eta^{1/p}/\mu_{\min}$ for a tolerance $\eta$ .
Compare. If $h_{\mathrm{acc}} \le h_{\mathrm{stab}}$ , the problem is not stiff for this method: keep the explicit method. If $h_{\mathrm{acc}} \gg h_{\mathrm{stab}}$ , it is stiff. By step 2 the two scale as $h_{\mathrm{acc}}/h_{\mathrm{stab}} \sim \sigma\,\eta^{1/p}/|z_{\mathrm{edge}}|$ : the wasted effort is proportional to $\sigma$ , but the tolerance factor $\eta^{1/p}/|z_{\mathrm{edge}}|$ in front is small, so the ratio is far below $\sigma$ itself. In the example below $\sigma = 1000$ while the ratio is $10$ .
Recognise it without eigenvalues. A reliable symptom is an adaptive explicit solver whose step size settles at a value far below what the solution's smoothness suggests, with the error estimator reporting that the rejected steps failed spectacularly rather than marginally. A second symptom is a computed solution that oscillates from step to step about a smooth exact one.
Switch to an implicit method, and choose it by what the transient is doing. If fast modes are still active and must be resolved, accuracy decides and A-stability is not the issue. If they must merely be suppressed, pick a method with $|R| \to 0$ at $-\infty$ : backward Euler if first order will do, a backward differentiation formula or an L-stable implicit Runge-Kutta method otherwise. Avoid the trapezoidal rule for a problem whose transient is not resolved.
Solve the stage equations by Newton, not by iteration. Fixed-point iteration on $w = y_n + hf(x_{n+1},w)$ contracts only when $h\,\mathrm{Lip}(f) < 1$ , which is the very restriction being escaped; so form $I - h\,\partial f/\partial y$ , factor it, and reuse the factorisation across steps until the Newton iteration starts needing more than two or three passes.
Exploit structure. If only part of the system is fast — a stiff linear term plus a mild nonlinearity — treat that part implicitly and the rest explicitly, or use the exact exponential of the linear part as in the exponential integrators of Nonhomogeneous systems and Duhamel's principle.

Example(A system with eigenvalues

-1

and

-1000

)Original

Take

A = \begin{pmatrix} -500.5 & 499.5 \\ 499.5 & -500.5\end{pmatrix} , \qquad y' = Ay , \qquad y(0) = \begin{pmatrix}2\\0\end{pmatrix} .

Since $A(1,1)^{\mathsf T} = -(1,1)^{\mathsf T}$ and $A(1,-1)^{\mathsf T} = -1000\,(1,-1)^{\mathsf T}$ , the eigenvalues are $-1$ and $-1000$ and the stiffness ratio is $\sigma = 1000$ . The initial vector is $(1,1)^{\mathsf T} + (1,-1)^{\mathsf T}$ , so

y(t) = e^{-t}\begin{pmatrix}1\\1\end{pmatrix} + e^{-1000t}\begin{pmatrix}1\\-1\end{pmatrix} ,

and by $t = 0.01$ the fast part is already $e^{-10} = 4.5\times10^{-5}$ . We ask for $y_1(1)$ , whose exact value is $e^{-1} + e^{-1000} = 0.367879441171$ to twelve figures. All runs below were done in $30$ -digit arithmetic, so nothing shown is round-off.

What accuracy alone would ask for. Both forward and backward Euler are first order, and on the surviving mode $y' = -y$ a step $h = 0.02$ gives $(1-0.02)^{50} = 0.3641697$ , a relative error of $1.008\%$ ; backward Euler at the same step gives a relative error of $0.992\%$ . So $h_{\mathrm{acc}} \approx 0.02$ for a one-percent answer.

What stability insists on. Part 2 requires $h\lambda \in [-2,0]$ for both eigenvalues, and the binding one is $\lambda = -1000$ : $h \le 0.002$ . The accuracy step is ten times too large, and this is exactly the definition of stiffness.

method	$h$	steps	$h\lambda_{\mathrm{fast}}$	$R$ at the fast mode	computed $y_1(1)$	relative error
forward Euler	$0.0025$	$400$	$-2.5$	$-1.5$	$2.732144\times10^{70}$	$7.4\times10^{70}$
forward Euler	$0.002$	$500$	$-2$	$-1$	$1.367511255$	$2.717$
forward Euler	$0.00125$	$800$	$-1.25$	$-0.25$	$0.3676493967$	$-6.253\times10^{-4}$
backward Euler	$0.05$	$20$	$-50$	$0.019608$	$0.3768894829$	$2.449\times10^{-2}$
backward Euler	$0.02$	$50$	$-20$	$0.047619$	$0.3715278821$	$9.917\times10^{-3}$
trapezoidal	$0.05$	$20$	$-50$	$-0.923077$	$0.5695269199$	$5.481\times10^{-1}$
trapezoidal	$0.02$	$50$	$-20$	$-0.818182$	$0.3679110807$	$8.601\times10^{-5}$
trapezoidal	$0.01$	$100$	$-10$	$-0.666667$	$0.3678763755$	$-8.333\times10^{-6}$

Four things in the table are worth reading carefully.

Outside $S$ , disaster is immediate. At $h = 0.0025$ the fast amplification factor is $-1.5$ , and $1.5^{400} = 2.732\times10^{70}$ . The computed answer is that number; the physics has vanished entirely. The transition is not gradual: $h = 0.0025$ is only $25\%$ above the stability limit.

On the boundary of $S$ , the method is useless but quiet. At $h = 0.002$ the fast factor is exactly $-1$ , so the fast component is reproduced as $(-1)^{500} = 1$ for ever instead of decaying to $e^{-1000}$ . The two contributions to $y_1(1)$ are $0.3675113$ from the slow mode, which is accurate, and $1$ from an undead transient; the answer is wrong by a factor of nearly four with nothing on the screen looking unstable.

Inside $S$ , with room to spare, the explicit method is correct and expensive. At $h = 0.00125$ the fast factor is $-0.25$ and $0.25^{800}$ is below $10^{-481}$ ; the error $6.3\times10^{-4}$ is pure first-order truncation on the slow mode. It cost $800$ steps to get three figures. Backward Euler is no more accurate per step — it is first order on the surviving mode too, and needs about $800$ steps of its own to reach the same $6\times10^{-4}$ — but its gain is that it may take the steps it likes: at $h = 0.02$ its $50$ steps return $0.3715$ , two correct figures, where forward Euler at that step size returns nothing usable at all.

A-stability is not L-stability. The trapezoidal rule at $h = 0.05$ is inside its stability region — the region is the whole left half-plane — and still gets the answer wrong by $55\%$ . Its fast amplification factor is $-0.923077$ , so after $20$ steps the spurious transient still has amplitude $0.923077^{20} = 0.2017$ , which is precisely the discrepancy: $0.3678028 + 0.2017241 = 0.5695269$ . Reducing the step to $0.02$ shrinks $0.818182^{50}$ to $4.4\times10^{-5}$ and the rule's second-order accuracy then shows: at that step its relative error is $8.601\times10^{-5}$ against backward Euler's $9.917\times10^{-3}$ at the same step, smaller by a factor of $115$ . Backward Euler would have to run at $h = 0.04$ to be compared at half the cost, and there its relative error is $1.967\times10^{-2}$ , worse by a factor of $229$ . The trapezoidal rule is the better method here once the transient is dead, and the worse one while it is alive.

Cross-check against the semi-discrete heat equation. Replacing $\partial_x^{2}$ by the centred second difference on a grid of spacing $\Delta x$ turns $u_t = u_{xx}$ into $U' = BU$ , and on the unit interval with Dirichlet ends and $\Delta x = 1/(N+1)$ , $B$ has eigenvalues

\lambda_j = -\frac{4}{(\Delta x)^{2}}\,\sin^{2}\Bigl(\frac{j\pi\,\Delta x}{2}\Bigr), \qquad j = 1,\dots,N ,

which fill out $[-4/(\Delta x)^{2},\,0)$ as the grid is refined — note the $\pi$ in the wavenumber $j\pi$ , without which the smallest eigenvalue comes out wrong. Euler's condition $h\lambda \ge -2$ is binding at $j = N$ , where $|\lambda| \to 4/(\Delta x)^{2}$ , and reads $\Delta t \le (\Delta x)^{2}/2$ , which is exactly the limit Stability, the CFL condition and von Neumann analysis obtains by Fourier analysis of the full scheme. The two calculations are the same calculation. For the stiffness ratio, the slowest mode is $j = 1$ , where $\sin(\pi\Delta x/2) \approx \pi\Delta x/2$ gives $|\lambda_1| \approx \pi^{2}$ — the decay rate of the first heat-equation mode, independent of the grid — so

\sigma = \frac{|\lambda|_{\max}}{|\lambda|_{\min}} \approx \frac{4/(\Delta x)^{2}}{\pi^{2}} = \frac{4}{\pi^{2}(\Delta x)^{2}} = O\bigl((\Delta x)^{-2}\bigr) :

refining the spatial grid makes a diffusion problem more stiff, without limit. (The quantity $4/(\Delta x)^{2}$ on its own is $|\lambda|_{\max}$ , not the ratio.)

ApplicationOriginal

Chemical kinetics is the original stiff problem. A reaction network contains rate constants spanning many orders of magnitude — radical recombination on nanoseconds, bulk conversion on hours — and the quantity of interest is almost always the slow one. Atmospheric chemistry modules inside air-quality and climate models solve tens to hundreds of coupled species equations with $\sigma$ in the range $10^{8}$ to $10^{12}$ ; an explicit method would need a step of microseconds to simulate a day. Every such module is implicit, and most use backward differentiation formulas with a reused Jacobian factorisation, exactly as in step 6 of the method above.

Circuit and control simulation. A circuit with a parasitic capacitance has a time constant thousands of times shorter than the signal it carries; a servo with a fast inner loop and a slow outer one is stiff by design. SPICE-class simulators are built on the trapezoidal rule and on backward differentiation formulas for this reason, and the trapezoidal rule's lack of damping is a known failure mode there, where it appears as point-to-point ringing on a step response.

Method of lines for parabolic equations. As the cross-check above shows, a diffusion operator discretised in space is stiff, with the stiffness growing as $(\Delta x)^{-2}$ . This is the single reason the implicit and semi-implicit schemes of Stability, the CFL condition and von Neumann analysis and The finite element method exist, and why a spectral discretisation (Spectral methods), whose second-derivative eigenvalues reach $-N^{2}/4$ , is treated implicitly in its linear part as a matter of course.

Singular perturbation problems. An equation $\epsilon y'' + ay' + by = 0$ with $\epsilon$ small has one slow root and one of size $a/\epsilon$ : the boundary-layer problems of Singular perturbations and boundary layers are stiff problems seen analytically, and the asymptotic layer analysis is a good guide to where the implicit solver should be allowed to take long steps and where it should not.

What stiffness is not. A chaotic system is hard to integrate for a different reason — the Lipschitz constant in the global error bound of Integrators for ODEs is positive, so errors grow exponentially no matter how stable the method. An implicit method does not help with that; nothing does. Stiffness is a cheap problem with an expensive cure, and chaos is the opposite.

Method · advanced

Symplectic integrators

Adapted

Add to PDF

Builds the Stormer-Verlet (leapfrog) step for a separable Hamiltonian as a sandwich of three shears, proves from $\Psi'^{\mathsf T}J\Psi' = J$ that it is symplectic for every step size, reversible and of order exactly $2$ , shows that for the harmonic oscillator it conserves a shifted energy exactly so that $|H(z_n) - H(z_0)| \le \omega^2h^2\tilde H_0/(4 - \omega^2h^2)$ for all $n$ , and confirms the bound numerically at $0.00125$ over $10^6$ steps against a steady fourth-order Runge-Kutta drift.

Prerequisites: Integrators for ODEs, The canonical equations, Taylor's theorem, Existence and uniqueness of solutions

Definition(Hamiltonian systems, flow maps and symplectic maps)AdaptedAdapted from Lee (2013), Ch. 22, Hamilton's equations (22.11) and Prop. 22.16(a), pp. 575-576; Hassani (2013), §26.4, Def. 26.4.1, p. 801 and the matrix criterion $A^{t}JA = J$, p. 804

Work on the phase space $\mathbb{R}^{2n}$ with coordinates $z = (q, p)$ , $q, p \in \mathbb{R}^{n}$ , and write

J \;=\; \begin{pmatrix} 0 & I_{n} \\ -I_{n} & 0 \end{pmatrix} ,

$I_{n}$ being the identity. For a $C^{2}$ function $H$ on $\mathbb{R}^{2n}$ , the canonical equations of The canonical equations read

\dot q \;=\; \frac{\partial H}{\partial p}, \qquad \dot p \;=\; -\,\frac{\partial H}{\partial q}, \qquad\text{that is}\qquad \dot z \;=\; J\,\nabla H(z) .

Lee writes the same system in Darboux coordinates on a symplectic manifold; we stay on $\mathbb{R}^{2n}$ , which is the case a program actually integrates, and leave the invariant picture to Symplectic manifolds and Hamiltonian vector fields. Write $\varphi_{t}$ for the flow map: $\varphi_{t}(z_{0})$ is the solution at time $t$ through $z_{0}$ , defined for $t$ in the maximal interval given by Existence and uniqueness of solutions. Two facts about $\varphi_{t}$ are the ones an integrator is asked to imitate. First, $H$ is constant along every trajectory, because $\tfrac{d}{dt}H(z) = \nabla H \cdot J \nabla H = 0$ by the antisymmetry of $J$ ; Lee states this as the fact that a function is constant along the integral curves of its own Hamiltonian vector field. Second, $\varphi_{t}$ preserves the symplectic structure: in Lee's language a Hamiltonian vector field is symplectic, so its flow leaves the $2$ -form $\omega = \sum_{i} dq^{i} \wedge dp_{i}$ invariant.

In coordinates that second property is a condition on a Jacobian matrix. Call a $C^{1}$ map $\Psi : U \to \mathbb{R}^{2n}$ , $U \subseteq \mathbb{R}^{2n}$ open, symplectic if

\Psi'(z)^{\mathsf T}\, J\, \Psi'(z) \;=\; J \qquad\text{for every } z \in U ,

where $\Psi'(z)$ is the $2n \times 2n$ Jacobian. This is Hassani's criterion for a matrix to represent a symplectic transformation, imposed pointwise on the derivative. Hassani also records that the symplectic transformations of a fixed symplectic vector space form a group; the only part of that we use is that a product of symplectic matrices is symplectic, which is proved in one line below.

A separable Hamiltonian is one of the form

H(q, p) \;=\; \tfrac12\, p^{\mathsf T} M^{-1} p \;+\; V(q) ,

with $M$ a constant symmetric positive definite mass matrix and $V$ the potential. Its canonical equations are $\dot q = M^{-1}p$ and $\dot p = F(q)$ , where $F = -\nabla V$ is the force. Every mechanical system with a velocity-independent potential and coordinates in which the kinetic energy is a fixed quadratic form is of this type.

A one-step method with step $h$ for such a system is a map $\Psi_{h} : \mathbb{R}^{2n} \to \mathbb{R}^{2n}$ with $z_{n+1} = \Psi_{h}(z_{n})$ , in the sense of Integrators for ODEs. It is a symplectic integrator if $\Psi_{h}$ is a symplectic map for every admissible $V$ and every $h$ .

Definition(The Stormer-Verlet (leapfrog) step)Original

Fix $h > 0$ and set $F = -\nabla V$ . The velocity Verlet or kick-drift-kick step takes $z_{n} = (q_{n}, p_{n})$ to $z_{n+1} = (q_{n+1}, p_{n+1})$ by

p_{n+1/2} \;=\; p_{n} + \tfrac{h}{2}\,F(q_{n}), \qquad q_{n+1} \;=\; q_{n} + h\,M^{-1} p_{n+1/2}, \qquad p_{n+1} \;=\; p_{n+1/2} + \tfrac{h}{2}\,F(q_{n+1}) .

Only one new force evaluation is needed per step: $F(q_{n+1})$ closes the current step and opens the next. In terms of the two elementary maps

K_{s}(q, p) \;=\; \bigl(q,\; p + s\,F(q)\bigr), \qquad D_{s}(q, p) \;=\; \bigl(q + s\,M^{-1}p,\; p\bigr) ,

the step is the sandwich

\Psi_{h} \;=\; K_{h/2} \circ D_{h} \circ K_{h/2} .

$K_{s}$ is the exact time- $s$ flow of the Hamiltonian $V(q)$ alone (which moves $p$ and leaves $q$ fixed), and $D_{s}$ the exact time- $s$ flow of $\tfrac12 p^{\mathsf T}M^{-1}p$ alone. The method is therefore a splitting: it integrates each half of $H$ exactly and interleaves the two.

Eliminating the momenta gives the form in which Stormer wrote the scheme, a three-term recurrence in the positions alone,

q_{n+1} - 2q_{n} + q_{n-1} \;=\; h^{2} M^{-1} F(q_{n}) ,

which holds exactly, not just to leading order, for the sequence $(q_{n})$ produced above. The name leapfrog comes from the staggered variant that stores the half-step momenta $p_{n+1/2}$ and advances positions and momenta alternately; it generates the same positions.

Dropping one of the two half kicks gives the first-order symplectic Euler step $p_{n+1} = p_{n} + hF(q_{n})$ , $q_{n+1} = q_{n} + hM^{-1}p_{n+1}$ , which is also symplectic by the argument below but is only of order $1$ and is not reversible.

Method(What the leapfrog step preserves)Original

Let $M$ be constant, symmetric and positive definite, and let $V \in C^{2}(\mathbb{R}^{n})$ . Let $\Psi_{h}$ be the velocity Verlet step.

Symplecticity. $\Psi_{h}'(z)^{\mathsf T} J \Psi_{h}'(z) = J$ for every $z$ and every $h$ , with no smallness condition on $h$ . Consequently $\det \Psi_{h}'(z) = 1$ , so $\Psi_{h}$ preserves Lebesgue measure on phase space exactly.
Reversibility. $\Psi_{-h} = \Psi_{h}^{-1}$ , and with the momentum flip $\rho(q, p) = (q, -p)$ also $\rho \circ \Psi_{h} \circ \rho \circ \Psi_{h} = \mathrm{id}$ .
Order. If in addition $V \in C^{4}$ , the method has order exactly $2$ . In one degree of freedom with $M = 1$ the local error is

\varphi_{h}(z) - \Psi_{h}(z) \;=\; h^{3}\begin{pmatrix} -\tfrac{1}{6}\,p\,V''(q) \\[2pt] \tfrac{1}{12}\bigl(p^{2}V'''(q) - V'(q)V''(q)\bigr) \end{pmatrix} \;+\; O(h^{4}) .

Hence, by the global-error theorem of Integrators for ODEs, the error after a fixed time is $O(h^{2})$ .

Energy for the harmonic oscillator. Take $n = 1$ , $M = 1$ and $V(q) = \tfrac12\omega^{2}q^{2}$ , and put $\kappa = 1 - \tfrac14\omega^{2}h^{2}$ . If $0 < \omega h < 2$ , then the shifted energy

\tilde H(q, p) \;=\; \tfrac12 p^{2} \;+\; \tfrac12\,\omega^{2}\kappa\, q^{2}

is conserved exactly by the step, $\tilde H(z_{n}) = \tilde H(z_{0})$ for all $n$ , and the true energy satisfies, for all $n \ge 0$ ,

\bigl|H(z_{n}) - H(z_{0})\bigr| \;\le\; \frac{\omega^{2}h^{2}}{4 - \omega^{2}h^{2}}\; \tilde H(z_{0}) .

The bound does not grow with $n$ : the energy error oscillates forever inside a band of width $O(h^{2})$ and never drifts. If $\omega h > 2$ the iteration is unstable and $\lVert z_{n}\rVert \to \infty$ geometrically.

Idea of the proofOriginal

Everything follows from the sandwich $\Psi_{h} = K_{h/2} \circ D_{h} \circ K_{h/2}$ . Each factor is a shear: it changes one half of the coordinates by a function of the other half and leaves that other half alone. The Jacobian of a shear is a triangular block matrix whose off-diagonal block is either the Hessian of $V$ or the matrix $M^{-1}$ — and both are symmetric. Symmetry of that block is exactly the condition for a shear to be symplectic, and symplectic matrices are closed under multiplication, so the chain rule hands us the first claim for free, at any step size.

Reversibility is the palindrome: reading $K_{h/2} D_{h} K_{h/2}$ backwards with $h$ replaced by $-h$ inverts each factor in the right order. A method that is its own inverse under $h \mapsto -h$ cannot carry a leading error term of even degree in $h$ : replacing $h$ by $-h$ multiplies that term by $(-1)^{p+1}$ , while inverting the two maps multiplies it by $-1$ , and the two operations give the same map, so $(-1)^{p+1} = -1$ and the order $p$ is even. A three-line expansion then shows the $h^{3}$ coefficient is not zero, pinning the order at $2$ . Finally, for the harmonic oscillator the step is a $2 \times 2$ matrix with determinant $1$ , and one looks for a quadratic form it preserves. There is one, it differs from $H$ by $O(h^{2})$ , and when $\omega h < 2$ its level sets are ellipses — compact, so the iterates are trapped and the energy cannot wander.

ProofOriginal

Step 1: a shear with a symmetric block is symplectic. Let $S, T$ be real $n \times n$ matrices and

A = \begin{pmatrix} I & 0 \\ S & I\end{pmatrix}, \qquad B = \begin{pmatrix} I & T \\ 0 & I\end{pmatrix} .

Multiplying out,

A^{\mathsf T} J A = \begin{pmatrix} I & S^{\mathsf T} \\ 0 & I\end{pmatrix}\begin{pmatrix} S & I \\ -I & 0\end{pmatrix} = \begin{pmatrix} S - S^{\mathsf T} & I \\ -I & 0 \end{pmatrix}, \qquad B^{\mathsf T} J B = \begin{pmatrix} 0 & I \\ -I & T^{\mathsf T} - T\end{pmatrix} .

Each product therefore equals $J$ precisely when the off-diagonal block is its own transpose: $A$ is symplectic exactly when $S = S^{\mathsf T}$ , and $B$ exactly when $T = T^{\mathsf T}$ . Both matrices are block triangular with identity diagonal blocks, so $\det A = \det B = 1$ .

Step 2: products. If $A^{\mathsf T}JA = J$ and $B^{\mathsf T}JB = J$ then $(AB)^{\mathsf T}J(AB) = B^{\mathsf T}(A^{\mathsf T}JA)B = B^{\mathsf T}JB = J$ . (This is the group property Hassani records for symplectic transformations; only closure is needed here.)

Step 3: the Jacobians of the three factors. The kick $K_{s}$ has Jacobian

K_{s}'(q,p) = \begin{pmatrix} I & 0 \\ s\,F'(q) & I \end{pmatrix} = \begin{pmatrix} I & 0 \\ -s\,\nabla^{2}V(q) & I \end{pmatrix} ,

and the block $-s\nabla^{2}V(q)$ is symmetric because it is a multiple of a Hessian and $V \in C^{2}$ , so mixed partials commute. The drift $D_{s}$ has Jacobian $\begin{pmatrix} I & sM^{-1} \\ 0 & I\end{pmatrix}$ , whose block $sM^{-1}$ is symmetric because $M$ is symmetric and invertible. By Step 1 both Jacobians are symplectic with determinant $1$ .

Step 4: claim 1. By the chain rule, $\Psi_{h}'(z) = K_{h/2}'(w_{2})\,D_{h}'(w_{1})\,K_{h/2}'(z)$ , where $w_{1} = K_{h/2}(z)$ and $w_{2} = D_{h}(w_{1})$ are the intermediate points. Each factor is symplectic by Step 3, so the product is symplectic by Step 2, and $\det \Psi_{h}' = 1 \cdot 1 \cdot 1 = 1$ . Preservation of Lebesgue measure is then the change-of-variables formula: $\lvert \Psi_{h}(E)\rvert = \int_{E} \lvert\det \Psi_{h}'\rvert \,dz = \lvert E \rvert$ for every measurable $E$ . No bound on $h$ entered.

Step 5: claim 2. $K_{s}$ and $D_{s}$ are one-parameter groups: $K_{s}K_{s'} = K_{s+s'}$ because $F(q)$ is evaluated at a $q$ that the kick does not move, and $D_{s}D_{s'} = D_{s+s'}$ for the same reason with the roles exchanged. Hence $K_{-s} = K_{s}^{-1}$ and $D_{-s} = D_{s}^{-1}$ , and

\Psi_{-h} = K_{-h/2} \circ D_{-h} \circ K_{-h/2} = K_{h/2}^{-1} \circ D_{h}^{-1} \circ K_{h/2}^{-1} = \bigl(K_{h/2} \circ D_{h} \circ K_{h/2}\bigr)^{-1} = \Psi_{h}^{-1} .

For the second identity, note $\rho K_{s} \rho = K_{-s}$ and $\rho D_{s}\rho = D_{-s}$ (flipping the momentum reverses the sign of each increment), so $\rho \Psi_{h}\rho = \Psi_{-h} = \Psi_{h}^{-1}$ , which rearranges to $\rho\Psi_{h}\rho\Psi_{h} = \mathrm{id}$ . This is the discrete version of the time-reversal symmetry discussed in Conservative and reversible systems.

Step 6: a reversible method has even order. Suppose the method has order exactly $p \ge 1$ , so that $\varphi_{h}(z) - \Psi_{h}(z) = C(z)h^{p+1} + O(h^{p+2})$ with $C$ continuous and not identically zero. First, inverting both maps negates the leading error: putting $z = \Psi_{h}^{-1}(w)$ in $\Psi_{h}(z) = \varphi_{h}(z) - C(z)h^{p+1} + O(h^{p+2})$ gives $\varphi_{h}(z) = w + C(z)h^{p+1} + O(h^{p+2})$ , and applying $\varphi_{h}^{-1}$ , whose derivative is $I + O(h)$ , yields

\Psi_{h}^{-1}(w) \;=\; \varphi_{h}^{-1}(w) + C(w)h^{p+1} + O(h^{p+2}) ,

where $C(z)$ could be replaced by $C(w)$ because $z = w + O(h)$ and $C$ is continuous. Second, $H$ is autonomous, so the flow is a group and $\varphi_{h}^{-1} = \varphi_{-h}$ ; by Step 5, $\Psi_{h}^{-1} = \Psi_{-h}$ . Therefore the displayed identity says $\varphi_{-h} - \Psi_{-h} = -C h^{p+1} + O(h^{p+2})$ . But the definition of the local error applied with step $-h$ says $\varphi_{-h} - \Psi_{-h} = C\,(-h)^{p+1} + O(h^{p+2})$ . Comparing the two, $(-1)^{p+1} = -1$ , so $p$ is even.

Step 7: the $h^{3}$ coefficient. Take $n = 1$ , $M = 1$ , $V \in C^{4}$ . The scheme gives, with no truncation at all,

q_{n+1} = q_{n} + hp_{n} - \tfrac{h^{2}}{2}V'(q_{n}), \qquad p_{n+1} = p_{n} - \tfrac{h}{2}V'(q_{n}) - \tfrac{h}{2}V'(q_{n+1}) .

Expanding $V'(q_{n+1})$ about $q_{n}$ with Taylor's theorem (Taylor's theorem) and writing $q, p$ for $q_{n}, p_{n}$ ,

p_{n+1} = p - hV' - \tfrac{h^{2}}{2}\,pV'' + \tfrac{h^{3}}{4}\bigl(V'V'' - p^{2}V'''\bigr) + O(h^{4}) ,

all derivatives at $q$ . The exact solution is expanded by repeated differentiation of $\dot q = p$ , $\dot p = -V'(q)$ : $\dddot q = -V''p$ and $\ddot p = -V'' p$ , $\dddot p = V'V'' - p^{2}V'''$ , so

q(h) = q + hp - \tfrac{h^{2}}{2}V' - \tfrac{h^{3}}{6}pV'' + O(h^{4}), \qquad p(h) = p - hV' - \tfrac{h^{2}}{2}pV'' + \tfrac{h^{3}}{6}\bigl(V'V'' - p^{2}V'''\bigr) + O(h^{4}) .

Subtracting gives the local error displayed in claim 3: the $h^{1}$ and $h^{2}$ terms cancel, and the $h^{3}$ terms are $-\tfrac16 pV''$ in the position and $\bigl(\tfrac16 - \tfrac14\bigr)(V'V'' - p^{2}V''') = \tfrac{1}{12}(p^{2}V''' - V'V'')$ in the momentum. This coefficient is not identically zero: already for $V(q) = \tfrac12\omega^{2}q^{2}$ and $p \neq 0$ the first component is $-\tfrac16 p\omega^{2} \neq 0$ . A local error whose leading term is of order $h^{3}$ and does not vanish is, by the definition of order in Integrators for ODEs, a method of order exactly $2$ . This is consistent with Step 6, which had already excluded every odd order.

Step 8: claim 4. With $V = \tfrac12\omega^{2}q^{2}$ the step is linear. Substituting $F(q) = -\omega^{2}q$ into the three substeps and writing $x = \omega^{2}h^{2}$ ,

\begin{pmatrix} q_{n+1} \\ p_{n+1}\end{pmatrix} = R_{h}\begin{pmatrix} q_{n} \\ p_{n}\end{pmatrix}, \qquad R_{h} = \begin{pmatrix} a & h \\ c & a \end{pmatrix}, \qquad a = 1 - \tfrac{x}{2}, \quad c = -\omega^{2}h\kappa, \quad \kappa = 1 - \tfrac{x}{4} .

Two identities drive everything. First,

a^{2} + \omega^{2}\kappa h^{2} = \Bigl(1 - \tfrac{x}{2}\Bigr)^{2} + x\Bigl(1 - \tfrac{x}{4}\Bigr) = 1 - x + \tfrac{x^{2}}{4} + x - \tfrac{x^{2}}{4} = 1 ,

which is also the statement $\det R_{h} = a^{2} - hc = 1$ . Second, $\operatorname{tr}R_{h} = 2a = 2 - x$ .

Now expand $\tilde H(z_{n+1}) = \tfrac12(cq + ap)^{2} + \tfrac12\omega^{2}\kappa(aq + hp)^{2}$ and collect:

the coefficient of $qp$ is $ca + \omega^{2}\kappa a h = a\bigl(-\omega^{2}h\kappa + \omega^{2}\kappa h\bigr) = 0$ ;
the coefficient of $\tfrac12 p^{2}$ is $a^{2} + \omega^{2}\kappa h^{2} = 1$ by the first identity;
the coefficient of $\tfrac12 q^{2}$ is $c^{2} + \omega^{2}\kappa a^{2} = \omega^{2}\kappa\bigl(\omega^{2}\kappa h^{2} + a^{2}\bigr) = \omega^{2}\kappa$ , again by the first identity.

So $\tilde H(z_{n+1}) = \tfrac12 p_{n}^{2} + \tfrac12\omega^{2}\kappa q_{n}^{2} = \tilde H(z_{n})$ : the shifted energy is conserved exactly, for every $h$ .

Assume now $0 < \omega h < 2$ , i.e. $x < 4$ , i.e. $\kappa > 0$ . Then $\tilde H$ is positive definite, so the level set $\{\tilde H = \tilde H_{0}\}$ is an ellipse and the iterates stay on it. The two energies differ by

H(z) - \tilde H(z) \;=\; \tfrac12\omega^{2}(1 - \kappa)q^{2} \;=\; \tfrac18\,\omega^{4}h^{2}q^{2} \;\ge\; 0 ,

and on the ellipse $\tfrac12\omega^{2}\kappa q^{2} \le \tilde H_{0}$ , so $q^{2} \le 2\tilde H_{0}/(\omega^{2}\kappa)$ and

0 \;\le\; H(z_{n}) - \tilde H_{0} \;\le\; \frac{\omega^{4}h^{2}}{8}\cdot\frac{2\tilde H_{0}}{\omega^{2}\kappa} \;=\; \frac{\omega^{2}h^{2}}{4\kappa}\,\tilde H_{0} \;=\; \frac{\omega^{2}h^{2}}{4 - \omega^{2}h^{2}}\,\tilde H_{0} .

Both $H(z_{n})$ and $H(z_{0})$ lie in that interval of length $\omega^{2}h^{2}\tilde H_{0}/(4 - \omega^{2}h^{2})$ , which gives the stated bound, uniformly in $n$ .

If instead $\omega h > 2$ then $\lvert\operatorname{tr}R_{h}\rvert = x - 2 > 2$ while $\det R_{h} = 1$ , so the eigenvalues are real and reciprocal, $\lambda$ and $1/\lambda$ with $\lvert\lambda\rvert > 1$ ; any initial condition with a component along the $\lambda$ -eigenvector grows like $\lvert\lambda\rvert^{n}$ . $\square$

Remark(Where the hypotheses are used)Original

$V \in C^{2}$ is used once, in Step 3, to make $\nabla^{2}V$ symmetric. That symmetry is the whole of symplecticity. If the kick uses a force field that is not a gradient — say $n = 2$ and $F(q) = (q_{2}, 0)$ , so that the Jacobian block is $\tfrac{h}{2}\begin{pmatrix}0 & 1\\ 0 & 0\end{pmatrix}$ — then $S - S^{\mathsf T} \neq 0$ and Step 1 fails at the first factor. Of course such a system is not Hamiltonian to begin with; the point is that the integrator inherits the structure from the problem rather than creating it.
$M$ symmetric is used in the same step for the drift. $M$ constant is what makes $D_{s}$ an exact flow and makes the drift explicit. A configuration-dependent mass matrix, which is what a Lagrangian in curvilinear coordinates produces, destroys both: the kinetic term then depends on $q$ as well as $p$ , the split is no longer into two exactly solvable pieces, and velocity Verlet in the form above does not apply.
$V \in C^{4}$ is needed only for the order statement in claim 3, because the Taylor expansions in Step 7 go to third order. Symplecticity and reversibility need no extra smoothness.
$H$ autonomous is used in Step 6, where $\varphi_{h}^{-1} = \varphi_{-h}$ is the group property of the flow. For a time-dependent Hamiltonian the parity argument has to be run on the extended system.
$\omega h < 2$ is the linear stability limit, and it is sharp: at $\omega h = 2$ the matrix $R_{h}$ has the double eigenvalue $-1$ and is not diagonalisable, so the amplitude grows linearly in $n$ ; beyond it, geometrically. Because leapfrog is explicit it has a bounded stability region, exactly as the explicit methods of Stability and stiff equations do, and it is a poor choice for a stiff problem. In a many-body system the relevant $\omega$ is the fastest vibration present, which is why a molecular-dynamics step is set by the stiffest bond.
Constant $h$ is essential to claim 4, since $\tilde H$ depends on $h$ . Varying the step destroys the bounded energy error even though every individual step is still symplectic; the experiment is in the second example below. A periodic step pattern is harmless, because the composition over one period is again a single fixed symplectic map; what hurts is a step chosen from the current state or from an error estimate.

Remark(The shifted energy in general: backward error analysis)Original

Claim 4 is the exactly solvable case of a general phenomenon. One asks for a modified Hamiltonian $\tilde H_{h}$ whose exact flow reproduces the numerical step. For velocity Verlet in one degree of freedom with $M = 1$ the first correction is

\tilde H_{h}(q,p) \;=\; \tfrac12 p^{2} + V(q) \;+\; h^{2}\Bigl(\tfrac{1}{12}\,p^{2}V''(q) \;-\; \tfrac{1}{24}\,V'(q)^{2}\Bigr) \;+\; O(h^{4}) ,

and we verified by symbolic expansion that the time- $h$ flow of this $\tilde H_{h}$ agrees with $\Psi_{h}$ through order $h^{3}$ , which is one order beyond the method itself. For $V = \tfrac12\omega^{2}q^{2}$ the correction is $h^{2}\bigl(\tfrac{1}{12}\omega^{2}p^{2} - \tfrac{1}{24}\omega^{4}q^{2}\bigr)$ , and expanding the product shows that $\tilde H_{h} = \bigl(1 + \tfrac16\omega^{2}h^{2}\bigr)\tilde H + O(h^{4})$ , where $\tilde H$ is the exactly conserved form of claim 4: the two differ only by the normalisation that makes the flow run at the right speed, since the leapfrog rotation angle per step is $\theta = \arccos\bigl(1 - \tfrac12\omega^{2}h^{2}\bigr) = \omega h\bigl(1 + \tfrac{1}{24}\omega^{2}h^{2} + \cdots\bigr)$ rather than $\omega h$ .

The series $H + h^{2}H_{2} + h^{4}H_{4} + \cdots$ does not converge in general. Truncating it at an $h$ -dependent optimal order gives the statement that for an analytic Hamiltonian and numerical trajectories confined to a compact set, a symplectic method of order $p$ keeps $\lvert H(z_{n}) - H(z_{0})\rvert = O(h^{p})$ over times exponentially long in $1/h$ . That result — backward error analysis, due to Benettin and Giorgilli and developed by Hairer, Lubich and Wanner — is not in any book in this handbook's library, so we state it only for orientation and prove nothing beyond the harmonic case of claim 4. The non-symplectic alternative has no such statement available: the energy error of a general method accumulates, and the examples below show it doing so linearly in $t$ .

Method(Integrating a separable Hamiltonian system with leapfrog)Original

Put the problem in separable form. Identify $q$ , $p$ , the constant symmetric positive definite $M$ , and $V$ with $F = -\nabla V$ . If the force you have is not a gradient, stop: leapfrog buys you nothing, and Integrators for ODEs applies instead.
Fix the step from the fastest frequency. Estimate the largest frequency $\omega_{\max}$ in the problem (for a potential well, $\omega_{\max}^{2}$ is the largest eigenvalue of $M^{-1}\nabla^{2}V$ over the region visited). The hard limit is $\omega_{\max}h < 2$ ; for a sensible phase error take $\omega_{\max}h \lesssim 0.1$ , i.e. at least about sixty steps per period of the fastest mode.
Initialise. Compute $F_{0} = F(q_{0})$ once.
Step, reusing the force: $p \leftarrow p + \tfrac{h}{2}F$ ; $q \leftarrow q + hM^{-1}p$ ; $F \leftarrow F(q)$ ; $p \leftarrow p + \tfrac{h}{2}F$ . The pair $(q, p)$ is synchronised at the end of every step, so output needs no interpolation.
Do not vary $h$ . Keep it fixed for the whole run; if the problem demands adaptivity (a close encounter, say), accept that the long-time energy bound is gone, or switch to a method designed for it.
Monitor $H$ as a diagnostic, not as a correction. Plot $H(z_{n}) - H(z_{0})$ . A bounded oscillation of size $O(h^{2})$ is the expected signature. A steady drift means either a coding error in the force (most often a force that is not the gradient of the potential being monitored) or a step size past the stability limit. Never rescale the velocities to force $H$ back: that is not a symplectic map and it throws away the property you chose the method for.
Check the order by halving $h$ on a short run with a known answer: the error should fall by a factor of $4$ .

Example(The harmonic oscillator: the energy band, and what the alternatives do)Original

Take $n = 1$ , $M = 1$ , $\omega = 1$ , so $H = \tfrac12 p^{2} + \tfrac12 q^{2}$ , with $q_{0} = 1$ , $p_{0} = 0$ and $h = 0.1$ . Then $H(z_{0}) = 0.5$ , $\kappa = 1 - 0.01/4 = 0.9975$ and

\tilde H(z_{0}) = \tfrac12(0)^{2} + \tfrac12(0.9975)(1)^{2} = 0.49875 .

One step by hand. $p_{1/2} = 0 + 0.05\cdot(-1) = -0.05$ ; $q_{1} = 1 + 0.1\cdot(-0.05) = 0.995$ ; $p_{1} = -0.05 + 0.05\cdot(-0.995) = -0.09975$ . Then

H(z_{1}) = \tfrac12(0.09975)^{2} + \tfrac12(0.995)^{2} = 0.004975031\,25 + 0.495012\,5 = 0.499987\,531\,25 ,

\tilde H(z_{1}) = \tfrac12(0.0099500625) + \tfrac12(0.9975)(0.990025) = 0.004975031\,25 + 0.493774\,968\,75 = 0.49875 ,

the shifted energy already back to its initial value to every digit, while the true energy has moved by $-1.2469 \times 10^{-5}$ . A second step gives $q_{2} = 0.98005$ , $p_{2} = -0.1985025$ , $H(z_{2}) = 0.499950622503125$ , $\tilde H(z_{2}) = 0.49875$ .

The predicted band. Claim 4 gives

\bigl|H(z_{n}) - H(z_{0})\bigr| \;\le\; \frac{0.01}{4 - 0.01}\times 0.49875 \;=\; \frac{0.0049875}{3.99} \;=\; 0.00125

for every $n$ . Running $10^{6}$ steps, to $t = 10^{5}$ (about $1.6 \times 10^{4}$ periods), the computed energy error stayed in $[-1.2499999999 \times 10^{-3},\, 9.83 \times 10^{-14}]$ — the lower end is the bound $-h^{2}/8 = -0.00125$ attained to ten digits, reached whenever the orbit passes $q \approx 0$ — and $\tilde H$ drifted by at most $1.14 \times 10^{-13}$ , which is accumulated rounding and not method error.

Order. Integrating to $t = 10$ and comparing $q_{N}$ with $\cos 10 = -0.839071529\ldots$ :

$h$	$\lvert q_{N} - \cos 10\rvert$	ratio
$0.2$	$9.225232\times10^{-3}$
$0.1$	$2.276602\times10^{-3}$	$4.052$
$0.05$	$5.673035\times10^{-4}$	$4.013$
$0.025$	$1.417106\times10^{-4}$	$4.003$
$0.0125$	$3.542044\times10^{-5}$	$4.001$

The ratios converge to $4$ , confirming order $2$ .

What the non-symplectic methods do. For a linear system each classical method is multiplication by a fixed matrix, and claim 1 says the determinant measures how fast phase area is created or destroyed. With $A = \begin{pmatrix}0 & 1\\ -\omega^{2} & 0\end{pmatrix}$ :

Explicit Euler, $I + hA$ , has $\det = 1 + \omega^{2}h^{2} = 1.01$ . For $\omega = 1$ this matrix is exactly $\sqrt{1.01}$ times a rotation, so $H(z_{n}) = (1.01)^{n}H(z_{0})$ identically. At $n = 10^{4}$ that predicts a factor $1.01^{10^{4}} = 1.6358287111890 \times 10^{43}$ ; the computed run gave $H(z_{n})/H(z_{0}) = 1.6358287111889\times10^{43}$ , agreeing to thirteen digits.
Heun's method, $I + hA + \tfrac{h^{2}}{2}A^{2}$ , has $\det = 1 + \tfrac14\omega^{4}h^{4} = 1.000025$ : still expanding, if far more slowly.
Classical fourth-order Runge-Kutta has $\det = 1 - \tfrac{1}{72}\omega^{6}h^{6} + \tfrac{1}{576}\omega^{8}h^{8} = 0.9999999861284722$ , just below $1$ , so it contracts phase area by a factor $0.98622$ over $10^{6}$ steps. The energy follows: the computed error was $-6.935\times10^{-5}$ , $-6.931\times10^{-4}$ and $-6.888\times10^{-3}$ after $10^{4}$ , $10^{5}$ and $10^{6}$ steps — a clean linear drift in $t$ , ten times larger for each factor of ten in the run length, against leapfrog's fixed band of $1.25\times10^{-3}$ at the same step size.

Varying the step. Repeating the $10^{6}$ -step leapfrog run but drawing each step independently and uniformly from $\{0.05, 0.1\}$ , the energy error grew to $1.84\times10^{-4}$ , $2.77\times10^{-3}$ and $5.41\times10^{-2}$ after $10^{4}$ , $10^{5}$ and $10^{6}$ steps (one realisation), a maximum of $5.6\times10^{-2}$ over the run, against the fixed-step $1.25\times10^{-3}$ that never moved. Alternating the two step sizes deterministically, by contrast, kept the error at $9.4\times10^{-4}$ : the two-step composition is itself one fixed symplectic map, and the argument of claim 4 applies to it.

Example(Kepler: a bounded energy error over two thousand orbits)Original

Take the planar two-body problem in relative coordinates with $GM = 1$ and unit reduced mass: $M = I_{2}$ and $V(q) = -1/\lVert q\rVert$ , so $F(q) = -q/\lVert q\rVert^{3}$ . Start at perihelion of the orbit with semi-major axis $a = 1$ and eccentricity $e = 1/2$ :

q_{0} = (0.5,\, 0), \qquad p_{0} = \Bigl(0,\, \sqrt{\tfrac{1+e}{1-e}}\Bigr) = (0,\, 1.7320508075688772) .

Then $H(z_{0}) = \tfrac12(3) - 1/0.5 = -0.5 = -1/(2a)$ , the angular momentum is $L_{0} = q_{1}p_{2} - q_{2}p_{1} = 0.8660254037844386$ , and the period is $2\pi a^{3/2} = 6.283185307179586$ . Take $h = 0.01$ , about $628$ steps per orbit, and run $2000$ orbits, that is $1\,256\,637$ steps.

Leapfrog. The energy error stayed in $[-1.86\times10^{-13},\; 1.359193\times10^{-4}]$ . The maximum $1.359193\times10^{-4}$ was already reached inside the first $250$ orbits and was still exactly that value at $500$ , $1000$ and $2000$ orbits: the error oscillates with the orbital phase (largest near perihelion, where the potential is steepest) and does not accumulate. The angular momentum was conserved to $7.18\times10^{-14}$ , which is rounding. That last fact is exact, not accidental: in a central field $F(q) = -g(\lVert q\rVert)\,q$ is parallel to $q$ , so a kick changes $L$ by $\tfrac{h}{2}(q_{1}F_{2} - q_{2}F_{1}) = 0$ , and a drift changes it by $h(p_{1}p_{2} - p_{2}p_{1}) = 0$ . Each of the three substeps conserves $L$ identically, so leapfrog does.

Fourth-order Runge-Kutta at the same step. The energy error was $-3.790447\times10^{-7}$ , $-7.580583\times10^{-7}$ , $-1.516216\times10^{-6}$ and $-3.032789\times10^{-6}$ at $250$ , $500$ , $1000$ and $2000$ orbits: it doubles every time the orbit count doubles, to four digits. At $2000$ orbits it is still about forty-five times smaller than leapfrog's oscillation, which is what one expects from a fourth-order method at $h = 0.01$ . But it is a drift, not a band: extrapolating the straight line, it reaches $1.359\times10^{-4}$ at about

2000 \times \frac{1.359193\times10^{-4}}{3.032789\times10^{-6}} \;\approx\; 9.0\times10^{4} \quad\text{orbits},

and keeps going, while the leapfrog figure is the same at $10^{4}$ orbits as at $10^{9}$ . Since the semi-major axis is $a = -1/(2H)$ , a secular energy error is a secular error in the orbital period, which is the error that matters over long integrations; leapfrog does not make it.

ApplicationAdaptedAdapted from Lee (2013), Ch. 22, Example 22.18, pp. 576-577

Lee sets out the $n$ -body problem as a Hamiltonian system: point masses $q_{1}, \dots, q_{n}$ in $\mathbb{R}^{3}$ , conservative forces derived from a potential energy $V$ , and the resulting Hamiltonian the sum of kinetic and potential energies. That is exactly the separable form above, with $M$ the diagonal matrix of masses, and it is the reason velocity Verlet is the integrator of molecular dynamics. A simulation of a protein in water runs $10^{8}$ to $10^{9}$ steps of $1$ or $2$ femtoseconds; the physical quantities wanted from it are thermodynamic averages, and these are wrong if the integrator heats or cools the system. A bounded energy error is therefore not an aesthetic preference but the condition under which the sampled ensemble means anything, and a thermostat can only be trusted on top of an integrator that is not itself pumping energy.

The same argument runs in three other places.

Solar system integration. Long-term studies of planetary orbits over $10^{8}$ to $10^{9}$ years use symplectic maps that split the Hamiltonian into Keplerian and interaction parts in the spirit of the splitting above. The question asked — is the system stable, do the eccentricities wander — is precisely a question about secular behaviour, so an integrator with a drifting energy and hence a drifting semi-major axis answers it wrongly, as the Kepler example shows in miniature.
Hybrid Monte Carlo in lattice field theory. A trajectory is generated by a leapfrog integration of a fictitious Hamiltonian system and then accepted with probability $\min\{1, e^{-\Delta H}\}$ . Exact reversibility and exact volume preservation (claims 1 and 2, which hold at any $h$ ) are what make this acceptance test unbiased; the $O(h^{2})$ bound on $\Delta H$ is what keeps the acceptance rate high and lets the step size be chosen on cost grounds alone.
Accelerator physics. Tracking a particle for $10^{8}$ turns of a storage ring to decide whether its motion stays bounded is a question about an invariant torus. A non-symplectic tracker slowly contracts or expands phase-space area, manufacturing damping or blow-up that the magnets do not produce.

The limits are as sharp as the virtues. Symplectic methods are for conservative, autonomous, non-stiff problems integrated at a fixed step. Add dissipation and there is no symplectic structure to preserve; add stiffness and the explicit step limit $\omega_{\max}h < 2$ bites, so the implicit methods of Stability and stiff equations are the right tool; add adaptivity and the long-time energy bound goes, as the first example shows.

Partial differential equations

Method

Finite-difference approximations

Adapted

Add to PDF

Replaces derivatives by difference quotients on a uniform grid, gives each quotient's error as an exact mean-value remainder, defines the local truncation error and the order of a scheme, and applies the recipe to the explicit schemes for the heat and wave equations, whose truncation errors are computed here term by term.

Prerequisites: Taylor's theorem, Classification of second-order equations, The one-dimensional wave equation

Definition(Grid, grid functions, difference quotients)AdaptedAdapted from Strauss (2008), §8.1, pp. 199-201, Eqs. (1)-(6)

Fix a step $h = \Delta x > 0$ and sample a function of one variable at the points $x_j = j\,\Delta x$ , $j \in \mathbb{Z}$ . A grid function is a doubly infinite list of numbers $u_j$ , thought of as approximating $u(x_j)$ . For a problem in space and time we also fix $\Delta t > 0$ , write $t_n = n\,\Delta t$ , and use $u_j^n$ for the approximation to $u(x_j, t_n)$ ; here $n$ is an index, not an exponent.

The four difference quotients used below are the forward, backward and centred first differences and the centred second difference,

D^{+}u_j = \frac{u_{j+1} - u_j}{\Delta x}, \qquad D^{-}u_j = \frac{u_j - u_{j-1}}{\Delta x},

D^{0}u_j = \frac{u_{j+1} - u_{j-1}}{2\,\Delta x}, \qquad \delta^{2}u_j = \frac{u_{j+1} - 2u_j + u_{j-1}}{(\Delta x)^{2}} .

The same four symbols are used in the time index when needed, so that for instance $(u_j^{n+1} - u_j^{n})/\Delta t$ is the forward difference in $t$ .

Method(Each quotient's error, with an exact remainder)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.5, Eqs. (27.57)-(27.59), p. 1019; Strauss (2008), §8.1, p. 200

Let $h > 0$ and let $f$ be real-valued on an interval containing $[x-h, x+h]$ .

If $f \in C^{2}([x, x+h])$ there is a $\xi \in (x, x+h)$ with

\frac{f(x+h) - f(x)}{h} = f'(x) + \frac{h}{2}\,f''(\xi) .

If $f \in C^{2}([x-h,x])$ there is a $\xi \in (x-h,x)$ with $\bigl(f(x) - f(x-h)\bigr)/h = f'(x) - \tfrac{h}{2}f''(\xi)$ .

If $f \in C^{3}([x-h, x+h])$ there is a $\xi \in (x-h, x+h)$ with

\frac{f(x+h) - f(x-h)}{2h} = f'(x) + \frac{h^{2}}{6}\,f'''(\xi) .

If $f \in C^{4}([x-h, x+h])$ there is a $\xi \in (x-h, x+h)$ with

\frac{f(x+h) - 2f(x) + f(x-h)}{h^{2}} = f''(x) + \frac{h^{2}}{12}\,f^{(4)}(\xi) .

In particular the one-sided quotients are first order in $h$ and the two centred ones are second order, with the leading error constants $h/2$ , $h^2/6$ and $h^2/12$ reported by Riley, Hobson and Bence.

Idea of the proofOriginal

Everything comes from the two Taylor expansions of $f$ about $x$ , one stepping forward by $h$ and one stepping back. Subtracting them cancels the even-order terms, so the odd-order derivative $f'$ survives at leading order and the first surviving error is $f'''$ ; adding them cancels the odd-order terms, so $f''$ survives and the first error is $f^{(4)}$ . That is the whole reason the centred quotients gain an order over the one-sided ones. The only extra care needed is that each expansion leaves its own remainder point, and the two must be merged into one: continuity of the top derivative lets the intermediate value theorem replace the average of two of its values by a single value.

ProofAdaptedAdapted from Strauss (2008), §8.1, p. 200; Riley, Hobson & Bence (2006), §27.5, Eq. (27.56), p. 1019

Step 1 (the one-sided quotients). Taylor's theorem with Lagrange remainder, applied at order one to $f \in C^{2}$ , gives a point $\xi \in (x, x+h)$ with $f(x+h) = f(x) + h f'(x) + \tfrac{1}{2}h^{2} f''(\xi)$ . Subtract $f(x)$ and divide by $h$ : this is claim 1. Applying the same expansion with $-h$ in place of $h$ gives the backward statement, the sign of the remainder flipping because dividing $\tfrac12(-h)^2 f''(\xi)$ by $-h$ produces $-\tfrac{h}{2}f''(\xi)$ .

Step 2 (two expansions to third order). Let $f \in C^{3}$ . Taylor's theorem at order two gives points $\xi_{+} \in (x, x+h)$ and $\xi_{-} \in (x-h, x)$ with

f(x \pm h) = f(x) \pm h f'(x) + \tfrac{1}{2}h^{2} f''(x) \pm \tfrac{1}{6}h^{3} f'''(\xi_{\pm}) .

Step 3 (subtract). Subtracting the two displays in Step 2 cancels $f(x)$ and the $f''$ term, leaving

f(x+h) - f(x-h) = 2h f'(x) + \tfrac{1}{6}h^{3}\bigl(f'''(\xi_{+}) + f'''(\xi_{-})\bigr) .

Since $f'''$ is continuous on $[\xi_{-}, \xi_{+}]$ , the number $\tfrac12\bigl(f'''(\xi_{+}) + f'''(\xi_{-})\bigr)$ lies between $f'''(\xi_-)$ and $f'''(\xi_+)$ , so by the intermediate value theorem it equals $f'''(\xi)$ for some $\xi$ in that interval, hence in $(x-h,x+h)$ . The bracket is therefore $2f'''(\xi)$ , and dividing by $2h$ gives claim 2.

Step 4 (add, one order higher). Let $f \in C^{4}$ . Taylor's theorem at order three gives $\eta_{\pm}$ with

f(x \pm h) = f(x) \pm h f'(x) + \tfrac{1}{2}h^{2}f''(x) \pm \tfrac{1}{6}h^{3}f'''(x) + \tfrac{1}{24}h^{4}f^{(4)}(\eta_{\pm}) .

Adding these cancels the two odd-order terms:

f(x+h) - 2f(x) + f(x-h) = h^{2}f''(x) + \tfrac{1}{24}h^{4}\bigl(f^{(4)}(\eta_{+}) + f^{(4)}(\eta_{-})\bigr) .

The continuity of $f^{(4)}$ and the same intermediate-value step turn the bracket into $2f^{(4)}(\xi)$ , and dividing by $h^{2}$ gives claim 3, the constant being $2/24 = 1/12$ . $\square$

Remark(Where the hypotheses are used)Original

Smoothness buys the order, and nothing else does. Each claim uses exactly the differentiability that appears in its remainder term, and the stated order fails as soon as that is dropped.

One degree short of $C^2$ . Let $f(x) = |x|^{3/2}$ , which is $C^{1}$ on $\mathbb{R}$ with $f'(x) = \tfrac32 |x|^{1/2}\operatorname{sgn} x$ and $f'(0) = 0$ , but has no second derivative at $0$ . The forward quotient at $x = 0$ is $h^{3/2}/h = h^{1/2}$ , so its error is $\sqrt{h}$ : it still tends to $0$ , but not at rate $h$ . Claim 1 is unavailable because there is no $f''$ to evaluate at $\xi$ .
One degree short of $C^4$ . Let $f(x) = |x|^{7/2}$ , which is $C^{3}$ on $\mathbb{R}$ : here $f''(x) = \tfrac{35}{4}|x|^{3/2}$ and $f'''(x) = \tfrac{105}{8}|x|^{1/2}\operatorname{sgn} x$ are continuous with $f''(0) = f'''(0) = 0$ , while $f^{(4)}(x) = \tfrac{105}{16}|x|^{-1/2}$ blows up at the origin. The centred second difference at $x = 0$ is $(h^{7/2} + h^{7/2})/h^{2} = 2h^{3/2}$ , an error of order $h^{3/2}$ , not $h^{2}$ . Dropping one more degree makes it worse still: for $f(x) = |x|^{3}$ , which is $C^{2}$ but not $C^{3}$ , the same quotient gives $2h$ . So the second-order rate of claim 3 is genuinely bought by the fourth derivative.
Continuity, not mere existence, of the top derivative. Steps 3 and 4 merge two remainder points into one by the intermediate value theorem, which needs $f'''$ (respectively $f^{(4)}$ ) to be continuous on the closed interval. Without continuity one can still bound the error by the supremum of the top derivative, but the clean mean-value identity is lost.
The identities are exact. Because the remainders are Lagrange remainders and not $O(\cdot)$ symbols, the three displays are equalities, not asymptotic statements. That is what makes the error constants testable at a single finite $h$ , as in the first example below.

Definition(Local truncation error, consistency, order)AdaptedAdapted from Strauss (2008), §8.1, p. 201

Write a linear evolution problem as $\mathcal{L}u = 0$ (for example $\mathcal{L} = \partial_t - \partial_x^2$ ), and let $\mathcal{L}_\Delta$ be the difference operator obtained by replacing every derivative in $\mathcal{L}$ by a difference quotient on the grid. The local truncation error of the scheme at a grid point is

\tau_j^{n} \;=\; \bigl(\mathcal{L}_\Delta u\bigr)(x_j, t_n),

the residual left when the exact solution $u$ of $\mathcal{L}u = 0$ is substituted into the difference equation. The scheme is consistent if $\tau_j^n \to 0$ as $\Delta x, \Delta t \to 0$ , and has order $(p, q)$ if $\tau_j^n = O\bigl((\Delta t)^{p} + (\Delta x)^{q}\bigr)$ uniformly on the region of interest.

Strauss separates three distinct errors, and so do we. The local truncation error above is the error made in the equation at one point. The global truncation error is the error $u(x_j,t_n) - u_j^n$ actually present in the computed solution, built up from all the local ones; he observes that following that accumulation in detail is usually out of reach, and the two can behave quite differently. Roundoff error is a separate effect, caused by the finite number of digits the machine keeps at each step, and is treated in Floating-point arithmetic and error.

Consistency alone says nothing about whether the computed solution is anywhere near the true one; what has to be added is stability, which is the subject of Stability, the CFL condition and von Neumann analysis.

Definition(The explicit scheme for the heat equation)AdaptedAdapted from Strauss (2008), §8.1, Eq. (7), p. 201; Strauss (2008), §8.2, Eqs. (1)-(2), p. 203; Strauss (2008), §8.2, Eq. (14), p. 208

For $u_t = u_{xx}$ , take the forward difference in $t$ and the centred second difference in $x$ :

\frac{u_j^{n+1} - u_j^{n}}{\Delta t} \;=\; \frac{u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}}{(\Delta x)^{2}} . \tag{H}

With the mesh ratio

s \;=\; \frac{\Delta t}{(\Delta x)^{2}}

this solves for the new row in terms of the old one,

u_j^{n+1} \;=\; s\,\bigl(u_{j+1}^{n} + u_{j-1}^{n}\bigr) + (1 - 2s)\,u_j^{n} ,

so the computation marches forward in time with no linear system to solve. A scheme with this property is called explicit. Scheme (H) is often labelled FTCS, for forward in time and centred in space.

On a bounded interval with Dirichlet data the end values $u_0^n$ and $u_J^n$ are prescribed. For a Neumann condition $u_x = g$ at an end, Strauss recommends a ghost point just outside the interval and a centred difference across the boundary, so that the boundary is treated to the same second-order accuracy as the interior; a one-sided quotient there would contaminate the computation with an $O(\Delta x)$ error by claim 1 above.

Definition(The explicit scheme for the wave equation)AdaptedAdapted from Strauss (2008), §8.3, Eqs. (1)-(2), p. 211 and Eqs. (5)-(6), p. 212; Strauss (2008), §8.4, Eqs. (1)-(2), pp. 218-219

For $u_{tt} = c^{2}u_{xx}$ with $c > 0$ , use the centred second difference in both variables:

\frac{u_j^{n+1} - 2u_j^{n} + u_j^{n-1}}{(\Delta t)^{2}} \;=\; c^{2}\,\frac{u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}}{(\Delta x)^{2}} . \tag{W}

Writing $\lambda = c\,\Delta t/\Delta x$ (so that Strauss's parameter is $s = \lambda^{2}$ ) and solving for the new row,

u_j^{n+1} \;=\; \lambda^{2}\bigl(u_{j+1}^{n} + u_{j-1}^{n}\bigr) + 2\bigl(1 - \lambda^{2}\bigr)u_j^{n} - u_j^{n-1} .

This is again explicit, but it reaches back two rows, so two starting rows are needed. Given $u(x,0) = \varphi(x)$ and $u_t(x,0) = \psi(x)$ , Strauss builds them by imposing the centred difference $(u_j^{1} - u_j^{-1})/(2\Delta t) = \psi_j$ and eliminating the fictitious row $u_j^{-1}$ with (W) at $n = 0$ , which gives

u_j^{0} = \varphi_j , \qquad u_j^{1} = \tfrac{1}{2}\lambda^{2}\bigl(\varphi_{j+1} + \varphi_{j-1}\bigr) + \bigl(1 - \lambda^{2}\bigr)\varphi_j + \psi_j\,\Delta t .

A one-sided start $u_j^1 = \varphi_j + \psi_j \Delta t$ would be only first order and would spoil the second-order accuracy of (W).

The corresponding elliptic discretisation, the five-point formula that makes $u_{j,k}$ the average $\tfrac14(u_{j+1,k} + u_{j-1,k} + u_{j,k+1} + u_{j,k-1})$ of its four neighbours, is treated in Relaxation methods for elliptic problems. Which of the three prototypes one is discretising is decided by the classification of Classification of second-order equations, and it controls whether a marching method exists at all: for the elliptic case there is no time direction to march in, and a large linear system must be solved instead.

Method(Writing down a scheme and finding its order)Original

Classify the equation and decide what kind of problem it is. For a hyperbolic or parabolic equation in time, the grid marches; for an elliptic equation there is no marching direction.
Choose the grid $x_j = j\,\Delta x$ , $t_n = n\,\Delta t$ , and name the mesh ratio that will appear ( $s = \Delta t/(\Delta x)^2$ for diffusion, $\lambda = c\,\Delta t/\Delta x$ for waves).
Replace each derivative by a quotient from the definition above. Use centred quotients wherever the equation's symmetry allows; a one-sided quotient in time is what makes a scheme explicit.
Substitute the exact solution into the resulting difference equation and expand every term by Taylor's theorem about the same grid point $(x_j, t_n)$ , keeping two orders beyond the one you expect to survive.
Use the differential equation itself to convert time derivatives into space derivatives (for $u_t = u_{xx}$ : $u_{tt} = u_{xxxx}$ , $u_{ttt} = u_{xxxxxx}$ ). Without this step the leading terms will not combine.
Collect powers of $\Delta x$ after substituting $\Delta t$ in terms of $\Delta x$ through the mesh ratio. The lowest surviving power is the order; a mesh ratio that annihilates the leading coefficient raises the order.
Match the boundary and starting formulas to the interior order, using ghost points for derivative conditions and a centred start for a second-order-in-time scheme.
Check the order numerically: evaluate the residual on a known exact solution at several step sizes and confirm that halving $\Delta x$ divides it by $2^{q}$ .

Example(The error constants are real numbers, and they are right)Original

Take $f(x) = e^{x}$ at $x = 1$ with $h = 0.1$ , so that every derivative at the point equals $e = 2.718281828459\ldots$ All values below were computed in $25$ -digit arithmetic and are quoted to the digits shown.

Centred second difference.

\frac{e^{1.1} - 2e^{1} + e^{0.9}}{0.01} = 2.720547818529 .

Its error is $2.720547818529 - 2.718281828459 = 2.265990\times10^{-3}$ . Claim 3 predicts a leading error $\tfrac{h^2}{12}f^{(4)}(\xi)$ , which for $\xi \approx 1$ is $\tfrac{0.01}{12}\,e = 2.265235\times10^{-3}$ . The next term of the expansion is $\tfrac{h^{4}}{360}e = 7.5508\times10^{-7}$ , and

2.265235\times10^{-3} + 7.5508\times10^{-7} = 2.265990\times10^{-3},

which reproduces the measured error in all seven figures shown.

Centred first difference. $(e^{1.1} - e^{0.9})/0.2 = 2.722814563947$ , an error of $4.532735\times10^{-3}$ . Claim 2 predicts $\tfrac{h^{2}}{6}e = 4.530470\times10^{-3}$ , and adding the next term $\tfrac{h^{4}}{120}e = 2.2652\times10^{-6}$ gives $4.532735\times10^{-3}$ : again all seven figures.

Forward difference. $(e^{1.1} - e^{1})/0.1 = 2.858841954874$ , an error of $1.405601\times10^{-1}$ against the predicted $\tfrac{h}{2}e = 1.359141\times10^{-1}$ . The agreement is only to two figures here, and that is as it should be: the remainder $\tfrac{h}{2}f''(\xi)$ is evaluated not at $1$ but at some $\xi \in (1, 1.1)$ , where $f''$ is larger by up to the factor $e^{h} = 1.1052$ . The measured ratio of error to prediction is $1.405601/1.359141 = 1.0342$ , so $f''(\xi) = e^{\xi}$ with $\xi = 1.0336$ , comfortably inside the interval the theorem allows. At this step size the one-sided quotient is worse than the centred one by a factor of $31$ .

Example(Truncation error of the heat scheme, and a mesh ratio that gains two orders)Original

Let $u$ solve $u_t = u_{xx}$ and let $\tau$ be the local truncation error of (H),

\tau \;=\; \frac{u(x, t+\Delta t) - u(x,t)}{\Delta t} \;-\; \frac{u(x+\Delta x, t) - 2u(x,t) + u(x-\Delta x, t)}{(\Delta x)^{2}} .

Step 1 (expand). By Taylor's theorem in $t$ , and by the expansion behind claim 3 carried one term further in $x$ ,

\frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = u_t + \frac{\Delta t}{2}u_{tt} + \frac{(\Delta t)^{2}}{6}u_{ttt} + O\bigl((\Delta t)^{3}\bigr),

\frac{u(x+\Delta x,t) - 2u(x,t) + u(x-\Delta x,t)}{(\Delta x)^{2}} = u_{xx} + \frac{(\Delta x)^{2}}{12}u_{xxxx} + \frac{(\Delta x)^{4}}{360}u_{xxxxxx} + O\bigl((\Delta x)^{6}\bigr),

all derivatives being taken at $(x,t)$ . Written with that explicit $O\bigl((\Delta x)^{6}\bigr)$ remainder, the second line needs $u$ to be $C^{8}$ in $x$ (the expansion is carried to order seven with a Lagrange remainder); $C^{6}$ would give the $(\Delta x)^{4}$ term only with an $o\bigl((\Delta x)^{4}\bigr)$ error. Solutions of the heat equation are $C^{\infty}$ in $x$ , so neither requirement bites here.

Step 2 (use the equation). The terms $u_t$ and $u_{xx}$ cancel. Differentiating $u_t = u_{xx}$ in $t$ gives $u_{tt} = u_{xxxx}$ , and once more $u_{ttt} = u_{xxxxxx}$ .

Step 3 (collect in $\Delta x$ ). Put $\Delta t = s(\Delta x)^{2}$ :

\tau \;=\; (\Delta x)^{2}\Bigl(\frac{s}{2} - \frac{1}{12}\Bigr)u_{xxxx} \;+\;(\Delta x)^{4}\Bigl(\frac{s^{2}}{6} - \frac{1}{360}\Bigr)u_{xxxxxx} \;+\;O\bigl((\Delta x)^{6}\bigr).

So (H) is consistent, of order one in $\Delta t$ and two in $\Delta x$ , hence of order two in $\Delta x$ at fixed $s$ — except at $s = 1/6$ , where the bracket $s/2 - 1/12$ vanishes and the scheme becomes fourth order in $\Delta x$ .

Step 4 (check it). Take the exact solution $u(x,t) = e^{-t}\sin x$ , for which $u_{xxxx} = u$ and $u_{xxxxxx} = -u$ , and evaluate $\tau$ at $x = 1$ , $t = 0$ , where $u = \sin 1 = 0.8414710$ .

$\Delta x$	$\tau$ at $s = 1/4$	ratio	$\tau$ at $s = 1/6$	ratio
$0.1$	$3.499706\times10^{-4}$		$-1.557074\times10^{-7}$
$0.05$	$8.761306\times10^{-5}$	$3.9945$	$-9.737364\times10^{-9}$	$15.991$
$0.025$	$2.191080\times10^{-5}$	$3.9986$	$-6.086735\times10^{-10}$	$15.998$
$0.0125$	$5.478170\times10^{-6}$	$3.9997$	$-3.804348\times10^{-11}$	$15.999$

The ratios tend to $4$ and to $16$ , confirming orders $2$ and $4$ . The predicted leading values agree as well. At $s = 1/4$ and $\Delta x = 0.1$ the formula gives $(0.01)\bigl(\tfrac18 - \tfrac1{12}\bigr)(0.8414710) = 3.506129\times10^{-4}$ against the measured $3.499706\times10^{-4}$ . At $s = 1/6$ the surviving coefficient is

\Bigl(\frac{s^{2}}{6} - \frac{1}{360}\Bigr)u_{xxxxxx} = \Bigl(\frac{1}{216} - \frac{1}{360}\Bigr)(-u) = -\frac{u}{540},

so $\tau \approx -(\Delta x)^{4}(0.8414710)/540 = -1.5583\times10^{-7}$ at $\Delta x = 0.1$ , against the measured $-1.5571\times10^{-7}$ .

A warning the next unit makes precise: $s = 1/6$ is a consistency bonus only. It happens to lie inside the stability range of (H), but nothing in this calculation says so.

Example(Truncation error of the wave scheme, and the step that is exact)Original

The same procedure applied to (W) is shorter. For a solution of $u_{tt} = c^{2}u_{xx}$ ,

\tau = u_{tt} + \frac{(\Delta t)^{2}}{12}u_{tttt} - c^{2}u_{xx} - \frac{c^{2}(\Delta x)^{2}}{12}u_{xxxx} + O\bigl((\Delta x)^{4}\bigr),

the first and third terms cancelling. Differentiating the equation twice in $t$ gives $u_{tttt} = c^{4}u_{xxxx}$ , and with $\Delta t = \lambda\,\Delta x/c$ ,

\tau \;=\; \frac{c^{2}(\Delta x)^{2}}{12}\bigl(\lambda^{2} - 1\bigr)u_{xxxx} \;+\; O\bigl((\Delta x)^{4}\bigr).

Taking $c = 2$ and the exact solution $u = \sin x\,\cos 2t$ , at $x = 1$ , $t = 0.3$ and $\lambda = 1/2$ the computed residuals are $-1.735517\times10^{-3}$ , $-4.340148\times10^{-4}$ , $-1.085122\times10^{-4}$ and $-2.712857\times10^{-5}$ for $\Delta x = 0.1, 0.05, 0.025, 0.0125$ , with ratios $3.9988$ , $3.9997$ and $3.9999$ . At this point $u_{xxxx} = \sin 1\,\cos 0.6 = 0.6944960$ (not $\sin 1$ : the example is evaluated at $t = 0.3$ , not at $t = 0$ ), so the predicted leading value at $\Delta x = 0.1$ is $\tfrac{4(0.01)}{12}\bigl(-\tfrac34\bigr)(0.6944960) = -1.736240\times10^{-3}$ , against the measured $-1.735517\times10^{-3}$ .

At $\lambda = 1$ the coefficient vanishes, and the computed residual drops to $5\times10^{-39}$ in $40$ -digit arithmetic, that is, to zero. This is not an accident of the leading term. With $c\,\Delta t = \Delta x$ the scheme reads $u_j^{n+1} = u_{j+1}^{n} + u_{j-1}^{n} - u_j^{n-1}$ ; substituting a travelling wave $u = f(x - ct)$ and writing $z = x_j - ct_n$ , the three terms on the right are $f(z + \Delta x)$ , $f(z - \Delta x)$ and $f(z + \Delta x)$ , whose combination is $f(z - \Delta x)$ , and that is exactly $u_j^{n+1}$ . The same holds for $g(x+ct)$ , so by d'Alembert's formula (The one-dimensional wave equation) the scheme reproduces every solution of the one-dimensional wave equation exactly when $\lambda = 1$ . Strauss observes numerically that $\lambda = 1$ gives an excellent approximation while $\lambda^2 = 2$ gives a useless one; the first half of that is explained here, the second half in Stability, the CFL condition and von Neumann analysis.

ApplicationOriginal

Finite differences are how most of the partial differential equations of geophysics and engineering are actually solved, and the error constants above are what set the cost of a run.

In an atmospheric or ocean model the horizontal grid spacing is chosen first, from the features that must be resolved, and the vertical diffusion of heat and momentum is then discretised by a scheme of the form (H). Because the order in $\Delta x$ is two, dividing the grid spacing by two to cut the discretisation error by four multiplies the number of columns by four in two horizontal dimensions; the companion stability constraint on $\Delta t$ , derived in Stability, the CFL condition and von Neumann analysis, multiplies the number of time steps as well, so the total cost grows by a factor of at least sixteen. That single arithmetic fact is why operational forecasting centres buy resolution in small increments, and why higher-order and spectral discretisations (Spectral methods) remain attractive even though each step of them costs more.

In seismic modelling the scheme (W) is the standard tool, and the exactness at $\lambda = 1$ in one dimension is the reason practitioners speak of a magic time step. In more than one space dimension no such step exists, but $\lambda$ is still pushed as close to its stability limit as the geometry allows, both for speed and because, as the formula for $\tau$ shows, the truncation error shrinks as $\lambda \to 1$ .

Theorem

Stability, the CFL condition and von Neumann analysis

Adapted

Add to PDF

Feeds a single discrete Fourier mode through a constant-coefficient difference scheme to read off its amplification factor, proves that this gives the limits $\Delta t \le (\Delta x)^2/2$ for the explicit heat scheme and $c\,\Delta t \le \Delta x$ for the explicit wave scheme, derives the Courant-Friedrichs-Lewy condition from the domains of dependence, and records the Lax equivalence theorem that converts stability plus consistency into convergence.

Prerequisites: Finite-difference approximations, The discrete Fourier transform and the FFT

Definition(Discrete modes and the amplification factor)AdaptedAdapted from Strauss (2008), §8.2, Eqs. (4)-(6) and (13), pp. 205-207; Strauss (2008), §8.3, p. 213

Work on the whole line, or on a periodic grid, with a constant-coefficient one-step scheme: there are finitely many real or complex numbers $a_m$ , independent of $j$ and $n$ , with

u_j^{n+1} \;=\; \sum_{m} a_m\, u_{j+m}^{\,n} .

Substitute a single discrete mode

u_j^{n} \;=\; \xi^{\,n}\,e^{i k j \Delta x},

a geometric sequence in time times a grid plane wave in space. The spatial factor comes out of the sum, and what is left is a number,

\boxed{\;\xi(\theta) \;=\; \sum_m a_m\,e^{i m \theta}\;} \qquad \theta = k\,\Delta x ,

the amplification factor of the scheme. The mode is then reproduced exactly at the next time level, multiplied by $\xi$ ; after $n$ steps its size is $|\xi|^{n}$ times its initial size.

Only $\theta = k\,\Delta x$ enters, and $e^{ikj\Delta x}$ is unchanged when $\theta$ moves by $2\pi$ : two wavenumbers whose $\theta$ differ by a multiple of $2\pi$ give literally the same grid function. This is aliasing, exactly as in The discrete Fourier transform and the FFT, and it means that $\theta$ ranges over $[-\pi, \pi]$ and no further. The shortest wave the grid can carry is $\theta = \pi$ , the sawtooth $(-1)^{j}$ with two grid points per wavelength.

The scheme satisfies the von Neumann condition if there is a constant $C$ , independent of $k$ and of the step sizes, with

|\xi(\theta)| \;\le\; 1 + C\,\Delta t \qquad\text{for all } \theta \in [-\pi,\pi]

and all sufficiently small $\Delta t$ . Strauss states the condition in this form, citing Richtmyer and Morton, and notes that the slack $C\,\Delta t$ is what is needed for problems whose exact solutions grow in time; when the exact solution does not grow, the sharper requirement $|\xi| \le 1$ is the one to use, and it is the one tested below. The slack is harmless because $(1 + C\Delta t)^{n} \le e^{CT}$ whenever $n\,\Delta t \le T$ , so a mode can grow by at most a fixed factor over a fixed time interval, whereas a violation of the condition by a fixed amount compounds without limit as the grid is refined.

Theorem(Stability limits of the two explicit model schemes)AdaptedAdapted from Strauss (2008), §8.2, Eqs. (9)-(10), pp. 206-207; Strauss (2008), §8.3, Eqs. (7)-(9), pp. 213-214

Let $\theta = k\,\Delta x \in [-\pi,\pi]$ .

(a) The explicit heat scheme. For scheme (H) of Finite-difference approximations, with mesh ratio $s = \Delta t/(\Delta x)^{2}$ , the amplification factor is

\xi(\theta) \;=\; 1 - 2s\bigl(1 - \cos\theta\bigr) \;=\; 1 - 4s\,\sin^{2}\!\frac{\theta}{2},

which is real. Then

|\xi(\theta)| \le 1 \ \text{ for every } \theta \in [-\pi,\pi] \qquad\Longleftrightarrow\qquad s \le \tfrac12, \quad\text{that is}\quad \Delta t \le \tfrac12 (\Delta x)^{2} .

If $s > 1/2$ then $\xi(\pi) = 1 - 4s < -1$ , and the sawtooth mode is multiplied by the fixed factor $4s - 1 > 1$ at every step, so the von Neumann condition fails for every $C$ .

(b) The explicit wave scheme. For scheme (W) of Finite-difference approximations, with $\lambda = c\,\Delta t/\Delta x$ , the mode $u_j^n = \xi^n e^{ikj\Delta x}$ satisfies the quadratic

\xi^{2} - 2(1+p)\,\xi + 1 = 0, \qquad p \;=\; \lambda^{2}\bigl(\cos\theta - 1\bigr) \;\le\; 0 .

Both roots satisfy $|\xi| \le 1$ for every $\theta \in [-\pi,\pi]$ if and only if

\lambda \le 1, \quad\text{that is}\quad c\,\Delta t \le \Delta x .

Whenever $-2 < p < 0$ the two roots form a complex-conjugate pair of modulus exactly $1$ , so those modes neither grow nor decay, which is what the wave equation itself does; the endpoint $p = -2$ , reached only when $\lambda = 1$ and $\theta = \pi$ , gives the double root $\xi = -1$ . When $\lambda > 1$ there are wavenumbers for which one root is real and less than $-1$ .

Idea of the proofOriginal

A constant-coefficient difference scheme treats every grid point alike, so shifting the grid by one point commutes with one step of the scheme. The grid plane waves $e^{ikj\Delta x}$ are the eigenvectors of that shift, hence eigenvectors of the scheme, and the scheme acts on each of them by a single number $\xi$ . Stability is then the question of whether any of those numbers exceeds $1$ in modulus, and since the whole dependence is through $\theta = k\Delta x$ in a bounded interval, it is a question about a function of one real variable on $[-\pi,\pi]$ . For the heat scheme that function is a real cosine and the worst case is the shortest wave the grid supports; for the wave scheme each $\theta$ gives a quadratic with constant term $1$ , so its roots multiply to $1$ and the only two possibilities are a conjugate pair on the unit circle or a real pair straddling it.

ProofAdaptedAdapted from Strauss (2008), §8.2, Eqs. (9)-(10), pp. 206-207; Strauss (2008), §8.3, Eqs. (7)-(8), p. 213

(a) Step 1 (substitute). The solved form of (H) is $u_j^{n+1} = s(u_{j+1}^n + u_{j-1}^n) + (1-2s)u_j^n$ , so the coefficients are $a_{\pm1} = s$ and $a_0 = 1-2s$ . By the boxed formula,

\xi(\theta) = s\,e^{i\theta} + s\,e^{-i\theta} + 1 - 2s = 1 - 2s(1 - \cos\theta),

and the half-angle identity $1 - \cos\theta = 2\sin^{2}(\theta/2)$ gives the second form. In particular $\xi$ is real because the stencil is symmetric.

Step 2 (one inequality is free). Since $s > 0$ and $\sin^{2}(\theta/2) \ge 0$ , we have $\xi(\theta) \le 1$ for every $\theta$ . So the condition $|\xi| \le 1$ reduces to $\xi(\theta) \ge -1$ .

Step 3 (find the worst mode). $\xi(\theta) \ge -1$ says $4s\sin^{2}(\theta/2) \le 2$ . As $\theta$ runs over $[-\pi,\pi]$ , $\sin^{2}(\theta/2)$ attains its largest value $1$ at $\theta = \pm\pi$ and no larger value is available, this being exactly where aliasing cuts the range off. The inequality therefore holds across the whole range of $\theta$ exactly when $4s \le 2$ , which is the stated bound $s \le 1/2$ .

Step 4 (failure is genuine). If $s > 1/2$ , take $\theta = \pi$ . Then $\xi = 1 - 4s$ and $|\xi| = 4s - 1 > 1$ is a constant that does not shrink with $\Delta t$ , so no bound of the form $1 + C\Delta t$ can hold for small $\Delta t$ . Over a fixed time $T$ the number of steps is $n = T/\Delta t = T/(s(\Delta x)^2)$ , so the sawtooth is amplified by $(4s-1)^{T/(s(\Delta x)^2)}$ , which tends to infinity as the grid is refined.

(b) Step 5 (substitute into the three-level scheme). Insert $u_j^n = \xi^{n}e^{ikj\Delta x}$ into

\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{(\Delta t)^2} = c^2\,\frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{(\Delta x)^2}

and divide by $\xi^{n}e^{ikj\Delta x}/(\Delta t)^{2}$ , which is non-zero. Using $\lambda^2 = c^2(\Delta t)^2/(\Delta x)^2$ and $e^{i\theta} + e^{-i\theta} = 2\cos\theta$ ,

\xi - 2 + \xi^{-1} \;=\; \lambda^{2}\bigl(2\cos\theta - 2\bigr) \;=\; 2p .

Multiplying by $\xi$ and rearranging gives $\xi^2 - 2(1+p)\xi + 1 = 0$ . Note $p \le 0$ always, and $p \ge -2\lambda^2$ with equality at $\theta = \pi$ .

Step 6 (the roots multiply to one). The product of the two roots is the constant term, $1$ . The discriminant of the quadratic is $4(1+p)^2 - 4 = 4p(p+2)$ .

Step 7 (case $-2 \le p \le 0$ ). Then $p(p+2) \le 0$ , so the discriminant is $\le 0$ and the roots are a complex-conjugate pair (or the double root $\xi = -1$ when $p = -2$ , and $\xi = 1$ when $p=0$ ). A conjugate pair whose product is $1$ consists of two numbers of modulus $1$ , since $\xi\,\overline{\xi} = |\xi|^2 = 1$ . So $|\xi| = 1$ for both roots.

Step 8 (case $p < -2$ ). Then $p(p+2) > 0$ , so the roots are real and distinct. Their product is $1$ and their sum is $2(1+p) < -2$ , so both are negative and they are reciprocals of each other; two distinct reciprocal numbers cannot both have modulus $1$ , so one has modulus strictly greater than $1$ . The scheme is unstable.

Step 9 (translate into $\lambda$ ). Over $\theta \in [-\pi,\pi]$ the quantity $p$ is smallest at $\theta = \pi$ , where it equals $-2\lambda^{2}$ . Step 7 therefore covers every $\theta$ exactly when $-2\lambda^{2} \ge -2$ , which is the condition $\lambda \le 1$ . Should instead $\lambda$ exceed $1$ , the sawtooth $\theta = \pi$ has $p = -2\lambda^{2} < -2$ , and Step 8 produces a root outside the unit circle. $\square$

Definition(The Courant-Friedrichs-Lewy condition)AdaptedAdapted from Strauss (2008), §8.3, p. 214 and Figure 4

The domain of dependence of the exact solution of $u_{tt} = c^2u_{xx}$ at a point $(x,t)$ is the interval $[x - ct,\, x+ct]$ on the initial line: by d'Alembert's formula (The one-dimensional wave equation) the value there is determined by the data on that interval and by nothing outside it. The numerical domain of dependence of $u_j^{n}$ under an explicit scheme whose stencil reaches one point to each side is the set of initial grid points within $n$ steps, namely $[x_j - n\Delta x,\, x_j + n\Delta x]$ .

The CFL condition, named for Courant, Friedrichs and Lewy, is the requirement that the numerical domain of dependence contain the exact one. For scheme (W) at $t_n = n\Delta t$ it reads $c\,n\,\Delta t \le n\,\Delta x$ , that is

\lambda = \frac{c\,\Delta t}{\Delta x} \le 1 .

Strauss puts the same thing as a comparison of speeds: $\Delta x/\Delta t$ is the speed at which the scheme propagates information, and it must be at least the physical speed $c$ .

Theorem(CFL is necessary for convergence)AdaptedAdapted from Strauss (2008), §8.3, p. 214

Let an explicit scheme for $u_{tt} = c^2u_{xx}$ on the whole line compute $u_j^{n+1}$ from grid values at the two previous levels lying within one grid point of $x_j$ , and refine the grid with $\lambda = c\,\Delta t/\Delta x$ held fixed. If $\lambda > 1$ , there is a pair of initial data for which the computed values at a fixed point $(x_\ast, t_\ast)$ , $t_\ast > 0$ , do not both converge to the corresponding exact solutions. Hence no such scheme can converge for all data unless $\lambda \le 1$ , whatever its truncation error.

ProofAdaptedAdapted from Strauss (2008), §8.3, p. 214

Step 1 (the numerical domain of dependence is too small). Suppose $\lambda > 1$ . Fix $t_\ast > 0$ and take the point $x_\ast = 0$ . Let $n$ be the number of steps to reach $t_\ast$ , so $n\,\Delta t = t_\ast$ . Each of the $n$ steps widens the stencil by one grid point, and the second starting row widens it by one more, so the computed value at $(0, t_\ast)$ depends only on the initial data at the grid points of $[-(n+1)\Delta x,\, (n+1)\Delta x]$ . Now

n\,\Delta x = n\,\frac{c\,\Delta t}{\lambda} = \frac{c\,t_\ast}{\lambda} \;=:\; r \;<\; c\,t_\ast ,

so that interval is contained in $[-(r+\Delta x),\, r+\Delta x]$ , which shrinks to the fixed interval $[-r,r]$ as the grid is refined with $\lambda$ held fixed. It is therefore eventually contained in any fixed neighbourhood of $[-r,r]$ inside the exact domain of dependence $[-ct_\ast, ct_\ast]$ .

Step 2 (choose the data). Pick $r'$ with $r < r' < ct_\ast$ . Take initial velocity $\psi \equiv 0$ in both cases, initial displacement $\varphi_1 \equiv 0$ , and for $\varphi_2$ a continuous bump supported in $(r', \infty)$ with $\varphi_2(ct_\ast) = 2$ . Since $\psi \equiv 0$ , d'Alembert's formula gives $u(0,t_\ast) = \tfrac12\bigl(\varphi(ct_\ast) + \varphi(-ct_\ast)\bigr)$ , so the exact values are $u_1(0,t_\ast) = 0$ and $u_2(0,t_\ast) = \tfrac12(2 + 0) = 1$ .

Step 3 (compare). Take any grid fine enough that $r + \Delta x < r'$ . On it, $\varphi_1$ and $\varphi_2$ agree (both vanish) at every grid point of $[-(r+\Delta x), r+\Delta x]$ , and $\psi$ is the same for both; by Step 1 the computed value at $(0,t_\ast)$ uses nothing else, so the two computations produce the same number. Their exact counterparts differ by $1$ . Hence at most one of the two numerical sequences can converge to its own exact value as the grid is refined, and convergence for all data fails. $\square$

Remark(Where the hypotheses are used)Original

Constant coefficients. The entire method rests on the grid plane waves being eigenvectors of the scheme, which needs the coefficients $a_m$ to be independent of $j$ . For a variable-coefficient problem the usual practice is to freeze the coefficients at each point, apply the analysis there, and impose the worst of the resulting limits; that is a heuristic, not a theorem.

No boundaries. The modes $e^{ikj\Delta x}$ live on the whole line or on a periodic grid. Strauss gets the same answer on the interval $[0,\pi]$ with Dirichlet data because the discrete sine functions $\sin(jk\Delta x)$ happen also to be eigenvectors of the scheme there; for a general boundary condition they are not, and a scheme that is von Neumann stable in the interior can still be destroyed by its boundary treatment.

Scalar, one step. For a scalar one-step scheme, $|\xi(\theta)| \le 1$ for all $\theta$ is both necessary and sufficient for the scheme to be a contraction in the grid $\ell^2$ norm, by the discrete Parseval relation of The discrete Fourier transform and the FFT. For a system, the amplification is a matrix $G(\theta)$ , and the spectral radius condition $\rho(G) \le 1$ is still necessary but no longer sufficient: powers of a non-normal matrix with spectral radius $1$ can grow polynomially, or transiently by a large factor. The multi-level scheme (W) was handled above by passing to the scalar quadratic, which is legitimate only because the two roots were examined individually. The same caution applies at the endpoint $\lambda = 1$ : there the sawtooth $\theta = \pi$ has the double root $\xi = -1$ , so the general solution of the recurrence for that mode is $(A + Bn)(-1)^{n}$ and grows linearly in $n$ , even though both roots have modulus $1$ . Starting rows built exactly from a d'Alembert solution have $B = 0$ , which is why the scheme is exact at $\lambda = 1$ (Finite-difference approximations); rounding error is not so obliging, and practice therefore keeps $\lambda$ strictly below $1$ .

Where the aliasing cut-off is used. Steps 3 and 9 both take the supremum of the relevant quantity over $\theta \in [-\pi,\pi]$ and over nothing larger. A continuum stability analysis, which would let $k$ run to infinity, would give no condition at all for the heat equation and the wrong one for the wave equation. The limits $s \le 1/2$ and $\lambda \le 1$ are properties of the grid, and they come from the shortest wave the grid can represent.

The slack $C\,\Delta t$ . For $u_t = u_{xx}$ the exact solution decays, so $|\xi| \le 1$ is the right test. For an equation such as $u_t = a u_{xx} + bu$ with $b > 0$ the exact solution grows like $e^{bt}$ , and a scheme that reproduced it would have $|\xi| \approx 1 + b\Delta t > 1$ ; insisting on $|\xi| \le 1$ would reject a perfectly good scheme. That is why Strauss states the condition with the $O(\Delta t)$ allowance.

Implicit schemes escape the limit. Strauss's $\theta$ scheme weights the second difference between the two time levels; writing the weight as $\mu \in [0,1]$ to keep $\theta$ for the phase, its amplification factor is

\xi(\theta) = \frac{1 - 2(1-\mu)s(1-\cos\theta)}{1 + 2\mu s(1-\cos\theta)} ,

and he shows that for $\tfrac12 \le \mu \le 1$ no restriction on $s$ is needed: the scheme is unconditionally stable, the case $\mu = \tfrac12$ being Crank-Nicolson. The price is that each step now requires the solution of a tridiagonal linear system. Whether that price is worth paying is the stiffness question of Stability and stiff equations.

Theorem (stated without proof)(Lax equivalence theorem)Original

Let a linear initial-value problem be well posed on a Banach space $B$ : for data $u_0$ in a dense subspace there is a unique solution $u(t) = E(t)u_0$ for $0 \le t \le T$ , and the solution operators $E(t)$ are uniformly bounded on $[0,T]$ . Let a difference scheme be written as $u^{n+1} = C(\Delta t)\,u^{n}$ with $C(\Delta t)$ a bounded operator on $B$ , and suppose the scheme is consistent with the problem. Call the scheme stable if the family

\bigl\{\,C(\Delta t)^{n} \;:\; 0 \le n\,\Delta t \le T,\ 0 < \Delta t \le \tau \,\bigr\}

is uniformly bounded in operator norm for some $\tau > 0$ . Then the scheme is convergent — $C(\Delta t)^{n}u_0 \to E(t)u_0$ for every $u_0 \in B$ as $\Delta t \to 0$ with $n\,\Delta t \to t$ — if and only if it is stable.

This is the theorem of Lax and Richtmyer. No book in this handbook's library proves it, and none is cited for it; it is stated here because it is what licenses the whole practice above, and it is the reason the von Neumann computation is worth doing at all. Strauss refers the reader to Richtmyer and Morton for the same circle of ideas when he states the von Neumann condition, and explicitly omits the convergence proof for the heat scheme.

NoteOriginal

How the three notions fit together, for a constant-coefficient scheme on a periodic grid. Consistency is checked by Taylor expansion (Finite-difference approximations) and says the scheme solves the right equation. Stability is checked by the amplification factor and says the scheme does not amplify anything without limit. Convergence — the thing one actually wants — follows from the two together, by Lax equivalence, and from neither alone. The mesh ratio $s = 1/6$ found in the previous unit raises the order of consistency and says nothing about stability; the condition $s \le 1/2$ found here constrains stability and says nothing about accuracy.

Method(Von Neumann analysis of a scheme)Original

Check the scheme is constant-coefficient and write it with all time levels on the left. Drop the boundaries: work on the whole line or a periodic grid.
Substitute $u_j^{n} = \xi^{n}e^{i\theta j}$ and divide out the common factor $\xi^{n}e^{i\theta j}$ (and a power of $\xi$ if the scheme is multi-level). What remains is a polynomial equation in $\xi$ whose coefficients depend on $\theta$ and on the mesh ratios.
Solve for $\xi$ . A one-step scheme gives $\xi$ directly; an $r$ -level scheme gives a polynomial of degree $r-1$ in $\xi$ with that many roots, every one of which must be tested.
Simplify with half angles. Replace $1 - \cos\theta$ by $2\sin^{2}(\theta/2)$ and $\sin\theta$ by $2\sin(\theta/2)\cos(\theta/2)$ ; the modulus is usually a quadratic in $\sigma = \sin^{2}(\theta/2) \in [0,1]$ .
Maximise $|\xi|$ over $\sigma \in [0,1]$ , not over all real $k$ . The maximum is usually at an endpoint: $\sigma = 1$ is the sawtooth $\theta = \pi$ and $\sigma = 0$ is the constant mode.
Impose $|\xi| \le 1$ (or $\le 1 + C\Delta t$ if the exact solution grows) and solve the resulting inequality for the mesh ratio. If it holds for every mesh ratio, the scheme is unconditionally stable.
Confirm numerically. Run the scheme on the offending mode for a few dozen steps and check the measured growth against $|\xi|^{n}$ .

Example(Three schemes for the advection equation)Original

Take $u_t + a u_x = 0$ with $a > 0$ , whose exact solution translates the initial profile rigidly to the right at speed $a$ (First-order equations and characteristics). Put $\nu = a\,\Delta t/\Delta x$ , the Courant number, and $\theta = k\,\Delta x$ .

Scheme 1: forward in time, centred in space. $u_j^{n+1} = u_j^{n} - \tfrac{\nu}{2}\bigl(u_{j+1}^n - u_{j-1}^n\bigr)$ . The coefficients are $a_0 = 1$ , $a_{\pm 1} = \mp \nu/2$ , so

\xi(\theta) = 1 - \frac{\nu}{2}\bigl(e^{i\theta} - e^{-i\theta}\bigr) = 1 - i\nu\sin\theta, \qquad |\xi|^{2} = 1 + \nu^{2}\sin^{2}\theta .

This exceeds $1$ for every $\theta$ that is not a multiple of $\pi$ , and the excess does not shrink with $\Delta t$ at fixed $\nu$ . The scheme is unstable for every choice of steps, even though its truncation error is $O(\Delta t) + O((\Delta x)^2)$ and it looks like the natural thing to write down.

Scheme 2: Lax-Friedrichs. Replace the $u_j^n$ on the right by the average of its neighbours: $u_j^{n+1} = \tfrac12(u_{j+1}^n + u_{j-1}^n) - \tfrac{\nu}{2}(u_{j+1}^n - u_{j-1}^n)$ . Now $a_{\pm1} = \tfrac12 \mp \tfrac{\nu}{2}$ , $a_0 = 0$ , so

\xi(\theta) = \cos\theta - i\nu\sin\theta, \qquad |\xi|^{2} = \cos^{2}\theta + \nu^{2}\sin^{2}\theta = 1 - (1-\nu^{2})\sin^{2}\theta .

So $|\xi| \le 1$ for all $\theta$ if and only if $\nu^{2} \le 1$ : the CFL condition again.

Scheme 3: upwind. Difference into the wind, $u_j^{n+1} = u_j^n - \nu\bigl(u_j^n - u_{j-1}^n\bigr)$ . Then $\xi = 1 - \nu + \nu e^{-i\theta}$ and

|\xi|^{2} = (1-\nu+\nu\cos\theta)^{2} + \nu^{2}\sin^{2}\theta = 1 - 2\nu(1-\nu)(1-\cos\theta),

after expanding and using $\cos^2 + \sin^2 = 1$ . Since $1 - \cos\theta \ge 0$ , this is at most $1$ precisely when $\nu(1-\nu) \ge 0$ , that is $0 \le \nu \le 1$ . (For $a < 0$ the mirror-image scheme, differencing forward in $x$ , is the stable one.)

The numbers. Take $\nu = \tfrac12$ and the mode with $\theta = \pi/2$ , that is four grid points per wavelength. On a periodic grid of $64$ points with $\Delta x = 2\pi/64$ this is the mode $e^{16 i x_j}$ .

Scheme 1: $|\xi| = \sqrt{1 + \tfrac14} = 1.11803399$ , so after $100$ steps the mode should be multiplied by $(5/4)^{50} = 7.00649232\times10^{4}$ . Marching the scheme $100$ steps in double precision and dividing by the initial mode gives $7.00649232\times10^{4}$ , agreeing in all nine figures.
Scheme 2: $|\xi|^2 = 0 + \tfrac14$ , so $|\xi| = \tfrac12$ and after $20$ steps the factor should be $2^{-20} = 9.5367431640625\times10^{-7}$ . The computed factor is $9.5367431640625\times10^{-7}$ , exact to the last bit, since $2^{-20}$ is a machine number.
Scheme 3: $|\xi|^{2} = 1 - 2(\tfrac12)(\tfrac12)(1 - 0) = \tfrac12$ , so after $20$ steps the factor should be $2^{-10} = 9.765625\times10^{-4}$ . The computed factor is $9.76562500\times10^{-4}$ .

Schemes 2 and 3 are stable here but heavily damped at this wavenumber, which is the other half of the story: stability prevents disaster, it does not deliver accuracy. Over the $20$ steps the exact solution merely translates the wave, losing no amplitude at all, while the computed one has lost three orders of magnitude of it. This numerical dissipation is the price of the averaging in Scheme 2 and of the one-sided difference in Scheme 3, and it is the reason schemes of higher order are used in practice.

And the heat scheme. The same experiment on (H) with $\theta = \pi$ : at $s = 0.55$ , just past the limit, $\xi = 1 - 4(0.55) = -1.2$ , and after $100$ steps the sawtooth is multiplied by $1.2^{100} = 8.28\times10^{7}$ — a complete loss of the solution, from a mesh ratio only ten per cent too large. At $s = 0.5$ exactly, $\xi = -1$ and the sawtooth survives forever with alternating sign; at $s = 0.25$ , $\xi = 0$ and it is annihilated in one step. Strauss reports exactly this behaviour in two runs with mesh ratios straddling $1/2$ : one on target, one oscillating wildly.

ApplicationOriginal

The CFL condition is the single hardest constraint on large-scale geophysical simulation, and it is a budget item rather than a mathematical curiosity.

A global atmospheric model that treats the full compressible equations explicitly must resolve sound waves, which travel at roughly $340\ \mathrm{m\,s^{-1}}$ . At a horizontal grid spacing of $10\ \mathrm{km}$ the CFL limit is about $\Delta t \le 10^{4}/340 \approx 29\ \mathrm{s}$ , so a one-year integration needs of the order of $10^{6}$ time steps. The meteorologically interesting signals travel at a tenth of that speed, so almost all of the work is spent keeping the acoustic modes stable. This is precisely the stiffness of Stability and stiff equations, and the standard answer is the one the $\theta$ scheme illustrates: treat the fast terms implicitly, pay for a linear solve each step, and choose $\Delta t$ from accuracy rather than from stability. Semi-implicit and split-explicit dynamical cores exist for exactly this reason.

The parabolic limit bites even harder. A vertical diffusion scheme of the form (H) with a diffusivity $\kappa$ requires $\Delta t \le (\Delta z)^2/(2\kappa)$ , which is quadratic in the spacing: refining the vertical grid by a factor of ten costs a hundred times as many steps. Boundary-layer schemes in operational models are therefore implicit almost without exception.

The same arithmetic governs explicit finite-difference codes for seismic wave propagation, for electromagnetic simulation by the finite-difference time-domain method, and for the diffusion step of a reactor or heat-exchanger calculation. In each case the engineering question is not whether the scheme is stable but how much the stability limit costs, which is why unconditionally stable implicit schemes and the iterative solvers of Relaxation methods for elliptic problems are so central to practice.

Method

Relaxation methods for elliptic problems

Adapted

Add to PDF

Sets up the five-point discrete Poisson problem, writes the Jacobi, Gauss-Seidel and successive over-relaxation sweeps for it, computes the spectral radius of each error propagator exactly on the unit square ( $\cos \pi h$ , $\cos^2 \pi h$ and $(1-\sin\pi h)/(1+\sin\pi h)$ ), and turns those into iteration counts and operation counts that are checked against a run.

Prerequisites: Finite-difference approximations, Solving linear systems numerically

Definition(The discrete Poisson problem)AdaptedAdapted from Strauss (2008), §8.4, Eqs. (1)-(2), pp. 218-219

Let $\Omega = (0,1)^{2}$ , let $h = 1/(m+1)$ , and place the grid points $(x_j, y_k) = (jh, kh)$ for $0 \le j,k \le m+1$ . The $N = m^{2}$ points with $1 \le j,k \le m$ are interior; the rest lie on the boundary. We solve

-\,u_{xx} - u_{yy} = f \ \text{ in } \Omega, \qquad u = 0 \ \text{ on } \partial\Omega

by replacing each second derivative with the centred second difference of Finite-difference approximations. Writing $U_{j,k}$ for the approximation to $u(x_j,y_k)$ and taking equal steps in both directions, the five-point scheme is

\frac{4U_{j,k} - U_{j+1,k} - U_{j-1,k} - U_{j,k+1} - U_{j,k-1}}{h^{2}} \;=\; f_{j,k}, \qquad 1 \le j,k \le m, \tag{P}

with $U_{j,k} = 0$ whenever $(j,k)$ is a boundary index. Collecting one equation per interior node turns (P) into the matrix problem $A\,U = F$ , in which $A$ holds the left-hand coefficients, $F$ holds the values $f_{j,k}$ , and both have size $N$ .

For $f \equiv 0$ , dividing (P) by $4/h^2$ puts it in the form Strauss uses,

U_{j,k} \;=\; \tfrac14\bigl(U_{j+1,k} + U_{j-1,k} + U_{j,k+1} + U_{j,k-1}\bigr),

the statement that each value is the average of its four neighbours. This is the exact discrete counterpart of the mean value property of harmonic functions (Harmonic functions), and it forces a discrete maximum principle: a non-constant solution cannot attain its largest or smallest value at an interior point, since such a point could not be the average of its neighbours. Strauss derives the uniqueness of the solution of (P) from that principle, and then nonsingularity of $A$ from uniqueness.

The decisive structural point is that there is no marching direction. For the heat and wave equations the new row was computed from the old; an elliptic problem couples all $N$ unknowns at once, and $N$ is large — a grid of $100$ intervals in each direction already gives $N \approx 10^{4}$ . Direct elimination is possible (Solving linear systems numerically) but its cost grows faster than linearly in $N$ , which is what motivates the iterations below.

Definition(Splitting, and the three sweeps)AdaptedAdapted from Stakgold & Holst (2011), §10.6, Eqs. (10.6.28)-(10.6.30), p. 799; Strauss (2008), §8.4, Eqs. (3)-(5), pp. 220-221

Split $A = D - L - U$ , with $D$ the diagonal part, $-L$ the strictly lower triangular part and $-U$ the strictly upper triangular part, the ordering of the unknowns being fixed once and for all. For the five-point matrix, $D = (4/h^{2})I$ . An iteration of the form

U^{(n+1)} \;=\; U^{(n)} + B\bigl(F - A\,U^{(n)}\bigr)

for an approximate inverse $B$ of $A$ has error propagation matrix $E = I - BA$ : subtracting the identity $U^{\ast} = U^{\ast} + B(F - AU^{\ast})$ satisfied by the exact discrete solution gives $U^{(n+1)} - U^{\ast} = E\,(U^{(n)} - U^{\ast})$ , so the error after $n$ sweeps is $E^{n}$ applied to the initial error. The three classical choices, in Stakgold and Holst's notation, are

B_{\mathrm{J}} = D^{-1}, \qquad B_{\mathrm{GS}} = (D-L)^{-1}, \qquad B_{\omega} = \omega\,(D - \omega L)^{-1} ,

giving Jacobi, Gauss-Seidel and successive over-relaxation (SOR) with relaxation factor $\omega$ . On the grid, with $g_{j,k} = h^{2}f_{j,k}$ , these read

Jacobi. $\displaystyle U^{(n+1)}_{j,k} = \tfrac14\bigl(U^{(n)}_{j+1,k} + U^{(n)}_{j-1,k} + U^{(n)}_{j,k+1} + U^{(n)}_{j,k-1} + g_{j,k}\bigr)$ : every value on the right is taken from the previous sweep, so the order of the nodes does not matter and the whole old array must be kept.
Gauss-Seidel. The same formula, but sweeping the nodes in a fixed order (say row by row from the bottom, left to right within a row) and using each new value the moment it is available. With the lexicographic order this uses the updated $U_{j-1,k}$ and $U_{j,k-1}$ and the old $U_{j+1,k}$ and $U_{j,k+1}$ . No second array is needed.
SOR. Compute the Gauss-Seidel value $\widetilde U_{j,k}$ and then move past it:

U^{(n+1)}_{j,k} \;=\; U^{(n)}_{j,k} + \omega\bigl(\widetilde U_{j,k} - U^{(n)}_{j,k}\bigr) .

$\omega = 1$ is Gauss-Seidel; $\omega > 1$ is over-relaxation.

A convention to watch. Strauss writes the SOR step as $U^{(n+1)}_{j,k} = U^{(n)}_{j,k} + \omega\bigl(\text{sum of the four neighbours} - 4U^{(n)}_{j,k}\bigr)$ , in which Gauss-Seidel is the case $\omega = 1/4$ . His $\omega$ is therefore one quarter of the $\omega$ used here and by Stakgold and Holst. Everything below uses the convention in which Gauss-Seidel is $\omega = 1$ ; a factor quoted from Strauss must be multiplied by four before it is substituted into any formula here.

Method(Convergence, and the rates for the model problem)AdaptedAdapted from Stakgold & Holst (2011), §10.6, Thm 10.6.4 and the spectral radii, pp. 799-800

Let $E$ be the error propagation matrix of one of the three iterations above, applied to the five-point matrix $A$ of (P) on the unit square with $h = 1/(m+1)$ .

Convergence. $U^{(n)} \to U^{\ast}$ for every starting vector if and only if the spectral radius satisfies $\rho(E) < 1$ , and then the error is reduced asymptotically by the factor $\rho(E)$ per sweep. Stakgold and Holst record that for an M-matrix $A$ the splitting obtained by setting the off-diagonal entries to zero is regular and the resulting iteration converges, and that Jacobi and Gauss-Seidel both arise from regular splittings of such a matrix.
Jacobi. The eigenvalues of $E_{\mathrm{J}} = I - D^{-1}A$ are

\mu_{p,q} = \tfrac12\bigl(\cos p\pi h + \cos q\pi h\bigr), \qquad 1 \le p,q \le m,

so $\rho(E_{\mathrm{J}}) = \cos \pi h$ .

Gauss-Seidel. For a consistently ordered matrix with Property A — which the five-point matrix with the lexicographic ordering is — the Gauss-Seidel spectral radius is the square of the Jacobi one, so $\rho(E_{\mathrm{GS}}) = \cos^{2}\pi h$ .
SOR. $\rho(E_{\omega}) \ge |1-\omega|$ for every matrix with non-zero diagonal, so convergence requires $0 < \omega < 2$ . For this model problem the optimal factor and the rate it achieves are

\omega_{\mathrm{opt}} = \frac{2}{1 + \sin\pi h}, \qquad \rho\bigl(E_{\omega_{\mathrm{opt}}}\bigr) = \omega_{\mathrm{opt}} - 1 = \frac{1 - \sin \pi h}{1 + \sin\pi h} .

Rates as $h \to 0$ . Expanding the logarithms,

-\ln\rho(E_{\mathrm{J}}) \sim \tfrac12\pi^{2}h^{2}, \qquad -\ln\rho(E_{\mathrm{GS}}) \sim \pi^{2}h^{2}, \qquad -\ln\rho(E_{\omega_{\mathrm{opt}}}) \sim 2\pi h .

Gauss-Seidel is asymptotically twice as fast as Jacobi; optimal SOR is faster than Gauss-Seidel by the factor $2/(\pi h)$ , which grows without bound as the grid is refined. All three rates tend to zero with $h$ : Stakgold and Holst emphasise that $\rho(E) \to 1$ in every case, so every classical method slows down as the problem is refined.

Idea of the proofOriginal

Two observations do all the work. First, the five-point matrix on a rectangle has explicitly known eigenvectors — the products of discrete sines, the same ones that separation of variables produces for the continuous problem — so the Jacobi error propagator, which is a polynomial in $A$ , is diagonalised at once and its spectral radius can simply be read off. The slowest mode is the smoothest one, $p = q = 1$ , whose eigenvalue is only $O(h^{2})$ below $1$ : Jacobi averages neighbours, and averaging barely changes a function that is already nearly constant over a few grid spacings. Second, Gauss-Seidel and SOR are not polynomials in $A$ and cannot be diagonalised this way, but SOR's determinant can be computed from triangularity alone, and that single determinant already forces $\omega$ into $(0,2)$ and sets a floor $|1-\omega|$ under the spectral radius. Over-relaxation works by pushing $\omega$ up until the floor $|1-\omega|$ rises to meet the falling spectral radius, and the two meet at $\omega_{\mathrm{opt}}$ .

ProofAdaptedAdapted from Stakgold & Holst (2011), §10.6, p. 800

Part 1 (convergence criterion). The error satisfies $e^{(n)} = E^{n}e^{(0)}$ . A standard fact about matrix powers is that $E^{n} \to 0$ for every starting vector if and only if $\rho(E) < 1$ , and that $\|E^{n}\|^{1/n} \to \rho(E)$ ; see Eigenvalues and eigenvectors and Solving linear systems numerically. The M-matrix statement is Stakgold and Holst's Theorem 10.6.4, reported here as they give it.

Part 2 (Jacobi eigenvalues). For $1 \le p,q \le m$ define the grid function

v^{(p,q)}_{j,k} \;=\; \sin(p\pi j h)\,\sin(q\pi k h) .

It vanishes at $j = 0$ and at $j = m+1$ , because $\sin\bigl(p\pi(m+1)h\bigr) = \sin(p\pi) = 0$ , and likewise in $k$ ; so it respects the homogeneous boundary condition built into (P). Using $\sin\bigl((j\pm1)\alpha\bigr) = \sin j\alpha\cos\alpha \pm \cos j\alpha\sin\alpha$ , the two terms with $\cos j\alpha$ cancel on addition and

v^{(p,q)}_{j+1,k} + v^{(p,q)}_{j-1,k} = 2\cos(p\pi h)\, v^{(p,q)}_{j,k},

with the same identity in $k$ . Hence

A\,v^{(p,q)} = \lambda_{p,q}\,v^{(p,q)}, \qquad \lambda_{p,q} = \frac{4 - 2\cos p\pi h - 2\cos q\pi h}{h^{2}} .

These $m^{2} = N$ grid functions are linearly independent (they are the columns of a tensor product of two discrete sine transforms), so they are a full eigenbasis. Since $D = (4/h^{2})I$ , the Jacobi propagator is $E_{\mathrm{J}} = I - (h^{2}/4)A$ , a polynomial in $A$ with the same eigenvectors, and its eigenvalues are

1 - \frac{h^{2}}{4}\lambda_{p,q} = 1 - \frac{4 - 2\cos p\pi h - 2\cos q\pi h}{4} = \tfrac12\bigl(\cos p\pi h + \cos q\pi h\bigr),

which is claim 2. The largest in modulus occurs at $p = q = 1$ , giving $\cos\pi h$ , and at $p = q = m$ , giving $\cos(\pi - \pi h) = -\cos\pi h$ ; both have modulus $\cos \pi h$ , and every other pair gives something smaller. So $\rho(E_{\mathrm{J}}) = \cos\pi h < 1$ , and the iteration converges.

Part 4a (the floor under the SOR spectral radius). From $B_\omega = \omega(D-\omega L)^{-1}$ and $A = D - L - U$ ,

E_{\omega} = I - \omega(D-\omega L)^{-1}A = (D-\omega L)^{-1}\bigl[(D - \omega L) - \omega(D - L - U)\bigr] = (D-\omega L)^{-1}\bigl[(1-\omega)D + \omega U\bigr].

The matrix $D - \omega L$ is lower triangular with diagonal $D$ , and $(1-\omega)D + \omega U$ is upper triangular with diagonal $(1-\omega)D$ . The determinant of a triangular matrix is the product of its diagonal entries, so

\det E_{\omega} = \frac{(1-\omega)^{N}\det D}{\det D} = (1-\omega)^{N} .

The determinant is also the product of the $N$ eigenvalues, each of modulus at most $\rho(E_\omega)$ , so $\rho(E_{\omega})^{N} \ge |1-\omega|^{N}$ , that is $\rho(E_{\omega}) \ge |1-\omega|$ . Convergence needs $\rho < 1$ , hence $|1-\omega| < 1$ , hence $0 < \omega < 2$ .

Part 5 (the rates). For small $h$ , $\cos\pi h = 1 - \tfrac12\pi^{2}h^{2} + O(h^{4})$ , so $-\ln\cos\pi h = \tfrac12\pi^{2}h^{2} + O(h^{4})$ and the Gauss-Seidel rate is twice that, since the logarithm of a square is twice the logarithm. For SOR write $\sigma = \sin\pi h$ ; then

-\ln\frac{1-\sigma}{1+\sigma} = \ln\frac{1+\sigma}{1-\sigma} = 2\Bigl(\sigma + \frac{\sigma^{3}}{3} + \cdots\Bigr) = 2\pi h + O(h^{3}) . \qquad\square

Remark(Where the hypotheses are used, and what is not proved here)Original

The eigenvector computation needs the rectangle and the uniform grid. Part 2 is exact only because the domain is a product of intervals, the grid is uniform, the stencil is the same at every interior node, and the boundary condition is homogeneous Dirichlet. On an L-shaped domain, or with a variable coefficient, the sines are no longer eigenvectors and the spectral radius has to be estimated rather than computed. The qualitative conclusion $\rho = 1 - O(h^2)$ survives, because it is driven by the smallest eigenvalue of $A$ , which is $O(1)$ while the largest is $O(h^{-2})$ : the condition number of the discrete Laplacian is $O(h^{-2})$ whatever the domain, and that is the real culprit.

Gauss-Seidel depends on the ordering, Jacobi does not. Permuting the unknowns changes $L$ and $U$ and so changes $E_{\mathrm{GS}}$ ; it leaves $E_{\mathrm{J}} = I - D^{-1}A$ similar to itself. Claim 3 is stated for a consistently ordered matrix with Property A, which is the hypothesis Stakgold and Holst attach to it; with a red-black ordering of the grid the identity still holds, but with an arbitrary permutation it need not. This is also why Strauss is careful to prescribe the sweep order when he writes the Gauss-Seidel formula.

Claim 4 is only half proved here. The floor $\rho(E_\omega) \ge |1-\omega|$ , and hence the necessity of $0 < \omega < 2$ , is proved above from the determinant alone and holds for any matrix with non-zero diagonal. The value of $\omega_{\mathrm{opt}}$ is a separate classical result of Young's, which requires the consistent-ordering hypothesis and relates the SOR eigenvalues $\xi$ to the Jacobi eigenvalues $\mu$ through $(\xi + \omega - 1)^{2} = \xi\,\omega^{2}\mu^{2}$ . No book in this handbook's library proves it; Stakgold and Holst state only that $\rho(E_{\mathrm{SOR}}) \approx 1 - O(h)$ , and Strauss says explicitly that he leaves the choice of relaxation factor to specialist texts. The formula is quoted here because it is what practitioners use, and it is verified numerically in the example below rather than taken on trust.

Spectral radius is an asymptotic statement. $\rho(E)$ governs the error reduction per sweep eventually. $E_{\mathrm{J}}$ is symmetric for this problem, so for Jacobi the bound bites immediately; $E_{\mathrm{GS}}$ and $E_{\omega}$ are not symmetric, and for SOR near $\omega_{\mathrm{opt}}$ the eigenvalues are nearly defective, so there is a transient during which the error falls more slowly than $\rho^{n}$ . The example below shows exactly that: Jacobi and Gauss-Seidel take the predicted number of sweeps, SOR takes about a quarter more.

Jacobi is a disguised heat-equation solve. Strauss points out that one Jacobi sweep is one step of the explicit scheme for $v_t = v_{xx} + v_{yy}$ with $\Delta x = \Delta y = h$ and $\Delta t = h^{2}/4$ , so that relaxing to a steady state is integrating the heat equation to $t = \infty$ . The two-dimensional von Neumann analysis of Stability, the CFL condition and von Neumann analysis gives the amplification factor $1 - 4s_1\sin^{2}(\theta_1/2) - 4s_2\sin^{2}(\theta_2/2)$ with $s_1 = s_2 = 1/4$ , and on the Dirichlet grid the admissible phases are $\theta_i = p\pi h$ with $1 \le p \le m$ . At $p = q = 1$ this is $1 - 2\sin^{2}(\pi h/2) = \cos\pi h$ , which is exactly $\rho(E_{\mathrm{J}})$ from Part 2 — the two calculations are the same calculation. It also shows that the Jacobi sweep sits precisely at the two-dimensional stability limit $s_1 + s_2 = 1/2$ , so it is the largest explicit step available: Jacobi is slow not because the step is timid but because the heat equation itself takes a time of order $1$ to reach equilibrium, needing $O(h^{-2})$ steps of size $h^2/4$ .

Why the classical methods stall, and what replaces them. Stakgold and Holst show that a Gauss-Seidel sweep damps the high-frequency components of the error quickly while the smooth components barely move — which is just the statement that $\mu_{p,q}$ is near $1$ for small $p,q$ . Representing the smooth error on a coarser grid, where it is no longer smooth, and correcting there, is the multigrid idea, and it removes the $h$ -dependence of the rate altogether.

Method(Solving a discrete elliptic problem by relaxation)Original

Discretise on a uniform grid with the five-point stencil, and fold the boundary values into the right-hand side, so that the unknowns are the interior values only.
Scale the equation so that the diagonal is $1$ after dividing: each node's update is the average of its four neighbours plus $h^{2}f_{j,k}/4$ .
Pick the sweep. Gauss-Seidel dominates Jacobi in both speed and storage and should be preferred unless the sweeps are to be run in parallel, where Jacobi (or a red-black Gauss-Seidel) is the natural choice because all nodes of one colour are independent.
Choose $\omega$ . For a rectangle use $\omega_{\mathrm{opt}} = 2/(1+\sin\pi h)$ . Otherwise estimate $\rho(E_{\mathrm{J}})$ by monitoring the ratio of successive Jacobi residual norms, then set $\omega = 2/\bigl(1 + \sqrt{1 - \rho(E_{\mathrm{J}})^{2}}\bigr)$ . Erring on the high side of $\omega_{\mathrm{opt}}$ is far safer than erring low: above the optimum the spectral radius is exactly $\omega - 1$ , so it rises with unit slope, while below it the rise is much steeper. On the grid of the example, $\omega = 1.50$ gives $\rho = 0.500$ and $\omega = 1.60$ gives $\rho = 0.600$ , but dropping only to $\omega = 1.44$ already gives $\rho = 0.511$ and $\omega = 1.40$ gives $\rho = 0.611$ .
Stop on the residual, not on the step: iterate until $\|F - AU^{(n)}\|_\infty$ falls below a tolerance set by the discretisation error, which for the five-point scheme is $O(h^{2})$ . Iterating far past that is wasted work, since the discrete solution is not the one wanted anyway.
Predict the cost before running. Sweeps needed $\approx |\ln\varepsilon| / \bigl(-\ln\rho(E)\bigr)$ ; work per sweep $\approx 5N$ floating-point operations for Jacobi or Gauss-Seidel and $\approx 8N$ for SOR, with $N$ the number of interior nodes.
Check the answer against the grid refinement, not against the iteration count: halve $h$ and confirm the solution changes by a factor of four less than the previous change, which is the second-order convergence of the discretisation.

Example(A problem with an exactly known discrete solution)Original

Solve $-u_{xx}-u_{yy} = f$ on $(0,1)^{2}$ with $u = 0$ on the boundary and

f(x,y) = 2\pi^{2}\sin\pi x\,\sin\pi y ,

whose exact solution is $u(x,y) = \sin\pi x\,\sin\pi y$ , with $u(\tfrac12,\tfrac12) = 1$ .

Step 1 (the discrete solution in closed form). The grid values of $f$ are $2\pi^{2}v^{(1,1)}_{j,k}$ , a multiple of the first eigenvector of Part 2 of the proof. So the discrete solution is the same multiple of that eigenvector divided by the eigenvalue:

U^{\ast} = \frac{2\pi^{2}}{\lambda_{1,1}}\,v^{(1,1)}, \qquad \lambda_{1,1} = \frac{4\bigl(1 - \cos\pi h\bigr)}{h^{2}} .

At the centre, where $v^{(1,1)} = 1$ , this gives

$m$	$h$	$\lambda_{1,1}$	$U^{\ast}_{\text{centre}}$	error	ratio
$3$	$1/4$	$18.7451660$	$1.05302929$	$5.302929\times10^{-2}$
$7$	$1/8$	$19.4868397$	$1.01295075$	$1.295075\times10^{-2}$	$4.0947$
$15$	$1/16$	$19.6758729$	$1.00321896$	$3.218964\times10^{-3}$	$4.0233$

The ratios tend to $4$ : the discretisation is second order, as Finite-difference approximations says it should be. Note that $\lambda_{1,1} \to 2\pi^{2} = 19.7392088$ , the eigenvalue of the continuous operator.

Step 2 (relax on the $h = 1/8$ grid). There are $N = 49$ unknowns. The predicted rates are

\rho(E_{\mathrm{J}}) = \cos\tfrac{\pi}{8} = 0.92387953, \qquad \rho(E_{\mathrm{GS}}) = 0.85355339, \qquad \omega_{\mathrm{opt}} = \frac{2}{1+\sin(\pi/8)} = 1.44646269 ,

the last giving $\rho(E_{\omega_{\mathrm{opt}}}) = \omega_{\mathrm{opt}} - 1 = 0.44646269$ . Forming the three propagators explicitly and computing their eigenvalues confirms all three to ten figures, and scanning $\omega$ over $[1, 1.99]$ in steps of $0.001$ puts the numerical minimum of $\rho(E_{\omega})$ at $\omega = 1.447$ with value $0.4470$ , which is Young's formula to the resolution of the scan.

Step 3 (count the sweeps). Starting from $U^{(0)} = 0$ and iterating until the maximum-norm error $\|U^{(n)} - U^{\ast}\|_{\infty}$ has fallen by a factor $10^{-10}$ :

method	predicted $\ln(10^{-10})/\ln\rho$	sweeps actually taken
Jacobi	$290.8$	$291$
Gauss-Seidel	$145.4$	$146$
SOR at $\omega_{\mathrm{opt}}$	$28.6$	$35$

Jacobi and Gauss-Seidel land within one sweep of the asymptotic prediction, and Gauss-Seidel takes almost exactly half as many as Jacobi, as claim 5 says. SOR takes about a quarter more than the asymptotic count, the transient predicted in the remark above; it is still better than Gauss-Seidel by a factor of $4.2$ , against the asymptotic $2/(\pi h) = 5.09$ .

Step 4 (the cost, in operations). At $5$ operations per node per sweep, Gauss-Seidel used about $5 \times 49 \times 146 \approx 3.6\times10^{4}$ operations and SOR about $8 \times 49 \times 35 \approx 1.4\times10^{4}$ . Extrapolating with claim 5: on a grid with $h = 1/1000$ , so $N \approx 10^{6}$ , Gauss-Seidel would need $|\ln\varepsilon|/(\pi^{2}h^{2}) \approx 2.3\times10^{6}$ sweeps for the same tolerance and SOR about $3.7\times10^{3}$ . In two dimensions $h \asymp N^{-1/2}$ , so the sweep counts are $O(N)$ and $O(N^{1/2})$ and the total costs $O(N^{2})$ and $O(N^{3/2})$ ; the two-dimensional column of Stakgold and Holst's complexity table gives $O(N^{2}\ln N)$ and $O(N^{3/2}\ln N)$ , the extra logarithm coming from taking the tolerance down to discretisation-error level, $\varepsilon = O(h^{2}) = O(N^{-1})$ , as the grid is refined. The same table puts multilevel methods at $O(N\ln N)$ and nested multilevel and domain decomposition at $O(N)$ .

ApplicationOriginal

Every problem whose steady state is governed by a potential leads to the system (P): electrostatics with a prescribed charge density, steady heat conduction in a slab, Darcy flow in an aquifer, the stream function of a two-dimensional incompressible flow, and the gravitational potential of a mass distribution (Fundamental solution of Laplace's equation).

The place where these solvers are run hardest is the pressure projection of an incompressible flow solver. Enforcing $\nabla\cdot\mathbf{u} = 0$ at each time step requires the solution of a Poisson equation for the pressure on the whole grid, and because this happens every time step, the relaxation cost, not the physics, dominates the run. The numbers in the example are the reason: at $h = 1/1000$ , plain Gauss-Seidel would need millions of sweeps per time step, which is simply not available, while optimally relaxed SOR needs thousands and a multigrid solver tens. The progression Jacobi $\to$ Gauss-Seidel $\to$ SOR $\to$ multigrid is a progression in the exponent of $N$ , and in a production code only the last two are usable.

Relaxation also has a long pre-computer history: Southwell's hand-relaxation of structural and fluid problems in the 1930s and 1940s is literally the Gauss-Seidel sweep, performed with a pencil on a drawing of the grid, with the residual at each node written down and the largest one "relaxed" first. The over-relaxation trick, of overshooting the local correction, was discovered in that practice before it was analysed.

Finally, the smoothing property noted in the remark — that a sweep kills rough error and leaves smooth error — is used on its own, outside any solver, as an image-processing operation and as a pre-conditioner inside Krylov methods, where a few sweeps are applied not to solve the system but to make it easier for something else to solve.

Method

Spectral methods

Adapted

Add to PDF

Represents the unknown by its values on a Fourier or Chebyshev grid rather than by a difference stencil, so that differentiation becomes multiplication by $\mathrm{i}n$ or a dense differentiation matrix; proves that the interpolation error is at most twice the tail of the coefficient series, so the error falls like $N^{-k}$ for a $C^{k}$ function and geometrically for an analytic one; and identifies aliasing as the single mechanism behind both the error bound and the dealiasing rule for nonlinear terms.

Prerequisites: The discrete Fourier transform and the FFT, Chebyshev polynomials, Smoothness and decay of Fourier coefficients, Cauchy's theorem

Definition(Collocation: the two standard grids)AdaptedAdapted from Stakgold & Holst (2011), §2.3, Eqs. (2.3.30)–(2.3.32), pp. 140–142; Riley, Hobson & Bence (2006), §18.4, Eq. (18.55), p. 595; Riley, Hobson & Bence (2006), §27.4.3, Eq. (27.45), p. 1009; Arfken, Weber & Harris (2013), §18.4, Eqs. (18.104)–(18.106), p. 904

A collocation method (or pseudo-spectral method) stores the unknown as its values at $N$ chosen points, interprets those values as the values of a single global interpolant drawn from a fixed $N$ -dimensional space of smooth functions, and differentiates by differentiating the interpolant exactly. Two choices of space and grid cover almost all practice.

The Fourier grid (periodic problems). Take $N = 2M+1$ odd and the equally spaced nodes

x_j \;=\; \frac{2\pi j}{N},\qquad j = 0,1,\dots,N-1 ,

on the circle, and let $\mathcal{T}_M$ be the $N$ -dimensional space of trigonometric polynomials $\sum_{\lvert n\rvert\le M} c_n \mathrm{e}^{\mathrm{i}nx}$ . Stakgold and Holst prove that for any data $(f_j)$ there is exactly one $p\in\mathcal{T}_M$ with $p(x_j)=f_j$ for every $j$ , and that its coefficients are the discrete Fourier coefficients of the data. Writing $I_N f$ for the interpolant of the samples of a function $f$ , this says

I_N f(x) \;=\; \sum_{\lvert n\rvert\le M}\hat F(n)\,\mathrm{e}^{\mathrm{i}nx}, \qquad \hat F(n) \;=\; \frac{1}{N}\sum_{j=0}^{N-1} f(x_j)\,\mathrm{e}^{-2\pi\mathrm{i}nj/N},

with $\hat F$ the transform of The discrete Fourier transform and the FFT. The Fourier differentiation operator is

\mathcal{D}_N f \;=\; \bigl(I_N f\bigr)' , \qquad\text{so}\qquad \widehat{\mathcal{D}_N f}(n) = \mathrm{i}n\,\hat F(n)\ \ (\lvert n\rvert \le M) .

In words: transform, multiply by $\mathrm{i}n$ , transform back. An odd $N$ is taken so that the retained frequencies $-M,\dots,M$ are symmetric and no lone mode sits at $n=N/2$ ; this is Stakgold and Holst's reason for writing the node count as $2N+1$ . With $N$ even one must decide what to do with the Nyquist mode $n=N/2$ , and the usual choice is to set its derivative to zero, which keeps $\mathcal{D}_N$ real and skew-symmetric.

The Chebyshev grid (non-periodic problems on an interval). On $[-1,1]$ replace trigonometric polynomials by algebraic ones, $\mathcal{P}_N = \operatorname{span}\{1,x,\dots,x^{N}\}$ , and take the Chebyshev–Gauss–Lobatto (CGL) nodes

x_j \;=\; \cos\frac{\pi j}{N},\qquad j=0,1,\dots,N ,

which include both endpoints $x_0=1$ and $x_N=-1$ . These are the images under $x=\cos\theta$ of the equally spaced angles $\theta_j = \pi j/N$ . Because $T_k(\cos\theta)=\cos k\theta$ (Chebyshev polynomials), a polynomial $\sum_{k\le N} a_k T_k(x)$ is a cosine polynomial in $\theta$ evaluated at equally spaced angles: Chebyshev collocation on $[-1,1]$ is Fourier collocation on the circle, restricted to even data. The Chebyshev differentiation matrix $D_N\in\mathbb{R}^{(N+1)\times(N+1)}$ is defined by differentiating the Lagrange interpolant,

(D_N)_{jk} \;=\; \ell_k'(x_j), \qquad \ell_k(x)=\prod_{m\neq k}\frac{x-x_m}{x_k-x_m} ,

so that $(D_N u)_j$ is the derivative at $x_j$ of the unique polynomial of degree $\le N$ through the data. Unlike its Fourier counterpart, $D_N$ is a full matrix.

A warning about "Chebyshev points". Two different node sets carry that name. The CGL nodes above are the extrema of $T_N$ together with the endpoints. The Chebyshev–Gauss nodes are the zeros of $T_n$ , $x_i=\cos\bigl((i-\tfrac12)\pi/n\bigr)$ , $i=1,\dots,n$ ; Riley, Hobson and Bence give them, with the equal weights $w_i=\pi/n$ , as the points of the Gauss–Chebyshev quadrature rule for the weight $(1-x^{2})^{-1/2}$ (see Numerical quadrature). Both sets cluster near the endpoints with spacing $O(N^{-2})$ , and both are used; only the Lobatto set contains the endpoints, which is what makes it the convenient one for boundary conditions.

Method(Aliasing and spectral accuracy)Original

Let $f$ be $2\pi$ -periodic and continuous with an absolutely convergent Fourier series, $\sum_{n}\lvert\hat f(n)\rvert<\infty$ , and let $N=2M+1$ .

1. Aliasing. The interpolation coefficients are the exact coefficients with every congruent frequency folded in:

\hat F(n) \;=\; \sum_{m \,\equiv\, n \ (\mathrm{mod}\ N)} \hat f(m), \qquad \lvert n\rvert\le M .

2. The interpolation error is at most twice the tail.

\bigl\lVert f - I_N f\bigr\rVert_{\infty} \;\le\; 2\!\!\sum_{\lvert m\rvert > M}\!\!\bigl\lvert \hat f(m)\bigr\rvert .

3. Smoothness buys algebraic order. If $f\in C^{k}$ with $f^{(k)}$ of bounded variation, so that $\lvert\hat f(m)\rvert\le C\lvert m\rvert^{-k-1}$ , then

\bigl\lVert f - I_N f\bigr\rVert_{\infty} \;\le\; \frac{4C}{k}\,M^{-k} \;=\; O\!\left(N^{-k}\right) .

The order of the method is not fixed in advance by a stencil: it is whatever the smoothness of the solution allows. This is what is meant by spectral accuracy.

4. Analyticity buys geometric convergence. If $f$ extends holomorphically to the strip $\lvert\operatorname{Im}z\rvert<b$ and satisfies $\lvert f\rvert\le A$ on $\lvert \operatorname{Im}z\rvert\le a$ for some $a<b$ , then $\lvert\hat f(m)\rvert\le A\mathrm{e}^{-a\lvert m\rvert}$ and

\bigl\lVert f - I_N f\bigr\rVert_{\infty} \;\le\; \frac{4A\,\mathrm{e}^{-aM}}{1-\mathrm{e}^{-a}} .

5. The Chebyshev grid aliases in the same way. On the CGL nodes $x_j=\cos(\pi j/N)$ ,

T_{2N-k}(x_j) = T_{2N+k}(x_j) = T_k(x_j)\qquad\text{for all } j \text{ and all } k\ge 0 ,

so the Chebyshev interpolant of $f$ on those nodes has coefficients $\tilde a_k=\sum a_m$ summed over the $m$ that fold onto $k$ , and parts 2–4 hold verbatim for the Chebyshev coefficients $a_k$ of $f$ in place of $\hat f(n)$ , with $\lVert\cdot\rVert_\infty$ taken on $[-1,1]$ .

Idea of the proofOriginal

A finite grid cannot tell two frequencies apart when they differ by a multiple of $N$ : the exponentials $\mathrm{e}^{\mathrm{i}nx}$ and $\mathrm{e}^{\mathrm{i}(n+N)x}$ agree at every node. So the interpolant carries each true coefficient plus the coefficients of all its impostors. That is part 1, and it is already proved in The discrete Fourier transform and the FFT.

Part 2 follows by splitting the error into two pieces that are each controlled by the same tail: the coefficients that were thrown away because their frequency was too high, and the pollution that those discarded coefficients left behind in the retained ones. No coefficient appears in both piles, so the two sums together are at most twice $\sum_{\lvert m\rvert>M}\lvert\hat f(m)\rvert$ .

Parts 3 and 4 are then nothing but the decay rate of Fourier coefficients, already established in Smoothness and decay of Fourier coefficients for the smooth case and obtained here by moving the contour of the coefficient integral into the complex strip for the analytic case. The whole content of spectral accuracy is that the error of the method is the tail of a series, not a fixed power of the mesh. Part 5 is part 1 read through the substitution $x=\cos\theta$ , which carries the CGL grid to an equally spaced grid in $\theta$ and the Chebyshev polynomials to cosines.

Proof(Parts 1–5)Original

1. This is the aliasing formula of The discrete Fourier transform and the FFT: substituting $f(x_j)=\sum_m\hat f(m)\mathrm{e}^{\mathrm{i}mx_j}$ into the definition of $\hat F(n)$ and exchanging the two sums, which the hypothesis $\sum_m\lvert\hat f(m)\rvert<\infty$ permits, leaves $\sum_m \hat f(m)\,\delta$ , where $\delta$ is $1$ when $m\equiv n$ modulo $N$ and $0$ otherwise, by the discrete orthogonality relation.

2. Absolute convergence of the coefficient series makes the Fourier expansion of $f$ uniformly convergent, so for every $x$ we may write

f(x)-I_Nf(x) =\underbrace{\sum_{\lvert m\rvert>M}\hat f(m)\,\mathrm{e}^{\mathrm{i}mx}}_{\text{(T) truncation}} \;-\;\underbrace{\sum_{\lvert n\rvert\le M}\bigl(\hat F(n)-\hat f(n)\bigr)\mathrm{e}^{\mathrm{i}nx}}_{\text{(A) aliasing}} .

Bounding each exponential by $1$ gives $\lvert (\mathrm{T})\rvert\le\sum_{\lvert m\rvert>M}\lvert\hat f(m)\rvert$ , which is the first copy of the tail.

For (A), part 1 gives $\hat F(n)-\hat f(n)=\sum \hat f(m)$ , the sum taken over those $m\equiv n$ modulo $N$ with $m\ne n$ . Every such $m$ satisfies $\lvert m\rvert\ge N-M=M+1$ , because the nearest congruent frequency to an $n$ with $\lvert n\rvert\le M$ is $n\pm N$ and $\lvert n\pm N\rvert\ge N-M$ . Moreover, as $n$ runs over $-M,\dots,M$ these index sets are disjoint: two frequencies congruent modulo $N$ to different $n$ in a window of $N$ consecutive integers are different. Hence

\sum_{\lvert n\rvert\le M}\bigl\lvert \hat F(n)-\hat f(n)\bigr\rvert \;\le\; \sum_{\lvert m\rvert>M}\bigl\lvert\hat f(m)\bigr\rvert ,

the second copy of the tail. Adding the two bounds gives part 2. $\square$ (for part 2)

3. With $\lvert\hat f(m)\rvert\le C\lvert m\rvert^{-k-1}$ the tail is bounded by comparison with an integral:

\sum_{\lvert m\rvert>M}\lvert\hat f(m)\rvert \;\le\; 2C\sum_{m=M+1}^{\infty}m^{-k-1} \;\le\; 2C\int_{M}^{\infty}t^{-k-1}\,\mathrm{d}t \;=\;\frac{2C}{k}\,M^{-k} ,

the middle step because $m^{-k-1}\le\int_{m-1}^{m}t^{-k-1}\mathrm{d}t$ for $m\ge M+1$ . Doubling, by part 2, gives the stated constant. The hypothesis that $f^{(k)}$ has bounded variation is what produces the decay rate $\lvert m\rvert^{-k-1}$ ; this is the content of Smoothness and decay of Fourier coefficients, where repeated integration by parts transfers each derivative onto a factor $\mathrm{i}m$ .

4. Fix $m>0$ and consider the rectangle with corners $0,\ 2\pi,\ 2\pi-\mathrm{i}a,\ -\mathrm{i}a$ . The integrand $f(z)\mathrm{e}^{-\mathrm{i}mz}$ is holomorphic inside and on it, since $a<b$ , so by Cauchy's theorem the contour integral around it vanishes. The two vertical sides cancel because $f(z)\mathrm{e}^{-\mathrm{i}mz}$ is $2\pi$ -periodic in $z$ — here periodicity is used a second time, now in the complex plane. Therefore the integral over $[0,2\pi]$ equals the integral over the shifted segment $[-\mathrm{i}a,2\pi-\mathrm{i}a]$ , and

\hat f(m)=\frac{1}{2\pi}\int_{0}^{2\pi}f(t-\mathrm{i}a)\,\mathrm{e}^{-\mathrm{i}m(t-\mathrm{i}a)}\,\mathrm{d}t , \qquad \bigl\lvert \mathrm{e}^{-\mathrm{i}m(t-\mathrm{i}a)}\bigr\rvert=\mathrm{e}^{-ma} ,

so $\lvert\hat f(m)\rvert\le A\mathrm{e}^{-ma}$ . For $m<0$ shift the contour the other way, to $+\mathrm{i}a$ , and the same estimate appears with $\lvert m\rvert$ . Summing the geometric series,

\sum_{\lvert m\rvert>M}\lvert\hat f(m)\rvert\le 2A\sum_{m=M+1}^{\infty}\mathrm{e}^{-ma} =\frac{2A\,\mathrm{e}^{-a(M+1)}}{1-\mathrm{e}^{-a}}\le \frac{2A\,\mathrm{e}^{-aM}}{1-\mathrm{e}^{-a}} ,

and part 2 doubles it.

5. At the node $x_j=\cos(\pi j/N)$ we have $T_k(x_j)=\cos(k\pi j/N)$ by the defining identity $T_k(\cos\theta)=\cos k\theta$ . Hence

T_{2N\pm k}(x_j)=\cos\!\Bigl(\frac{(2N\pm k)\pi j}{N}\Bigr)=\cos\!\Bigl(2\pi j\pm\frac{k\pi j}{N}\Bigr)=\cos\frac{k\pi j}{N}=T_k(x_j),

using $\cos(2\pi j+\varphi)=\cos\varphi$ for integer $j$ and the evenness of the cosine for the minus sign. So degrees $k$ and $2N\pm k$ are indistinguishable on the grid, exactly as frequencies $n$ and $n\pm N$ are on the Fourier grid. For the error bound, set $g(\theta)=f(\cos\theta)$ ; $g$ is $2\pi$ -periodic and even, its Fourier cosine coefficients are the Chebyshev coefficients $a_k$ of $f$ , the CGL nodes are the points $\theta_j=\pi j/N$ , and $\sup_{[-1,1]}\lvert f-I_Nf\rvert=\sup_{\theta}\lvert g - I g\rvert$ because $\theta\mapsto\cos\theta$ maps $[0,\pi]$ onto $[-1,1]$ . Parts 2–4 applied to $g$ therefore give the claim. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Periodicity is not a technicality. Part 2 needs the coefficients of a periodic $f$ to be summable; part 4 uses periodicity again, to cancel the vertical sides of the contour. Apply a Fourier method to a function that is smooth on $[0,2\pi]$ but has $f(0)\ne f(2\pi)$ and the periodic extension has a jump, so $\lvert \hat f(m)\rvert$ decays only like $1/\lvert m\rvert$ , the tail is $O(1)$ , and the interpolant oscillates near the ends with an overshoot that does not shrink: this is exactly The Gibbs phenomenon, and no amount of resolution removes it. The Chebyshev grid exists precisely to deal with this case: $g(\theta)=f(\cos\theta)$ is automatically periodic and even, whatever $f$ does at the endpoints.

Smoothness sets the order, not the grid. In part 3 the order $k$ comes from the function and not from the method. A solution with a shock or a corner gives first-order accuracy from a spectral method just as from a difference scheme, and worse: the error is spread over the whole domain rather than localised at the singularity, because the basis functions are global. Compare Finite-difference approximations, whose order is fixed by the stencil and is therefore insensitive, for better and for worse, to how smooth the solution is.

Clustering of the nodes is essential for Chebyshev. Taking $N+1$ equally spaced nodes on $[-1,1]$ and interpolating by a polynomial of degree $N$ is not a spectral method at all: the nodal polynomial $\omega(x)=\prod_j (x-x_j)$ , which carries the interpolation error (Polynomial and spline interpolation), is then exponentially larger near the ends than in the middle, and the interpolant of a perfectly analytic function can diverge as $N$ grows. The CGL nodes are, up to the endpoints, the extrema of $T_N$ , and Kreyszig's Theorem 6.4-3 says that $2^{1-N}T_N$ is the monic polynomial of degree $N$ with the smallest maximum modulus on $[-1,1]$ , namely $2^{1-N}$ . No node set can make the nodal polynomial smaller than that; the Chebyshev nodes very nearly achieve it.

The summability hypothesis. That $\sum_m\lvert\hat f(m)\rvert$ is finite is what lets us exchange the sums in part 1 and bound the exponentials termwise in part 2. It holds under the hypotheses of parts 3 (for $k\ge1$ ) and 4. Without it the interpolant still exists — interpolation is a finite linear-algebra problem — but the bound in part 2 carries no information.

Sharpness. Part 2 is sharp to within a small factor, not merely an order of magnitude: the example below finds the predicted bound exceeding the measured error by $7\%$ at $N=9$ and by $4\%$ at $N=17$ .

Method(Running a spectral solver)Original

Fourier collocation (periodic domain).

Choose the grid. $N$ points $x_j=2\pi j/N$ on $[0,2\pi)$ , rescaled to the physical period $L$ by $x\mapsto 2\pi x/L$ , which multiplies every wavenumber by $2\pi/L$ . Use $N$ a power of two if the FFT library wants one, and remember the Nyquist convention of the definition; use $N$ odd if you prefer the symmetric mode set.
Transform. Compute $\hat F(n)$ for $\lvert n\rvert\le M$ with the FFT, at a cost of $O(N\log N)$ (The discrete Fourier transform and the FFT).
Apply the operator in frequency space. Each $\partial_x$ is a multiplication by $\mathrm{i}n$ , each $\partial_x^2$ by $-n^{2}$ . A constant-coefficient linear problem $L(\partial_x)u=f$ is therefore solved by one division per mode: $\hat U(n)=\hat F(n)/L(\mathrm{i}n)$ , with the solvability condition $L(0)\neq0$ or else an orthogonality condition on $\hat F(0)$ .
Evaluate nonlinear terms in physical space. To form $u\partial_x u$ : multiply by $\mathrm{i}n$ , transform back, multiply pointwise by the grid values of $u$ , transform forward again. Two forward and two inverse transforms per step, so $O(N\log N)$ — against $O(N^{2})$ for the same operation done as a convolution in frequency space.
Dealias. A quadratic nonlinearity of modes with $\lvert n\rvert\le K$ produces modes up to $2K$ , and the grid folds a mode $p$ with $\lvert p\rvert>M$ back onto $p\mp N$ . The folded modes first land inside the retained band when $N-2K\le K$ . Keeping $3K<N$ , that is zeroing all modes with $\lvert n\rvert>N/3$ before and after each product, removes the contamination; this is the two-thirds rule.
Step in time with any integrator from Integrators for ODEs, bearing in mind that the diagonal entries $-n^{2}$ of the second-derivative operator reach $-M^{2}\approx -N^{2}/4$ , so an explicit scheme for a diffusion problem is restricted to $\Delta t = O(N^{-2})$ and the problem is stiff in the sense of Stability and stiff equations. Treating the linear part implicitly costs nothing, because it is diagonal.

Chebyshev collocation (non-periodic interval).

Map the interval. Send $[\alpha,\beta]$ to $[-1,1]$ by $x\mapsto (2x-\alpha-\beta)/(\beta-\alpha)$ ; every derivative acquires a factor $2/(\beta-\alpha)$ .
Build $D_N$ once, from $(D_N)_{jk}=\ell_k'(x_j)$ at the CGL nodes, and form $D_N^{2}$ for a second-order problem. Cost $O(N^{2})$ to store, $O(N^{3})$ to form the square or factorise it.
Impose boundary conditions by deletion. For $u(\pm1)=0$ , strike out the first and last rows and columns of $D_N^{2}$ ; the remaining $(N-1)\times(N-1)$ block is invertible and acts on the interior values. For a Neumann condition, keep the row and replace it by the corresponding row of $D_N$ .
Solve the dense system by LU, at a cost of $\tfrac23 N^{3}+O(N^{2})$ operations once and $O(N^{2})$ per right-hand side thereafter (Solving linear systems numerically).
If $N$ is large, transform instead. The map between nodal values and Chebyshev coefficients is a discrete cosine transform, computable in $O(N\log N)$ , and the coefficient recurrence for differentiation is $O(N)$ ; this replaces the dense matrix when $N$ runs into the hundreds.
Watch the conditioning. The CGL nodes are spaced $O(N^{-2})$ near the endpoints, so $\lVert D_N^{2}\rVert$ grows like $N^{4}$ and the condition number of the interior block grows at the same rate (Matrix norms and condition numbers). Spectral accuracy is therefore eventually limited by rounding: in the example below the error stops improving at about $5\times10^{-15}$ .

Example(A periodic Helmholtz problem, and a Chebyshev two-point problem)Original

(a) How large is the interpolation error, really? Take $f(x)=\mathrm{e}^{\sin x}$ , which is entire, hence analytic in every strip. Its Fourier coefficients are Bessel functions. Put $t=\mathrm{e}^{\mathrm{i}x}$ in the generating function $\mathrm{e}^{(z/2)(t-1/t)}=\sum_n J_n(z)t^{n}$ of Bessel functions of the first kind: the exponent becomes $\mathrm{i}z\sin x$ , so the choice $z=-\mathrm{i}$ gives $\mathrm{e}^{\sin x}=\sum_n J_n(-\mathrm{i})\,\mathrm{e}^{\mathrm{i}nx}$ . Since $J_n(-z)=(-1)^nJ_n(z)$ and $I_n(x)=\mathrm{i}^{-n}J_n(\mathrm{i}x)$ (Modified Bessel functions), this is $\hat f(n)=(-\mathrm{i})^{n}I_n(1)$ , so $\lvert\hat f(n)\rvert=I_n(1)$ . Numerically

I_0(1)=1.266066,\quad I_1(1)=0.565159,\quad I_2(1)=0.135748,\quad I_3(1)=2.216842\times10^{-2},

I_4(1)=2.737120\times10^{-3},\quad I_5(1)=2.714632\times10^{-4},\quad I_6(1)=2.248866\times10^{-5} .

For $N=9$ we have $M=4$ , so part 2 predicts

\lVert f-I_9f\rVert_\infty\;\le\;4\sum_{n\ge5}I_n(1) = 4\bigl(2.714632\times10^{-4}+2.248866\times10^{-5}+1.599218\times10^{-6}+\cdots\bigr) = 1.18263\times10^{-3},

the factor $4$ being the $2$ of part 2 times the $2$ from $\pm n$ . Evaluating the interpolant on a fine grid of $4001$ points gives a measured maximum error of $1.109873\times10^{-3}$ : the bound is high by $7\%$ . At $N=17$ the same computation predicts $2.32268\times10^{-8}$ and measures $2.223491\times10^{-8}$ , high by $4\%$ . At $N=33$ the bound is $8.9\times10^{-20}$ and the measurement is $1.8\times10^{-15}$ , which is where double-precision rounding takes over: the error of the method has long since fallen below the error of the arithmetic.

(b) Solving $-u''+u=f$ on the circle. Let $u(x)=\mathrm{e}^{\sin x}$ again. Then $u'=\cos x\,\mathrm{e}^{\sin x}$ and $u''=(\cos^{2}x-\sin x)\mathrm{e}^{\sin x}$ , so the right-hand side is

f(x)=-u''+u=\bigl(1-\cos^{2}x+\sin x\bigr)\mathrm{e}^{\sin x}=\bigl(\sin^{2}x+\sin x\bigr)\mathrm{e}^{\sin x} .

Step 3 of the method is one line: $\hat U(n)=\hat F(n)/(n^{2}+1)$ . Carrying it out and comparing with $u$ at the nodes:

$N$	spectral, $\max_j\lvert U_j-u(x_j)\rvert$	three-point differences, same $N$
$9$	$3.984146\times10^{-4}$	$7.788135\times10^{-2}$
$17$	$3.176116\times10^{-9}$	$2.160161\times10^{-2}$
$33$	$4.996004\times10^{-16}$	$5.716368\times10^{-3}$
$65$	$1.054712\times10^{-15}$	$1.472297\times10^{-3}$

The second column falls by a factor of about $4$ each time $N$ doubles — second order, as Finite-difference approximations predicts for the stencil $(U_{j+1}-2U_j+U_{j-1})/h^{2}$ . The first column falls by five orders of magnitude from $N=9$ to $N=17$ and by seven more from $N=17$ to $N=33$ , where it meets machine precision and stops. To reach $5\times10^{-16}$ with the difference scheme at second order would need $N\approx 33\sqrt{5.716\times10^{-3}/(5\times10^{-16})}\approx 1.1\times10^{8}$ points.

(c) Chebyshev collocation for $u''=f$ on $[-1,1]$ , $u(\pm1)=0$ . Take $u(x)=(1-x^{2})\mathrm{e}^{x}$ , which vanishes at both ends, and compute

u'=(1-2x-x^{2})\mathrm{e}^{x},\qquad u''=f(x)=-(x^{2}+4x+1)\,\mathrm{e}^{x} .

With $N=4$ the CGL nodes are $x_j=\cos(\pi j/4)$ , that is

x_0=1,\quad x_1=\tfrac{\sqrt2}{2},\quad x_2=0,\quad x_3=-\tfrac{\sqrt2}{2},\quad x_4=-1 .

Differentiating the five Lagrange basis polynomials and evaluating at the nodes gives

D_4=\begin{pmatrix} \tfrac{11}{2} & -4-2\sqrt2 & 2 & -4+2\sqrt2 & \tfrac12\\[2pt] 1+\tfrac{\sqrt2}{2} & -\tfrac{\sqrt2}{2} & -\sqrt2 & \tfrac{\sqrt2}{2} & -1+\tfrac{\sqrt2}{2}\\[2pt] -\tfrac12 & \sqrt2 & 0 & -\sqrt2 & \tfrac12\\[2pt] 1-\tfrac{\sqrt2}{2} & -\tfrac{\sqrt2}{2} & \sqrt2 & \tfrac{\sqrt2}{2} & -1-\tfrac{\sqrt2}{2}\\[2pt] -\tfrac12 & 4-2\sqrt2 & -2 & 4+2\sqrt2 & -\tfrac{11}{2} \end{pmatrix} .

(Two checks: every row sums to zero, because $D_4$ differentiates the constant function to zero; and the centre row, $\bigl(-\tfrac12,\sqrt2,0,-\sqrt2,\tfrac12\bigr)$ , is antisymmetric, as it must be because the node set is symmetric about the origin.) Squaring and deleting the first and last rows and columns leaves the interior operator — remarkably, an integer matrix:

\bigl(D_4^{2}\bigr)_{\text{int}}=\begin{pmatrix}-14 & 6 & -2\\ 4 & -6 & 4\\ -2 & 6 & -14\end{pmatrix}, \qquad F=\begin{pmatrix}f(\sqrt2/2)\\ f(0)\\ f(-\sqrt2/2)\end{pmatrix} =\begin{pmatrix}-8.778548\\ -1\\ \phantom{-}0.655006\end{pmatrix} .

Solving the $3\times3$ system gives, against the exact values,

node	$U_j$ (collocation)	$u(x_j)$ exact	error
$\sqrt2/2$	$1.025786$	$1.014057$	$+1.17\times10^{-2}$
$0$	$1.010295$	$1$	$+1.03\times10^{-2}$
$-\sqrt2/2$	$0.239657$	$0.246534$	$-6.88\times10^{-3}$

Five points, three unknowns, one percent. Repeating with larger $N$ :

$N$	$4$	$8$	$12$	$16$	$20$
max error	$1.17\times10^{-2}$	$4.85\times10^{-7}$	$2.41\times10^{-12}$	$5.22\times10^{-15}$	$1.33\times10^{-15}$

Roughly five digits gained for every four points added, until the $N^{4}$ growth of $\operatorname{cond}(D_N^{2})$ stops the descent at about $5\times10^{-15}$ — the limit anticipated in step 6 of the method.

(d) Aliasing on the Chebyshev grid, by hand. With $N=4$ , part 5 says $T_{2N-1}=T_7$ should agree with $T_1(x)=x$ at all five nodes. At $x=\sqrt2/2$ , using $T_7(x)=64x^{7}-112x^{5}+56x^{3}-7x$ ,

64\Bigl(\tfrac{\sqrt2}{2}\Bigr)^{7}-112\Bigl(\tfrac{\sqrt2}{2}\Bigr)^{5}+56\Bigl(\tfrac{\sqrt2}{2}\Bigr)^{3}-7\cdot\tfrac{\sqrt2}{2} = \tfrac{\sqrt2}{2}\bigl(64\cdot\tfrac18-112\cdot\tfrac14+56\cdot\tfrac12-7\bigr) =\tfrac{\sqrt2}{2}(8-28+28-7)=\tfrac{\sqrt2}{2} ,

which is $x$ itself. A degree- $7$ oscillation is invisible on a five-point Chebyshev grid, exactly as a frequency- $7$ oscillation is invisible on a five-point Fourier grid.

Application(Where spectral methods are used)Original

Turbulence and direct numerical simulation. The reference simulations of homogeneous isotropic turbulence are pseudo-spectral: the Navier–Stokes equations in a periodic box, with the viscous term diagonal in Fourier space, the nonlinear term evaluated on the grid, and the two-thirds rule applied at every step. The appeal is that the energy spectrum — the quantity the physics is about — is a primary variable of the method rather than something reconstructed afterwards, and that the dissipation range is resolved without the numerical viscosity that a low-order scheme adds.

Weather and climate models. Global atmospheric models have long used spherical harmonics as the horizontal basis: the Laplacian is diagonal, which makes semi-implicit time stepping cheap, and there is no pole problem of the kind a latitude–longitude difference grid suffers. The transform between grid-point and spectral representation is the dominant cost, and it is a Fourier transform in longitude combined with a Legendre transform in latitude (Spherical harmonics).

Quantum dynamics. The split-step Fourier method evolves the Schrödinger equation (The Schrodinger equation) by alternating an exact kinetic step, diagonal in momentum space because $\hat{p}=-\mathrm{i}\hbar\partial_x$ is multiplication by $\hbar k$ , with an exact potential step, diagonal in position space. Both steps are unitary, so the norm is conserved to rounding error, and each costs one pair of transforms. The same structure underlies plane-wave electronic-structure codes, where the basis is the set of Fourier modes commensurate with a crystal lattice and the cutoff is quoted as an energy $\hbar^{2}k_{\max}^{2}/2m$ .

Why not always? Complex geometry. A Fourier basis needs a periodic box and a Chebyshev basis needs a tensor-product interval, so a domain with a curved or re-entrant boundary needs either a mapping, a decomposition into spectral elements, or the finite-element method of The finite element method, which gives up spectral accuracy and gains complete geometric freedom. Shocks are the other obstacle: for the conservation laws of Conservation laws and shocks the solution leaves the smoothness class in finite time, the tail of the coefficient series stops decaying, and the global basis spreads the resulting oscillations across the whole domain.

Method · advanced

The finite element method

Adapted

Add to PDF

Replaces the weak formulation of an elliptic boundary-value problem by the same equation tested only against piecewise-linear hat functions on a triangulation, which turns it into a sparse symmetric positive-definite linear system; proves Galerkin orthogonality and Cea's quasi-optimality bound, states the interpolation estimate that converts quasi-optimality into an $O(h)$ energy-norm error, and shows on a worked mesh why the method is the one that copes with awkward geometry.

Prerequisites: Weak derivatives and Sobolev spaces, The Lax-Milgram theorem and the Dirichlet problem, Quadratic forms and the Rayleigh quotient

Definition(Triangulation, hat functions and the Galerkin system)AdaptedAdapted from Strauss (2008), §8.5, Figures 1–2 and Eqs. (1)–(5), pp. 222–224; Stakgold & Holst (2011), §10.5, Assumptions 10.5.1–10.5.3 and Definition 10.5.1, pp. 737–739

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded polygon. The model problem is the one Strauss takes: Poisson's equation on $\Omega$ with the value of the unknown prescribed to be zero all along the boundary,

-\Delta u = f \ \text{ in }\Omega,\qquad u=0 \ \text{ on }\partial\Omega ,

in the weak form of Weak derivatives and Sobolev spaces: find $u\in H^{1}_{0}(\Omega)$ with

a(u,v)\;:=\;\int_\Omega \nabla u\cdot\nabla v \,\mathrm{d}x \;=\;\int_\Omega f v\,\mathrm{d}x\;=:\;\ell(v) \qquad\text{for every } v\in H^{1}_{0}(\Omega) .

Triangulation. A triangulation $\mathcal{T}_h$ of $\Omega$ is a finite family of closed triangles whose union is $\overline\Omega$ , with disjoint interiors, and such that the intersection of two distinct triangles is empty, a common vertex, or a common whole edge. Stakgold and Holst call the last condition conformity and state it for elements of any shape in any dimension; it is what forbids a vertex of one triangle from landing in the middle of a neighbour's edge. Write $h_K$ for the diameter of the smallest ball containing $K$ and $\rho_K$ for the diameter of the largest ball inside $K$ , put $h=\max_K h_K$ , and call the family shape regular with constant $\sigma$ if

\frac{h_K}{\rho_K}\le\sigma\qquad\text{for every } K\in\mathcal{T}_h ,

which says no triangle is allowed to become arbitrarily thin as the mesh is refined.

The space and its basis. Let $V_1,\dots,V_N$ be the vertices of the triangulation that lie in the interior of $\Omega$ , and let

V_h=\bigl\{\,v\in C(\overline\Omega)\ :\ v\ \text{is affine on each } K\in\mathcal{T}_h,\ v=0 \text{ on }\partial\Omega\,\bigr\} .

For each interior vertex $V_i$ let $v_i\in V_h$ be the unique element with

v_i(V_j)=\delta_{ij} \qquad (j=1,\dots,N)

and $v_i=0$ at the boundary vertices. Each $v_i$ is determined on each triangle by its three vertex values, so it exists and is unique; its graph is a pyramid of height $1$ with apex over $V_i$ , supported on the triangles that touch $V_i$ . These are the hat functions (Strauss's trial functions). They are a basis of $V_h$ , and the expansion of any $v\in V_h$ reads off its nodal values: $v=\sum_j v(V_j)\,v_j$ .

That $V_h\subset H^{1}_0(\Omega)$ is a small theorem, not a definition: a function that is a polynomial on each element and merely continuous across the edges has a weak gradient, because in the integration by parts the jump terms on each interior edge appear twice with opposite normals and cancel. Stakgold and Holst prove exactly this.

The discrete problem. The Galerkin approximation is: find $u_h\in V_h$ with

a(u_h,v_h)=\ell(v_h)\qquad\text{for every } v_h\in V_h .

Writing $u_h=\sum_i U_i v_i$ and testing against each basis function in turn turns this into the $N\times N$ linear system

\sum_{i=1}^{N} A_{ji}\,U_i = F_j \quad (j=1,\dots,N), \qquad A_{ji}=\int_\Omega \nabla v_i\cdot\nabla v_j\,\mathrm{d}x, \qquad F_j=\int_\Omega f\,v_j\,\mathrm{d}x .

$A$ is the stiffness matrix and $F$ the load vector; when a time derivative is present the matrix $\int_\Omega v_iv_j$ also appears and is called the mass matrix. Two properties are immediate. $A$ is symmetric, and it is positive definite because $\xi^{\mathsf T}A\xi=\int_\Omega\lvert\nabla(\sum_i\xi_iv_i)\rvert^{2}\ge0$ with equality only if $\sum\xi_iv_i$ is constant, hence zero by the boundary condition (Quadratic forms and the Rayleigh quotient). And $A$ is sparse: $A_{ji}=0$ unless $V_i$ and $V_j$ are joined by an edge, because otherwise the supports of $v_i$ and $v_j$ meet in a set of measure zero. Finally $U_k=u_h(V_k)$ , so the unknowns are the nodal values themselves.

A warning about names. Strauss writes the matrix with entries $\int\nabla v_i\cdot\nabla v_j$ as $M$ and calls it the stiffness matrix, and writes the matrix with entries $\int v_iv_j$ as $K$ and calls it the mass matrix. The names are his and are standard; the letters are the reverse of the usual engineering ones, in which $K$ is the stiffness and $M$ the mass. This unit keeps Strauss's names and uses neutral letters $A$ for the stiffness matrix and $\ell$ , $F$ for the load.

Method(Galerkin orthogonality, quasi-optimality, and the rate)AdaptedAdapted from Stakgold & Holst (2011), §10.2, Eqs. (10.2.16) and (10.2.19)–(10.2.22), Theorem 10.2.9 and Lemma 10.2.3, pp. 677–679; Stakgold & Holst (2011), §10.5, Theorem 10.5.2, pp. 744–745; Evans (2010), §6.2.1, Theorem 1, pp. 297–299; Evans (2010), §6.3.2, Theorem 4, p. 315

Let $H$ be a real Hilbert space, let $a:H\times H\to\mathbb{R}$ be bilinear with

\lvert a(u,v)\rvert\le M\lVert u\rVert\,\lVert v\rVert\quad (M>0), \qquad a(u,u)\ \ge\ m\lVert u\rVert^{2}\quad (m>0),

and let $\ell$ be a bounded linear functional on $H$ . Let $V_h\subset H$ be a closed subspace.

1. Both problems are uniquely solvable. There is exactly one $u\in H$ with $a(u,v)=\ell(v)$ for all $v\in H$ , and exactly one $u_h\in V_h$ with $a(u_h,v_h)=\ell(v_h)$ for all $v_h\in V_h$ . This is the Lax–Milgram theorem of The Lax-Milgram theorem and the Dirichlet problem applied twice: the hypotheses on $a$ are inherited by the subspace, which is why no separate condition on $V_h$ is needed. (For a non-symmetric Petrov–Galerkin method, where the trial and test spaces differ, the inherited condition fails and must be imposed separately; Stakgold and Holst call it the Ladyzhenskaya–Babuska–Brezzi condition.)

2. Galerkin orthogonality.

a(u-u_h,\,v_h)=0\qquad\text{for every } v_h\in V_h .

3. Quasi-optimality (Cea).

\lVert u-u_h\rVert \;\le\; \frac{M}{m}\,\inf_{v_h\in V_h}\lVert u-v_h\rVert .

The Galerkin solution is, up to the factor $M/m$ , the best approximation to $u$ available in $V_h$ . If $a$ is in addition symmetric, then $a$ is an inner product, $u_h$ is the $a$ -orthogonal projection of $u$ onto $V_h$ , and in the energy norm $\lVert v\rVert_a=a(v,v)^{1/2}$ the constant is exactly $1$ : $\lVert u-u_h\rVert_a=\inf_{v_h}\lVert u-v_h\rVert_a$ .

4. The interpolation estimate. Let $\mathcal{T}_h$ be shape regular with constant $\sigma$ and let $I_h w$ denote the element of $V_h$ that agrees with $w$ at every vertex. There is a constant $C$ depending only on $\sigma$ with

\lVert w-I_hw\rVert_{L^{2}(\Omega)}\le C h^{2}\lvert w\rvert_{H^{2}(\Omega)}, \qquad \lvert w-I_hw\rvert_{H^{1}(\Omega)}\le C h\,\lvert w\rvert_{H^{2}(\Omega)},

for every $w\in H^{2}(\Omega)\cap H^{1}_0(\Omega)$ . This is the case $k=1$ , $p=q=2$ , $m=0$ and $m=1$ of the general estimate of Stakgold and Holst, which for an affine-equivalent family reads

\lvert w-I_Kw\rvert_{m,p,K}\ \le\ C\bigl(\mathrm{meas}(\Omega)\bigr)^{(1/q)-(1/p)}\,\frac{h_K^{\,k+1}}{\rho_K^{\,m}}\,\lvert w\rvert_{k+1,p,K}

on each element $K$ ; at $p=q$ the measure factor is $1$ , and shape regularity converts $h_K^{k+1}/\rho_K^{m}$ into $\sigma^{m}h_K^{k+1-m}$ .

5. The a priori error bound. Combining 3 and 4 with $H=H^{1}_0(\Omega)$ and the Dirichlet form $a$ : if the exact solution lies in $H^{2}(\Omega)$ then

\lVert u-u_h\rVert_{H^{1}(\Omega)}\ \le\ C' h\,\lvert u\rvert_{H^{2}(\Omega)} .

That $u\in H^{2}$ is a hypothesis about the domain, not a free gift: Evans's elliptic regularity theorem supplies it when $f\in L^{2}$ and $\partial\Omega$ is of class $C^{2}$ , and it is also true for a convex polygon, but it fails at a re-entrant corner.

Idea of the proofOriginal

Everything turns on one observation. The exact solution satisfies the weak equation against every test function; the approximation satisfies it against the test functions we kept. Subtracting, the error is invisible to the whole of $V_h$ : that is Galerkin orthogonality, and it is the statement that $u_h$ is a projection of $u$ along $V_h$ .

A projection cannot be much worse than the best approximation, and in an inner-product space it is the best approximation. Coercivity and boundedness say that $a$ is a distorted inner product, distorted by no more than the factor $M/m$ , so the Galerkin solution is within that factor of best. This is Cea's lemma, and it splits the problem cleanly in two: an analysis question (how good is the projection?) answered once and for all, and an approximation-theory question (how well can piecewise-linear functions approximate $u$ ?) answered by Taylor's theorem on each triangle.

The second question is where the rate comes from. On a triangle of size $h$ the linear interpolant of a $C^{2}$ function is wrong by $O(h^{2})$ and its gradient by $O(h)$ ; shape regularity is exactly the condition that stops the gradient estimate from degenerating when a triangle is thin. Multiply the two answers and the energy-norm error is $O(h)$ .

Proof(Parts 2 and 3; part 4 by Taylor's theorem)AdaptedAdapted from Stakgold & Holst (2011), §10.2, proof of Theorem 10.2.9 and Lemma 10.2.3, pp. 678–679

2. The exact solution satisfies $a(u,v)=\ell(v)$ for every $v\in H$ , and in particular for every $v_h\in V_h$ , because $V_h\subset H$ . The discrete solution satisfies $a(u_h,v_h)=\ell(v_h)$ for the same $v_h$ . Subtracting and using bilinearity in the first slot gives $a(u-u_h,v_h)=0$ . The hypothesis used is only the conformity of the method, $V_h\subset H$ ; nothing else.

3. Let $v_h\in V_h$ be arbitrary. Coercivity applied to the error, then part 2 with the test function $u_h-v_h\in V_h$ , then boundedness:

m\lVert u-u_h\rVert^{2} \ \le\ a(u-u_h,\,u-u_h) \ =\ a\bigl(u-u_h,\,u-u_h\bigr)+a\bigl(u-u_h,\,u_h-v_h\bigr)

=\ a\bigl(u-u_h,\,u-v_h\bigr) \ \le\ M\,\lVert u-u_h\rVert\,\lVert u-v_h\rVert .

The added term vanishes by part 2 since $u_h-v_h\in V_h$ ; the two arguments then combine by bilinearity. If $\lVert u-u_h\rVert=0$ the claim is trivial, and otherwise dividing by $m\lVert u-u_h\rVert$ gives $\lVert u-u_h\rVert\le (M/m)\lVert u-v_h\rVert$ . Taking the infimum over $v_h$ gives part 3.

Stakgold and Holst reach the constant $1+M/m$ by bounding the norm of the generalized projection $P_h$ and writing $u-u_h=(I-P_h)(u-v_h)$ , then remove the leading $1$ in the Hilbert-space case by the Kato–Xu–Zikatanov lemma, which says that a non-trivial idempotent $P\ne 0,I$ on a Hilbert space satisfies $\lVert P\rVert=\lVert I-P\rVert$ . The direct argument above gives $M/m$ at once in the Hilbert setting, and is the one we use.

Symmetric case. If $a(v,w)=a(w,v)$ then coercivity and boundedness make $a$ an inner product whose norm $\lVert\cdot\rVert_a$ is equivalent to $\lVert\cdot\rVert$ . Part 2 says $u-u_h\perp_a V_h$ , so for every $v_h$ , by Pythagoras in that inner product,

\lVert u-v_h\rVert_a^{2}=\lVert u-u_h\rVert_a^{2}+\lVert u_h-v_h\rVert_a^{2}\ \ge\ \lVert u-u_h\rVert_a^{2} ,

with equality exactly at $v_h=u_h$ . So in the energy norm the Galerkin solution is the best approximation; the factor $M/m$ in part 3 is the price of measuring in a different norm.

4, in one dimension and in outline. On an interval $[x_{j},x_{j+1}]$ of length $h$ let $e=w-I_hw$ , which vanishes at both endpoints. By Rolle's theorem $e'$ vanishes at some interior $\xi$ , so for any $x$ in the interval $e'(x)=\int_\xi^x w''$ , and by Cauchy–Schwarz $\lvert e'(x)\rvert^{2}\le h\int\lvert w''\rvert^{2}$ ; integrating over the interval gives $\int\lvert e'\rvert^{2}\le h^{2}\int\lvert w''\rvert^{2}$ . Since $e$ vanishes at $x_j$ , $\lvert e(x)\rvert^2\le h\int\lvert e'\rvert^{2}$ in the same way, and a second integration gives $\int\lvert e\rvert^{2}\le h^{4}\int\lvert w''\rvert^{2}$ . Summing over the elements gives part 4 with $C=1$ in one dimension. In two dimensions the same two powers of $h$ appear, by mapping each triangle affinely to a fixed reference triangle, estimating there by the Bramble–Hilbert argument, and mapping back; the Jacobian factors produce $h_K^{k+1}/\rho_K^{m}$ , which is where shape regularity enters. Stakgold and Holst give this chain as Lemmas 10.5.1–10.5.3 and Theorem 10.5.2 and attribute the detailed proof to Ciarlet. $\square$

Remark(Where the hypotheses are used, and what goes wrong without them)Original

Conformity ( $V_h\subset H^{1}_0$ ) is the only hypothesis in Galerkin orthogonality, and it is why continuity across edges matters: a piecewise-linear function with a jump across an edge is not in $H^1$ , its weak gradient picks up a delta on the edge, and the whole argument collapses. Non-conforming elements are used in practice, but they need their own error analysis with an extra consistency term.

Coercivity is what makes $A$ invertible and the error bound finite. For the Dirichlet form it comes from the Poincare inequality (The Lax-Milgram theorem and the Dirichlet problem), which uses the boundary condition; drop the boundary condition and $a(1,1)=0$ , the constant function is in the kernel, and the pure Neumann problem is solvable only up to a constant and only for $f$ with mean zero.

Shape regularity is used only in part 4, and the bound genuinely fails without it. Take the flat triangle $K$ with vertices $(0,0)$ , $(h,0)$ , $(h/2,\varepsilon)$ , whose largest angle tends to $\pi$ as $\varepsilon\to0$ , and interpolate $w(x,y)=x^{2}$ . The affine function through the three vertex values $0$ , $h^{2}$ , $h^{2}/4$ is $p(x,y)=hx-\tfrac{h^{2}}{4\varepsilon}y$ , so $\partial_y(w-p)=h^{2}/(4\varepsilon)$ on all of $K$ and

\lvert w-I_Kw\rvert_{H^1(K)}^{2}\ \ge\ \frac{h^{4}}{16\varepsilon^{2}}\cdot\mathrm{meas}(K)=\frac{h^{5}}{32\varepsilon}, \qquad h_K^{2}\,\lvert w\rvert^{2}_{H^{2}(K)}=h^{2}\cdot 4\,\mathrm{meas}(K)=2h^{3}\varepsilon .

Their ratio is $h^{2}/(64\varepsilon^{2})$ , which is unbounded: the clean estimate $\lvert w-I_Kw\rvert_{H^1}\le Ch_K\lvert w\rvert_{H^2}$ is false for a degenerating family. The estimate that survives is the one Stakgold and Holst actually state, with $h_K^{2}/\rho_K$ in place of $h_K$ ; here $\rho_K$ is of order $\varepsilon$ , that factor is of order $h^{2}/\varepsilon$ , and the bound it gives, of order $h^{5}/\varepsilon$ for the square of the seminorm, matches the true error. So the error is not mysteriously large — it is exactly as large as the general theorem says, and shape regularity is the hypothesis that turns the general theorem into a rate in $h$ . A mesh refined by flattening rather than shrinking its triangles need not converge in $H^1$ at all. Quasi-uniformity, the stronger condition $h_{K_1}/h_{K_2}\le\nu$ that Stakgold and Holst state alongside, is not needed for the bound above; it is deliberately broken by adaptive methods, which refine only where an error indicator is large.

$u\in H^{2}$ is used in part 5 and nowhere else. On a convex polygon it holds; on an L-shaped domain the solution behaves like $r^{2/3}\sin(2\theta/3)$ near the re-entrant corner, which is in $H^1$ but not in $H^2$ , and the uniform-mesh rate drops from $h$ to $h^{2/3}$ in the energy norm. Recovering the full rate there requires grading the mesh towards the corner — which is precisely the freedom that shape regularity without quasi-uniformity leaves open.

What quasi-optimality does not say. Part 3 bounds the error by the best approximation error; it says nothing about whether the best approximation is any good. All of the approximation power of the method is in part 4, and all of the geometry is in the triangulation. This separation is the reason the finite element method can be analysed at all on a domain whose shape no coordinate system fits.

Method(Assembling and solving)Original

Triangulate. Produce a conforming, shape-regular mesh of $\overline\Omega$ ; record the vertex coordinates and, for each triangle, the indices of its three vertices. Mark the vertices on $\partial\Omega$ as constrained.
Compute the element matrices. On a triangle $K$ with vertices $P_1,P_2,P_3$ , the hat function $\lambda_r$ associated with $P_r$ is affine with $\lambda_r(P_s)=\delta_{rs}$ , so $\nabla\lambda_r$ is a constant vector on $K$ . Explicitly, with $\mathrm{area}(K)=\tfrac12\lvert\det[P_2-P_1,\,P_3-P_1]\rvert$ , $\nabla\lambda_r=\frac{1}{2\,\mathrm{area}(K)}\,\bigl(P_{r+2}-P_{r+1}\bigr)^{\perp},$ where the indices are read modulo $3$ in $\{1,2,3\}$ , $(a,b)^{\perp}=(-b,a)$ , and $P_1,P_2,P_3$ are listed anticlockwise. (Check on the triangle $(0,0),(1,0),(0,1)$ : $\lambda_1=1-x-y$ , and the formula gives $\bigl((0,1)-(1,0)\bigr)^{\perp}=(-1,-1)$ .) The $3\times3$ element stiffness matrix is then $A^K_{rs}=\mathrm{area}(K)\,\nabla\lambda_r\cdot\nabla\lambda_s$ , exactly, because the integrand is constant.
Assemble. Initialise $A=0$ , $F=0$ and loop over triangles, adding $A^{K}_{rs}$ into the global entry $(i_r,i_s)$ for each pair of local vertices and $\int_K f\lambda_r$ into $F_{i_r}$ . Skip rows and columns belonging to constrained vertices (or keep them and impose $U=0$ there). The whole loop is $O(\#\mathcal{T}_h)=O(N)$ .
Integrate the load. $\int_K f\lambda_r$ is rarely available in closed form; the usual choice is a quadrature rule exact for the polynomial degree in play (Numerical quadrature). For linear elements, the one-point rule $\int_K g\approx \mathrm{area}(K)\,g(\text{centroid})$ is exact for affine $g$ and enough to preserve the $O(h)$ rate. For $f\equiv1$ it gives the exact value $\int_K\lambda_r=\mathrm{area}(K)/3$ .
Solve. $A$ is symmetric positive definite and sparse, with about $7$ non-zeros per row for a typical triangulation. Use the Cholesky factorisation of a banded or reordered matrix for moderate $N$ , or the conjugate gradient method preconditioned so that the iteration count does not grow with $N$ (Solving linear systems numerically, Relaxation methods for elliptic problems). An unpreconditioned iteration is slow because $\operatorname{cond}(A)=O(h^{-2})$ , the same growth as for the difference Laplacian.
Check the rate. Solve on a mesh and on its uniform refinement and compare: the energy-norm error should fall by $2$ and the $L^{2}$ error by $4$ . A rate below that is the usual sign that the solution is not in $H^{2}$ , that the quadrature in step 4 is too crude, or that the mesh is degenerating.
For a time-dependent problem keep the same space and leave the time variable continuous: $K\,\dot U = -\kappa A\,U + F(t)$ with $K$ the mass matrix $\int v_iv_j$ , then step in time by any method of Integrators for ODEs. Strauss writes out both the forward and the backward Euler versions; the backward one requires solving with $K+\kappa\,\Delta t\,A$ at each step, which is the price of avoiding the $\Delta t=O(h^{2})$ restriction of Stability and stiff equations.

Example(One dimension exactly, two dimensions on four triangles)Original

(a) A one-dimensional problem solved completely. Take

-u''=f \ \text{ on } (0,1),\qquad u(0)=u(1)=0,\qquad f(x)=\pi^{2}\sin\pi x,

whose exact solution is $u(x)=\sin\pi x$ . Use the uniform mesh $x_j=jh$ , $h=1/n$ , with the tent functions $v_j$ of the definition. Two integrals are needed.

Stiffness. $v_j'=1/h$ on $(x_{j-1},x_j)$ and $-1/h$ on $(x_j,x_{j+1})$ , so

A_{jj}=\int_0^1 (v_j')^{2}=2h\cdot\frac1{h^{2}}=\frac2h, \qquad A_{j,j+1}=\int_{x_j}^{x_{j+1}}\!\!\Bigl(-\frac1h\Bigr)\Bigl(\frac1h\Bigr)=-\frac1h ,

and $A_{ji}=0$ for $\lvert i-j\rvert\ge2$ . So $A=h^{-1}\mathrm{tridiag}(-1,2,-1)$ , which is Strauss's tridiagonal matrix with $2$ on the diagonal and $-1$ beside it, divided by $h$ .

Load. With $\int_{0}^{h}\frac{t}{h}\sin\bigl(\pi(x_j-h+t)\bigr)\mathrm{d}t$ and its mirror image combined,

F_j=\pi^{2}\!\int_0^1 \sin(\pi x)\,v_j(x)\,\mathrm{d}x=\frac{2\bigl(1-\cos\pi h\bigr)}{h}\,\sin(\pi x_j) .

Solving. Try $U_j=c\sin(\pi x_j)$ . The left side of the $j$ th equation is

\frac{c}{h}\bigl(2\sin\pi x_j-\sin\pi x_{j-1}-\sin\pi x_{j+1}\bigr) =\frac{c}{h}\bigl(2-2\cos\pi h\bigr)\sin\pi x_j ,

using $\sin(\theta\pm\pi h)=\sin\theta\cos\pi h\pm\cos\theta\sin\pi h$ , whose $\cos\theta$ terms cancel. Comparing with $F_j$ gives $c=1$ . So

U_j=\sin(\pi x_j)=u(x_j)\quad\text{exactly, at every node and for every }h .

With $h=\tfrac14$ : $A=4\,\mathrm{tridiag}(-1,2,-1)$ , $F=\bigl(8-4\sqrt2\bigr)\bigl(\tfrac{\sqrt2}{2},1,\tfrac{\sqrt2}{2}\bigr)=(1.656854,\,2.343146,\,1.656854)$ , and the solution of the $3\times3$ system is $U=(\tfrac{\sqrt2}{2},1,\tfrac{\sqrt2}{2})$ .

Why exact? Not an accident of this $f$ . Fix a node $x_k$ and let $G_k$ be the Green's function $G(\cdot,x_k)$ of the problem, so $u(x_k)=\int_0^1 G_k f$ and $a(w,G_k)=w(x_k)$ for every $w\in H^1_0$ (Green's function for a boundary-value problem). For $-\mathrm{d}^2/\mathrm{d}x^2$ on an interval, $G_k$ is piecewise linear with its only kink at $x_k$ — a node — so $G_k\in V_h$ . Galerkin orthogonality then gives $u(x_k)-u_h(x_k)=a(u-u_h,\,G_k)=0$ . In two dimensions the Green's function has a logarithmic singularity and is not in $V_h$ , and nodal exactness is lost.

The error between the nodes is not zero, and this is what carries the rate. Computing $\lVert u-u_h\rVert_{L^{2}}$ and $\lvert u-u_h\rvert_{H^{1}}$ exactly, element by element:

$n$	$h$	$\lVert u-u_h\rVert_{L^2}$	ratio	$\lvert u-u_h\rvert_{H^1}$	ratio
$4$	$0.25$	$3.928435\times10^{-2}$	—	$0.4985085$	—
$8$	$0.125$	$9.920920\times10^{-3}$	$3.960$	$0.2511818$	$1.9847$
$16$	$0.0625$	$2.486501\times10^{-3}$	$3.990$	$0.1258332$	$1.9962$

The ratios approach $4$ and $2$ : $O(h^{2})$ in $L^{2}$ and $O(h)$ in $H^{1}$ , as parts 4 and 5 predict. (Here $\lvert u\rvert_{H^{2}}=\pi^{2}/\sqrt2=6.9789$ , and the measured $H^1$ error at $h=1/16$ is $0.1258$ , so the constant $C'$ in part 5 is about $0.29$ .)

(b) Two dimensions on the coarsest possible mesh. Take $\Omega=(0,1)^{2}$ , $f\equiv1$ , $u=0$ on the boundary. Triangulate by drawing both diagonals of the square: four triangles, four boundary vertices, and one interior vertex $V_1=(\tfrac12,\tfrac12)$ . So $N=1$ and the system is a single equation.

On the lower triangle, with vertices $(0,0),(1,0),(\tfrac12,\tfrac12)$ , the hat function $v_1$ vanishes at the first two and equals $1$ at the third, and the affine function with those values is $v_1(x,y)=2y$ . Hence $\lvert\nabla v_1\rvert^{2}=4$ there, and the triangle has area $\tfrac14$ , contributing $4\cdot\tfrac14=1$ . By the fourfold symmetry of the mesh the other three contribute the same, so

A_{11}=\int_\Omega\lvert\nabla v_1\rvert^{2}=4 .

For the load, the integral of an affine function over a triangle is its area times the mean of its three vertex values, so each triangle contributes $\tfrac14\cdot\tfrac{0+0+1}{3}=\tfrac1{12}$ and

F_1=\int_\Omega v_1 = 4\cdot\tfrac1{12}=\tfrac13 .

Therefore

u_h\Bigl(\tfrac12,\tfrac12\Bigr)=U_1=\frac{F_1}{A_{11}}=\frac{1/3}{4}=\frac{1}{12}=0.0833333 .

The exact value, from the double sine series $u=\sum_{m,n\ \mathrm{odd}} \dfrac{16}{\pi^{4}mn(m^{2}+n^{2})}\sin m\pi x\,\sin n\pi y$ of Separation of variables in standard coordinates, is

u\Bigl(\tfrac12,\tfrac12\Bigr)=\sum_{m,n\ \mathrm{odd}}\frac{16\,(-1)^{(m+n)/2-1}}{\pi^{4}mn(m^{2}+n^{2})}=0.0736714 ,

so the single-element answer is high by $13.1\%$ . That is a respectable showing for one unknown.

Refining. Take instead the uniform mesh of $n^{2}$ squares each cut by its lower-left-to-upper-right diagonal. A short computation of the six triangles around an interior node (the same calculation as above, done on right triangles with legs $h$ ) gives the row

A_{ii}=4,\qquad A_{ij}=-1 \ \text{for the four axis neighbours},\qquad A_{ij}=0 \ \text{for the two diagonal ones},

with no $h$ in it, while $F_i=\int_\Omega v_i = h^{2}$ for $f\equiv1$ . So on this mesh the finite element equations are identical to the five-point difference scheme of Finite-difference approximations — the methods coincide here, and they diverge as soon as the mesh or the data stop being so symmetric. Solving:

$n$	$h$	$u_h(\tfrac12,\tfrac12)$	error	ratio
$2$	$0.5$	$0.0625000$	$-1.1171\times10^{-2}$	—
$4$	$0.25$	$0.0703125$	$-3.3589\times10^{-3}$	$3.33$
$8$	$0.125$	$0.0727826$	$-8.8872\times10^{-4}$	$3.78$
$16$	$0.0625$	$0.0734458$	$-2.2559\times10^{-4}$	$3.94$
$32$	$0.03125$	$0.0736147$	$-5.6616\times10^{-5}$	$3.98$

The nodal error falls by a factor approaching $4$ : second order, one power better than the $O(h)$ of part 5, which is a bound in the energy norm. The $L^{2}$ and pointwise errors of the finite element method on a convex domain are $O(h^{2})$ , by a duality argument that trades the extra power for the regularity of an auxiliary problem.

Application(Why engineering computes this way)Original

Structural analysis, where the names come from. The method was built by aeronautical and civil engineers in the 1950s to compute stresses in built-up structures, and the vocabulary is theirs: $A$ is the stiffness matrix because for the elasticity bilinear form $a(u,v)=\int \sigma(u):\varepsilon(v)$ it is literally the matrix of spring constants relating nodal displacements to nodal forces, and $F$ is the load. The energy-norm statement of part 3 then has a physical reading: among all displacement fields representable on the mesh, the computed one is the one of least potential energy, which is the Rayleigh–Ritz principle of The Rayleigh-Ritz method restricted to a finite-dimensional space. Strauss notes that the method was invented to deal with curved and irregular domains, where finite differences fit badly.

Why it wins over finite differences. Not accuracy on a square — the example above shows the two methods coinciding there — but geometry and boundary conditions. An aerofoil, a heart ventricle, a turbine blade and a reservoir have no coordinate system; a triangulation has no trouble with any of them, grades itself finely where the solution varies and coarsely elsewhere, and handles a Neumann or Robin condition for free, because such conditions appear as boundary terms in the weak form rather than as extra equations at ghost points.

Where the error estimate is used in anger. Part 5 bounds the error by a norm of the unknown solution, which cannot be evaluated. Practical codes therefore use a posteriori estimators, computable from $u_h$ alone — typically the element residual $\lVert f+\Delta u_h\rVert_{L^2(K)}$ together with the jumps of $\nabla u_h$ across edges — to decide which elements to refine, and iterate solve, estimate, mark, refine. Stakgold and Holst give exactly this loop as an algorithm. The theory of this unit is what guarantees that the loop is converging to something.

Beyond Poisson. The same three lines — weak form, conforming subspace, Galerkin — carry over to linear elasticity, to Stokes flow (where the velocity and pressure spaces must be chosen compatibly, or the inf–sup condition of part 1 fails and the pressure develops a checkerboard oscillation), to the eigenvalue problems of Sturm-Liouville eigenvalue problems in several dimensions, and to the time-dependent problems of step 7. The analysis of all of them is the same analysis: Cea's lemma plus an approximation estimate.

Method

Monte Carlo methods

Adapted

Add to PDF

Turns an integral into the mean of a random variable, so that the law of large numbers computes it and the central limit theorem supplies an error bar; shows that the root-mean-square error is $\sigma/\sqrt{n}$ whatever the dimension, which is why the method is used above four dimensions and not below; and derives importance sampling and the Metropolis algorithm as the two ways of changing the sampling distribution when the obvious one is too noisy or unavailable.

Prerequisites: Laws of large numbers, The central limit theorem, Convergence of Markov chains, Expectation, variance and covariance

Definition(The estimator, and the three ways of choosing the samples)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.4.4, Eqs. (27.48)–(27.49), pp. 1011–1012; Grimmett & Stirzaker (2020), §2.6, Example (3), p. 44; Grimmett & Stirzaker (2020), Exercise (4.11.12), pp. 142–143; Grimmett & Stirzaker (2020), §6.14, Eqs. (6.14.3)–(6.14.5) and Example (6.14.8), pp. 331–332

Every integral is an expectation. If $D\subseteq\mathbb{R}^{d}$ has finite volume $\lvert D\rvert$ and $X$ is uniform on $D$ , then

\theta=\int_D g(x)\,\mathrm{d}x=\lvert D\rvert\,\mathrm{E}\bigl[g(X)\bigr] ,

and more generally, for any probability density $p$ that is positive wherever $g$ is non-zero,

\theta=\int_D g(x)\,\mathrm{d}x=\mathrm{E}_p\!\left[\frac{g(Y)}{p(Y)}\right],\qquad Y\sim p .

A Monte Carlo method estimates $\theta$ by drawing samples and averaging. Three choices of sampling distribution give the three methods of this unit. Throughout, $X_1,X_2,\dots$ are independent and $\xi_1,\xi_2,\dots$ are independent and uniform on $[0,1)$ .

1. Crude Monte Carlo. Sample uniformly. Riley, Hobson and Bence first reduce the problem by a change of variable to the standard form $\theta=\int_0^1 f(x)\,\mathrm{d}x$ and then estimate it by

\hat\theta_n=\frac1n\sum_{i=1}^{n}f(\xi_i) .

They insist on two properties of any such scheme: the estimator must be unbiased, $\mathrm{E}[\hat\theta_n]=\theta$ , and the method must report the variance of the estimate — a statistical statement, not an interval inside which the answer is guaranteed to lie.

A variant, the hit-or-miss method, estimates the area under a curve by the proportion of random points in a bounding box that fall below it: if $0\le f\le1$ on $[0,1]$ and $(\xi_{2i-1},\xi_{2i})$ are uniform on the unit square, the fraction of $i$ with $\xi_{2i}\le f(\xi_{2i-1})$ estimates $\theta$ , because that event has probability exactly $\theta$ . It is simpler and, as the remark below shows, always noisier.

2. Importance sampling. Sample from a density $p$ of our choosing and correct by the likelihood ratio. For $\theta=\mathrm{E}[g(X)]$ with $X$ of density $q$ , choose a density $p$ that vanishes exactly where $q$ does and put

\hat\theta_n^{\,\mathrm{IS}}=\frac1n\sum_{i=1}^{n} w(Y_i)\,g(Y_i), \qquad w(y)=\frac{q(y)}{p(y)}, \qquad Y_i\sim p \text{ independent} .

Grimmett and Stirzaker call $p$ the importance density. The point is to put the samples where $\lvert g\rvert q$ is large; the price is that one must be able to sample from $p$ .

3. Markov chain Monte Carlo. When $\pi$ is a distribution on a finite set $\Theta$ known only up to its normalising constant, as a Boltzmann weight $\pi(\vartheta)=Z^{-1}\mathrm{e}^{-\beta E(\vartheta)}$ is, independent sampling from $\pi$ is out of reach. Instead build a Markov chain with stationary distribution $\pi$ and average along a single trajectory. Grimmett and Stirzaker's recipe has two ingredients: a proposal matrix $H=(h_{ij})$ , any stochastic matrix that is cheap to simulate, and an acceptance matrix $A=(a_{ij})$ with $0\le a_{ij}\le1$ . From state $i$ one draws $j$ with probability $h_{ij}$ and moves there with probability $a_{ij}$ , staying put otherwise, so that

p_{ij}=h_{ij}a_{ij}\ (i\ne j), \qquad p_{ii}=1-\sum_{k\ne i}h_{ik}a_{ik} .

The Hastings choice of acceptance probability is

a_{ij}=\min\left\{1,\ \frac{\pi_j h_{ji}}{\pi_i h_{ij}}\right\} ,

and when the proposal is symmetric, $h_{ij}=h_{ji}$ , this reduces to the Metropolis rule

a_{ij}=\min\left\{1,\ \frac{\pi_j}{\pi_i}\right\} .

Both depend on $\pi$ only through ratios, so the unknown constant $Z$ never appears. The estimator, allowing a burn-in of $M$ discarded steps, is

\hat\theta^{\,\mathrm{MC}}_{M,n}=\frac1n\sum_{r=M}^{M+n-1} g(X_r) .

Method(What the three estimators cost)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.4.4, pp. 1011–1012; Grimmett & Stirzaker (2020), Exercise (4.11.12), pp. 142–143; Grimmett & Stirzaker (2020), §6.14, Eqs. (6.14.3)–(6.14.5), pp. 331–332

1. Crude Monte Carlo. Let $f$ be measurable on $[0,1]$ with $\theta=\int_0^1f$ and $\sigma^{2}=\int_0^1 f^{2}-\theta^{2}<\infty$ . Then for every $n$

\mathrm{E}\bigl[\hat\theta_n\bigr]=\theta, \qquad \operatorname{Var}\bigl[\hat\theta_n\bigr]=\frac{\sigma^{2}}{n}, \qquad \sqrt{\mathrm{E}\bigl[(\hat\theta_n-\theta)^{2}\bigr]}=\frac{\sigma}{\sqrt n} .

Moreover $\hat\theta_n\to\theta$ almost surely; and if in addition $\sigma>0$ , then $\sqrt n(\hat\theta_n-\theta)/\sigma$ converges in distribution to a standard normal, so that

\Pr\Bigl(\lvert\hat\theta_n-\theta\rvert\le 1.96\,\frac{s_n}{\sqrt n}\Bigr)\longrightarrow 0.95, \qquad s_n^{2}=\frac{1}{n-1}\sum_{i=1}^{n}\bigl(f(\xi_i)-\hat\theta_n\bigr)^{2} .

(The hypothesis $\sigma>0$ is not cosmetic: for $f$ constant the ratio is $0/0$ , $s_n$ vanishes identically, and the probability on the left is $1$ for every $n$ , not $0.95$ .)

2. Dimension. None of this mentions the dimension. In $d$ dimensions, with $\theta=\int_{[0,1]^d}f$ and $\sigma^{2}=\operatorname{Var}f(X)$ for $X$ uniform on the cube, the same three statements hold verbatim. By contrast a product quadrature rule of order $p$ (Numerical quadrature) with $k$ nodes per axis uses $n=k^{d}$ points and has error $O(k^{-p})=O(n^{-p/d})$ . Hence

\text{Monte Carlo is asymptotically the better rate exactly when } d>2p .

For the trapezoidal rule ( $p=2$ ) the crossover is $d=4$ ; for Simpson's rule ( $p=4$ ) it is $d=8$ .

3. Importance sampling. With $q$ , $p$ , $w=q/p$ and $Y_i\sim p$ as in the definition, and assuming $\mathrm{E}_p[(wg)^{2}]<\infty$ ,

\mathrm{E}\bigl[\hat\theta^{\,\mathrm{IS}}_n\bigr]=\theta, \qquad \operatorname{Var}\bigl[\hat\theta^{\,\mathrm{IS}}_n\bigr]=\frac{1}{n}\left(\int \frac{g^{2}q^{2}}{p}-\theta^{2}\right),

and $\hat\theta^{\,\mathrm{IS}}_n\to\theta$ almost surely. The variance is minimised, over all densities $p$ that are positive wherever $gq\neq0$ , by

p^{\ast}(y)=\frac{\lvert g(y)\rvert q(y)}{\int\lvert g\rvert q} ,

for which the variance is $\bigl(\int\lvert g\rvert q\bigr)^{2}-\theta^{2}$ , and this is zero when $g$ does not change sign. The class competed over here is wider than the one the definition allows: $p^{\ast}$ vanishes wherever $g$ does, so if $g=0$ on a set where $q>0$ then $p^{\ast}$ is not a density that "vanishes exactly where $q$ does". That is harmless — the estimator needs $p>0$ only where $gq\neq0$ , since the summand $wg$ is zero elsewhere — but the two conditions are different, and it is the weaker one that the optimality is proved over. The ideal is unusable — it needs the answer — but it says what to aim at: make $w g$ nearly constant.

4. Markov chain Monte Carlo. Let $\pi$ be positive on a finite $\Theta$ , let $H$ be a stochastic matrix with $h_{ij}>0$ exactly when $h_{ji}>0$ , and let $P$ be built with the Hastings acceptance probabilities. Then $P$ satisfies the detailed balance equations $\pi_ip_{ij}=\pi_jp_{ji}$ , so $\pi$ is stationary for $P$ and the chain is reversible. If $P$ is in addition irreducible and aperiodic then, by Convergence of Markov chains,

\hat\theta^{\,\mathrm{MC}}_{M,n}\longrightarrow\sum_{\vartheta}g(\vartheta)\pi(\vartheta) \quad\text{almost surely as } n\to\infty,\ \text{for every } M \text{ and every starting state},

and the distribution of $X_n$ approaches $\pi$ geometrically at the rate $\lvert\lambda_2\rvert$ , the second largest eigenvalue modulus of $P$ . For large $n$ the variance of the estimator is

\operatorname{Var}\bigl[\hat\theta^{\,\mathrm{MC}}_{M,n}\bigr]\ \sim\ \frac{\gamma_0\,\tau}{n}, \qquad \tau=1+2\sum_{k=1}^{\infty}\rho_k, \qquad \rho_k=\frac{\gamma_k}{\gamma_0},\qquad \gamma_k=\operatorname{Cov}_\pi\bigl(g(X_0),g(X_k)\bigr) ,

so $n/\tau$ correlated samples are worth $n$ independent ones. The quantity $\tau$ is the integrated autocorrelation time; it exceeds $1$ when successive states are positively correlated, which is the usual case, and can be less than $1$ .

Idea of the proofOriginal

The law of large numbers is a numerical method. Write the quantity you want as an average over a probability distribution, draw from that distribution, and average: convergence is guaranteed by Laws of large numbers, and the central limit theorem (The central limit theorem) converts the convergence into an error bar computed from the same samples that produced the answer. Nothing in that argument knows how many variables there are, which is the whole reason the method exists: a product quadrature rule must pay $k^{d}$ points to resolve $k$ scales per axis, while a Monte Carlo estimate pays $\sigma^{2}/\epsilon^{2}$ points for accuracy $\epsilon$ no matter what $d$ is.

The price is the rate: $n^{-1/2}$ is slow, and nothing improves it, because it is the central limit theorem. What can be improved is $\sigma$ , and the two refinements do exactly that. Importance sampling changes the measure so that the quantity being averaged is more nearly constant, which is the only thing that makes a variance small. The Metropolis algorithm changes the measure in the opposite situation, where the target distribution is known only up to a constant: it gives up independence, building instead a chain whose equilibrium is the target, and buys that with a loss of efficiency measured by the autocorrelation time.

The acceptance formula is not guessed. Demand detailed balance — that the chain in equilibrium runs the same backwards as forwards — write out what that requires of $a_{ij}$ , and the condition is a single equation in the ratio $a_{ij}/a_{ji}$ . Taking each $a$ as large as that equation allows gives the minimum, and it involves $\pi$ only through ratios, which is why an unknown partition function is no obstacle.

Proof(Parts 1–4)AdaptedAdapted from Grimmett & Stirzaker (2020), Exercise (4.11.12), pp. 142–143; Grimmett & Stirzaker (2020), §6.14, Eqs. (6.14.3)–(6.14.5), p. 332

1. Each $f(\xi_i)$ has mean $\theta$ and variance $\sigma^{2}$ , so $\mathrm{E}[\hat\theta_n]=\theta$ by linearity. Since the $\xi_i$ are independent, the variances add (Expectation, variance and covariance):

\operatorname{Var}\bigl[\hat\theta_n\bigr]=\frac1{n^{2}}\sum_{i=1}^{n}\operatorname{Var}\bigl[f(\xi_i)\bigr]=\frac{\sigma^{2}}{n} .

Because the estimator is unbiased, its mean square error is its variance, which gives the third formula. Almost sure convergence is the strong law of Laws of large numbers, which needs only $\mathrm{E}\lvert f(\xi)\rvert<\infty$ — weaker than the variance hypothesis. The normal limit is The central limit theorem applied to the i.i.d. summands $f(\xi_i)$ ; replacing $\sigma$ by the sample standard deviation $s_n$ is legitimate because $s_n\to\sigma$ almost surely and Slutsky's lemma lets a convergent-in-probability factor be substituted into a weak limit.

2. The computation in part 1 used independence and the two moments of $f(\xi)$ and nothing else, so it is insensitive to the domain. For the comparison, a product rule of order $p$ has error $Ck^{-p}$ per axis with $k$ nodes per axis, hence $n=k^{d}$ and error $Cn^{-p/d}$ . Setting $n^{-p/d}=n^{-1/2}$ gives $p/d=1/2$ ; for $d>2p$ the exponent $p/d$ is smaller than $1/2$ and the deterministic rule is asymptotically worse.

3. Since $p>0$ wherever $q>0$ , the substitution is legitimate and

\mathrm{E}\bigl[w(Y)g(Y)\bigr]=\int \frac{q(y)}{p(y)}g(y)\,p(y)\,\mathrm{d}y=\int g\,q=\theta ,

so each summand is unbiased and so is the average. The summands are independent with common variance $\mathrm{E}_p[(wg)^{2}]-\theta^{2}=\int g^{2}q^{2}/p-\theta^{2}$ , and dividing by $n$ gives the stated variance; the strong law gives almost sure convergence. For the optimal density, apply the Cauchy–Schwarz inequality in the form

\left(\int \lvert g\rvert q\right)^{2}=\left(\int \frac{\lvert g\rvert q}{\sqrt p}\cdot\sqrt p\,\right)^{2}\le \int\frac{g^{2}q^{2}}{p}\cdot\int p=\int\frac{g^{2}q^{2}}{p} ,

using $\int p=1$ . So $\int g^{2}q^{2}/p$ is at least $(\int\lvert g\rvert q)^{2}$ for every $p$ that is positive wherever $gq\neq0$ (read $g^{2}q^{2}/p$ as $0$ on the set where $gq$ and $p$ both vanish, which is what the estimator does), and equality holds exactly when $\lvert g\rvert q/\sqrt p$ is a constant multiple of $\sqrt p$ , that is when $p\propto\lvert g\rvert q$ . Substituting $p=p^{\ast}$ gives the variance $(\int\lvert g\rvert q)^{2}-\theta^{2}$ , which vanishes when $\lvert\int gq\rvert=\int\lvert g\rvert q$ , that is when $g$ has one sign.

4. Detailed balance forces the acceptance rule. For $i\ne j$ the requirement $\pi_ip_{ij}=\pi_jp_{ji}$ reads

\pi_i h_{ij}a_{ij}=\pi_j h_{ji}a_{ji}, \qquad\text{that is}\qquad \frac{a_{ij}}{a_{ji}}=\frac{\pi_j h_{ji}}{\pi_i h_{ij}}=:R_{ij} ,

with $R_{ji}=1/R_{ij}$ . Any pair $(a_{ij},a_{ji})$ in $[0,1]^{2}$ with that ratio will do; taking the larger of the two equal to $1$ gives the largest admissible pair, namely $a_{ij}=\min\{1,R_{ij}\}$ and $a_{ji}=\min\{1,1/R_{ij}\}$ — exactly the Hastings formula, which is therefore the choice that rejects as rarely as possible. The diagonal entries satisfy detailed balance trivially. Since detailed balance implies stationarity (sum $\pi_ip_{ij}=\pi_jp_{ji}$ over $i$ and use $\sum_ip_{ji}=1$ ), $\pi$ is stationary and the chain is reversible; both statements, and the convergence of the $n$ -step probabilities and the geometric rate $\lvert\lambda_2\rvert$ , are proved in Convergence of Markov chains. The almost sure convergence of the time average is the ergodic theorem for an irreducible positive-recurrent chain; it holds for every starting state, so the burn-in $M$ affects the finite- $n$ bias but not the limit.

The autocorrelation formula. In equilibrium $\{g(X_r)\}$ is stationary (Stationary processes and spectral densities), so $\operatorname{Cov}(g(X_r),g(X_s))=\gamma_{\lvert r-s\rvert}$ and

\operatorname{Var}\Bigl[\frac1n\sum_{r=0}^{n-1}g(X_r)\Bigr] =\frac{1}{n^{2}}\sum_{r,s=0}^{n-1}\gamma_{\lvert r-s\rvert} =\frac{1}{n}\left(\gamma_0+2\sum_{k=1}^{n-1}\Bigl(1-\frac kn\Bigr)\gamma_k\right) ,

counting the $n-k$ pairs at lag $k$ for each $k\ge1$ . If $\sum_k\lvert\gamma_k\rvert<\infty$ , which geometric mixing guarantees, dominated convergence lets $n\to\infty$ inside and the bracket tends to $\gamma_0+2\sum_{k\ge1}\gamma_k=\gamma_0\tau$ . $\square$

Remark(Where the hypotheses are used, and the traps)Original

Finite variance is the hypothesis that matters. The strong law needs only $\mathrm{E}\lvert f\rvert<\infty$ , so the estimator converges; but if $\sigma^{2}=\infty$ the error bar is meaningless, the central limit theorem does not apply, and the sample standard deviation $s_n$ wanders without settling — while the running mean, deceptively, still looks as though it is converging. This is the characteristic failure of an importance density with tails lighter than $\lvert g\rvert q$ : then $w g$ has heavy tails, $\int g^{2}q^{2}/p$ diverges, and the estimator is worse than the crude one it was meant to improve, yet it reports a small error bar until a rare large $w$ arrives. The safe habit is to choose $p$ with tails at least as heavy as $q$ and to monitor $\max_i w(Y_i)/\sum_i w(Y_i)$ .

Unbiasedness is not automatic for a ratio. When $q$ is known only up to a constant, one uses the self-normalised estimator $\sum_i w_ig(Y_i)/\sum_i w_i$ . It is consistent but biased at order $1/n$ ; the bias is harmless beside the $n^{-1/2}$ noise, but it means two such estimates may not be averaged naively.

Hit-or-miss is always worse than crude. With $0\le f\le1$ , the hit-or-miss indicator has variance $\theta(1-\theta)=\mathrm{E}[f]-\theta^{2}\ \ge\ \mathrm{E}[f^{2}]-\theta^{2}=\sigma^{2}$ , since $f^{2}\le f$ on $[0,1]$ . It also spends two random numbers per sample instead of one. Its merits are elsewhere: it needs no bound on $f$ beyond the box, and it generalises immediately to a region of complicated shape, which is Riley, Hobson and Bence's octahedron example.

The crossover in part 2 is asymptotic and the constants matter. At $d=10$ with a smooth integrand, a sparse grid or a quasi-random (low-discrepancy) sequence will usually beat plain Monte Carlo in practice — but not for the reason the rates suggest, and the rate bound must not be offered as the explanation. The Koksma–Hlawka bound for a low-discrepancy sequence is $O\bigl(n^{-1}(\log n)^{d}\bigr)$ , and at $d=10$ that bound is worse than $n^{-1/2}$ throughout every range anyone computes in. Comparing the two, $(\log n)^{d}n^{-1/2}$ equals $1$ when $d\log\log n=\tfrac12\log n$ ; at $d=10$ this gives $\log n=89.995$ , that is $n=1.2\times10^{39}$ , and below that the quasi-random bound is the larger one — by a factor $4.4\times10^{7}$ at $n=10^{4}$ , worst of all ( $4.6\times10^{8}$ ) near $n=\mathrm{e}^{2d}=4.9\times10^{8}$ , and still $1.1$ at $n=10^{39}$ . What actually wins is that the bound is grossly pessimistic for an integrand of low effective dimension, where the observed error behaves like $n^{-1}$ with no power of $\log n$ in evidence. The rate bound therefore settles the question in neither direction at $d=10$ . The honest statement of the case for Monte Carlo is not the rate alone but the rate together with the fact that the method needs nothing from the integrand — no smoothness, no product structure, no description of the domain beyond the ability to test membership.

Burn-in and ergodicity are different questions. Part 4 says the time average converges for any $M$ , including $M=0$ : discarding a burn-in is about reducing bias at finite $n$ , not about correctness. The hard practical question is $\tau$ , not $M$ , and Grimmett and Stirzaker are explicit that deciding when a chain is near equilibrium is notoriously difficult. What makes $\tau$ large is a chain that moves slowly through the state space; the acceptance rule is the usual culprit, since a proposal that is too bold is almost always rejected and a proposal that is too timid moves nowhere.

Detailed balance is sufficient, not necessary. Part 4 builds a reversible chain because reversibility makes the construction easy and the spectral bound explicit. A chain can have $\pi$ stationary without being reversible, and such chains are deliberately used to reduce $\tau$ .

Method(Running a Monte Carlo calculation)Original

Write the target as an expectation. Identify $\theta=\mathrm{E}_p[h(Y)]$ for some density $p$ you can sample and some $h$ you can evaluate. For $\int_D g$ over a complicated $D$ inside a box $B$ , take $p$ uniform on $B$ and $h=\lvert B\rvert\,g\,\mathbf{1}_D$ .
Generate the uniforms. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.4.4, p. 1017) give the linear congruential recipe $x_i=ax_{i-1}+c \pmod m$ with $\xi_i=x_i/m$ , where $c$ is odd and $a=4k+1$ , which runs through all $m$ values before repeating; they take $m$ a large power of two. Any modern library generator is better, but fix and record the seed: a Monte Carlo result that cannot be reproduced cannot be debugged.
Transform to the distribution you want. If $F$ is the distribution function of the target and $F^{-1}$ is available, $Y=F^{-1}(\xi)$ has distribution $F$ , because $\Pr(F^{-1}(\xi)\le y)=\Pr(\xi\le F(y))=F(y)$ ; this is the inverse-transform rule of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.4.4, p. 1017), who also give the look-up table for a discrete or tabulated $F$ . When $F^{-1}$ has no closed form, use rejection, a look-up table, or a special-purpose transformation.
Run and accumulate. Keep running sums of $h$ and $h^{2}$ ; never store the samples. After $n$ draws report $\hat\theta_n\pm 1.96\,\frac{s_n}{\sqrt n}, \qquad s_n^{2}=\frac{n}{n-1}\left(\frac{1}{n}\sum h_i^{2}-\hat\theta_n^{2}\right) .$
Decide how many samples. Inverting the error formula, a target root-mean-square error $\epsilon$ needs $n=\sigma^{2}/\epsilon^{2}$ . Estimate $\sigma$ from a short pilot run. Halving the error costs four times the work; one more decimal digit costs a hundred times.
Reduce the variance before increasing $n$ . Factor out what you can integrate: if $g\approx \tilde g$ with $\int\tilde g$ known, estimate $\int(g-\tilde g)$ instead (a control variate). Or sample from $p\propto\tilde g$ (importance sampling). Either way the figure of merit is variance $\times$ cost per sample, since doubling the cost per sample to halve the variance is a wash.
If you cannot sample the target, build a chain. Choose a symmetric proposal $H$ , use $a_{ij}=\min\{1,\pi_j/\pi_i\}$ , and tune the proposal so that a sensible fraction of moves is accepted. Discard a burn-in, then estimate $\tau$ from the empirical autocorrelations and quote the error as $\sqrt{\gamma_0\tau/n}$ , not $\sqrt{\gamma_0/n}$ — the second is the commonest way to publish an error bar that is too small by a factor of several.

Example(An error function, an importance density, and a three-state Metropolis chain)Original

(a) Crude Monte Carlo and its exact variance. Estimate

\theta=\int_0^1 \mathrm{e}^{-x^{2}}\,\mathrm{d}x=\frac{\sqrt\pi}{2}\operatorname{erf}(1)=0.746824132812

(Error function and incomplete Gamma function). The crude estimator is $\hat\theta_n=n^{-1}\sum \mathrm{e}^{-\xi_i^{2}}$ . Its variance constant is exactly computable: writing $J(a)=\int_0^1\mathrm{e}^{-at^{2}}\mathrm{d}t$ ,

\sigma^{2}=J(2)-J(1)^{2}=0.598144006661-0.557746285351=0.040397721310, \qquad \sigma=0.200991844 .

So an accuracy of $10^{-3}$ needs $n=\sigma^{2}/10^{-6}=40398$ samples.

(b) An importance density, chosen and used. The integrand falls from $1$ to $\mathrm{e}^{-1}=0.3679$ , so sample more near $0$ . Take the linear density

p(y)=\tfrac23(2-y)\ \text{ on }[0,1], \qquad P(y)=\int_0^y p=\frac{4y-y^{2}}{3} ,

which is positive throughout $[0,1]$ — so the ratio $w=1/p$ is bounded — and invertible in closed form: solving $P(y)=\xi$ , that is $y^{2}-4y+3\xi=0$ , and taking the root in $[0,1]$ gives

y=2-\sqrt{4-3\xi} .

The weight is $w(y)g(y)=\mathrm{e}^{-y^{2}}\big/\bigl[\tfrac23(2-y)\bigr]=\dfrac{3\,\mathrm{e}^{-y^{2}}}{2(2-y)}$ , and by part 3 its variance is

\sigma^{2}_{\mathrm{IS}}=\int_0^1\frac{\mathrm{e}^{-2y^{2}}}{p(y)}\,\mathrm{d}y-\theta^{2} =0.562409079663-0.557746285351=0.004662794312, \qquad \sigma_{\mathrm{IS}}=0.068284656 .

The variance has dropped by the factor $0.040397721/0.004662794=8.664$ : the same accuracy now needs $n=4663$ samples instead of $40398$ . (The weight ranges only over $[0.5518,\,0.8064]$ on $[0,1]$ : its smallest value there is $wg(1)=0.551819$ and its largest is $wg(1-1/\sqrt2)=0.806444$ , at the stationary point $y=1-1/\sqrt2=0.292893$ of $\log wg$ . That near-constancy is the reason for the gain, since part 3 says the variance measures how far $wg$ is from constant.)

(c) The two estimators on the same ten numbers. Generate $\xi_i$ with the congruential recipe of step 2, taking $m=2^{31}$ , $a=1103515245$ , $c=12345$ and $x_0=1$ :

\xi=(0.5138701,\ 0.1757413,\ 0.3086515,\ 0.5345339,\ 0.9476279,\ 0.1717363,\ 0.7022312,\ 0.2264307,\ 0.4947734,\ 0.1247203).

Crude: averaging $\mathrm{e}^{-\xi_i^{2}}$ gives $\hat\theta_{10}=0.8104593$ , an error of $+0.063635$ . The predicted standard deviation is $\sigma/\sqrt{10}=0.0635592$ — the error landed at $1.001$ standard deviations.

Importance sampling, on the same ten uniforms: $y_i=2-\sqrt{4-3\xi_i}$ gives

y=(0.4320747,\ 0.1364614,\ 0.2467044,\ 0.4519695,\ 0.9243066,\ 0.1332405,\ 0.6240253,\ 0.1777190,\ 0.4139106,\ 0.0958364),

and averaging $3\mathrm{e}^{-y_i^{2}}/[2(2-y_i)]$ gives $\hat\theta^{\,\mathrm{IS}}_{10}=0.7675041$ , an error of $+0.020680$ against a predicted standard deviation $\sigma_{\mathrm{IS}}/\sqrt{10}=0.0215935$ — $0.958$ standard deviations. Ten samples is far too few to demonstrate anything statistically, and that is the point of quoting the predicted deviations: the two errors are where they should be, and the two predicted deviations differ by the factor $\sqrt{8.664}=2.94$ that part 3 forecast.

(d) A Metropolis chain small enough to do exactly. Take $\Theta=\{1,2,3\}$ with target $\pi\propto(1,2,3)$ , so $\pi=(\tfrac16,\tfrac13,\tfrac12)$ , and the symmetric proposal "pick one of the other two states at random", $h_{ij}=\tfrac12$ for $i\ne j$ . The Metropolis rule $a_{ij}=\min\{1,\pi_j/\pi_i\}$ gives, for instance, $a_{12}=\min\{1,2\}=1$ and $a_{21}=\min\{1,\tfrac12\}=\tfrac12$ , hence $p_{12}=\tfrac12$ and $p_{21}=\tfrac14$ . Carrying this out for every pair:

P=\begin{pmatrix} 0 & \tfrac12 & \tfrac12\\[2pt] \tfrac14 & \tfrac14 & \tfrac12\\[2pt] \tfrac16 & \tfrac13 & \tfrac12 \end{pmatrix} .

Checks. Rows sum to $1$ . Detailed balance: $\pi_1p_{12}=\tfrac16\cdot\tfrac12=\tfrac1{12}$ and $\pi_2p_{21}=\tfrac13\cdot\tfrac14=\tfrac1{12}$ ; $\pi_1p_{13}=\tfrac1{12}=\tfrac12\cdot\tfrac16=\pi_3p_{31}$ ; $\pi_2p_{23}=\tfrac13\cdot\tfrac12=\tfrac16=\tfrac12\cdot\tfrac13=\pi_3p_{32}$ . So $\pi P=\pi$ , as part 4 promises.

Rate. The characteristic polynomial of $P$ factorises as $\lambda(\lambda-1)(\lambda+\tfrac14)$ , so the spectrum is $\{1,-\tfrac14,0\}$ and the second largest eigenvalue modulus is $\lvert\lambda_2\rvert=\tfrac14$ . Starting from state $1$ , the exact $n$ -step row is

p_{1\cdot}(n)=\Bigl(\tfrac16+\tfrac{2}{3}\bigl(-\tfrac14\bigr)^{n},\ \ \tfrac13-\tfrac{2}{3}\bigl(-\tfrac14\bigr)^{n},\ \ \tfrac12\Bigr),\qquad n\ge1

— the third component is exactly right after a single step, because the third column of $P$ is constant. Hence, in the normalisation of Convergence of Markov chains,

d_{\mathrm{TV}}\bigl(p_{1\cdot}(n),\pi\bigr)=\tfrac43\,4^{-n} = \tfrac13,\ \tfrac1{12},\ \tfrac1{48},\ \tfrac1{192},\dots\ \text{ for } n=1,2,3,4 .

The reversible spectral bound of that unit gives, with $N=3$ , $\lvert\lambda_2\rvert=\tfrac14$ and the $\pi$ -orthonormal right eigenvectors $v_1=(1,1,1)$ , $v_2=(-2,1,0)$ , $v_3=(1,1,-1)$ (each already of unit norm for $\langle x,y\rangle=\sum\pi_ix_iy_i$ ), $\sup_r\lvert v_r(1)\rvert=2$ and so $d_{\mathrm{TV}}\le 6\cdot4^{-n}$ . True value $\tfrac43 4^{-n}$ : the bound has the rate exactly right and the prefactor too large by $4.5$ .

Efficiency. Estimate $\theta=\mathrm{E}_\pi[g]$ with $g(i)=i$ , so $\theta=\tfrac16+\tfrac23+\tfrac32=\tfrac73$ . Then $\gamma_0=\operatorname{Var}_\pi g=6-\tfrac{49}{9}=\tfrac59$ , and expanding $g$ in the eigenbasis gives $\gamma_k=\tfrac19(-\tfrac14)^{k}$ for $k\ge1$ (the $\lambda=0$ mode contributes only at $k=0$ ). Hence

\tau=1+2\sum_{k\ge1}\frac{\gamma_k}{\gamma_0} =1+2\cdot\frac95\cdot\frac19\cdot\frac{-\tfrac14}{1+\tfrac14} =1-\frac{2}{25}=\frac{23}{25}=0.92 ,

so the predicted variance of the time average is $\gamma_0\tau/n=\tfrac{23}{45n}=0.511111/n$ . Computing the variance exactly from the finite- $n$ formula of the proof gives $0.0514667$ at $n=10$ and $0.00511467$ at $n=100$ , against the predictions $0.0511111$ and $0.00511111$ . Note $\tau<1$ : successive states are negatively correlated, because a Metropolis move from state $1$ never stays put, so this chain is very slightly more efficient than independent sampling. That is unusual, and it happens here only because the chain mixes in essentially one step.

Application(Where the method earns its keep)Original

Statistical mechanics. The canonical average $\langle A\rangle=Z^{-1}\sum_{\vartheta}A(\vartheta)\mathrm{e}^{-\beta E(\vartheta)}$ over the $2^{L^{d}}$ configurations of a lattice model cannot be enumerated or sampled directly, and $Z$ is not known. The Metropolis rule needs neither: a single spin flip changes the energy by a local amount $\Delta E$ and is accepted with probability $\min\{1,\mathrm{e}^{-\beta\Delta E}\}$ . This is the original application — the algorithm was published for exactly this problem in the early 1950s — and it remains the standard tool for lattice gauge theory, for polymer and protein conformations, and for phase diagrams. Near a critical point $\tau$ diverges with system size, which is the critical slowing down discussed in Convergence of Markov chains.

Particle transport and radiation. Neutron transport, radiotherapy dose planning and detector response are all high-dimensional integrals over particle histories, with geometry too complicated for a mesh. Each simulated history is one sample, variance reduction is done by splitting and Russian roulette in the regions that matter, and the reported dose carries a statistical error bar computed exactly as in step 4. Grimmett and Stirzaker record that the name and the method date from the Los Alamos work of Metropolis, von Neumann and Ulam in the 1940s, where this was the problem.

Financial valuation. The price of a path-dependent option is an expectation under a risk-neutral measure over paths of the stochastic differential equation of Stochastic differential equations; the dimension is the number of time steps, often hundreds, so no quadrature rule is possible. Importance sampling is used in precisely the form of part 3 to price deep out-of-the-money options, where the payoff is non-zero on an event of probability $10^{-4}$ , so that $10^{8}$ crude paths would produce only $10^{4}$ non-zero ones and the other $99.99\%$ of the work would be wasted.

Bayesian inference. The posterior $\pi(\vartheta\mid x)\propto f(x\mid\vartheta)\pi(\vartheta)$ is known up to its normalising constant by construction, which is the situation part 4 was built for; Markov chain Monte Carlo is what made Bayesian analysis of large structured models practical, and Grimmett and Stirzaker open their treatment with this example. The quantity a practitioner must report alongside a posterior mean is the effective sample size $n/\tau$ , which is the content of part 4's variance formula.

And where it does not. For a smooth one-dimensional integral, Monte Carlo is absurd: Gauss–Legendre quadrature (Numerical quadrature) gets twenty digits from twenty points, while Monte Carlo would need some $10^{39}$ samples (step 5 with $\epsilon=7.5\times10^{-21}$ gives $n=\sigma^{2}/\epsilon^{2}=7.2\times10^{38}$ ). The method is for problems whose dimension or whose geometry defeats a grid, and for nothing else.