Numerical analysis basics

Errors, roots, interpolation, quadrature and linear algebra.

In this chapter 6

Tick a unit to add it to your PDF basket.

Basic numerical methods

Basic numerical methods

Definition

Floating-point arithmetic and error

Adapted

Add to PDF

Defines the floating-point numbers, machine epsilon and the unit roundoff, proves the relative-error bound that every rounded operation obeys, and separates the conditioning of a problem (which no algorithm can repair) from the stability of an algorithm (which a better algorithm can).

NoteOriginal

What it means. A computer does not do arithmetic in $\mathbb{R}$ . It works in a finite set of numbers, and every operation lands on the nearest member of that set. The single fact worth carrying away is quantitative: in IEEE double precision each operation commits a relative error of at most $1.11 \times 10^{-16}$ , and that is all the arithmetic itself does to you. Everything else — the catastrophic loss of digits, the answer that is wrong in the third place — comes from how those tiny errors are propagated by the formula you chose, or from the problem itself being hypersensitive to its data. Keeping those two apart is the whole content of this unit.

Definition(Floating-point numbers, machine epsilon, unit roundoff)Original

Fix a base $\beta \ge 2$ , a precision $t \ge 1$ and exponent bounds $e_{\min} \le e_{\max}$ . The set of normalised floating-point numbers is

\mathbb{F} = \{0\} \cup \Bigl\{ \pm m \,\beta^{\,e-t+1} \ : \ \beta^{\,t-1} \le m \le \beta^{\,t}-1,\ m \in \mathbb{Z},\ e_{\min} \le e \le e_{\max} \Bigr\} .

Writing $m$ in base $\beta$ displays the familiar form $\pm (d_0 . d_1 d_2 \cdots d_{t-1})_\beta \times \beta^{e}$ with the leading digit $d_0 \neq 0$ ; the condition $m \ge \beta^{t-1}$ is exactly that normalisation. A non-zero $x \in \mathbb{F}$ with exponent $e$ satisfies $\beta^{e} \le |x| < \beta^{e+1}$ , and the spacing between neighbouring elements of $\mathbb{F}$ in that range — one unit in the last place, $\operatorname{ulp}$ — is $\beta^{\,e-t+1}$ . The spacing is therefore not uniform: it is proportional to the magnitude, which is the defining feature of a floating-point system and the reason relative error, not absolute error, is the natural currency here.

Two numbers measure the fineness of $\mathbb{F}$ near $1$ :

the machine epsilon $\varepsilon_{M} = \beta^{\,1-t}$ , the gap between $1$ and the next larger element of $\mathbb{F}$ ;
the unit roundoff $u = \tfrac{1}{2}\beta^{\,1-t} = \tfrac{1}{2}\varepsilon_M$ , half that gap.

Usage differs between books and between software: some call $u$ the machine epsilon, and some programming environments report $\varepsilon_M$ under that name. This unit uses the two symbols above throughout, and every statement below is in terms of $u$ .

Let $\mathrm{fl} : \mathbb{R} \to \mathbb{F}$ be round to nearest: $\mathrm{fl}(x)$ is the element of $\mathbb{F}$ closest to $x$ , ties broken by any fixed rule. Say $x$ is in range if $\beta^{\,e_{\min}} \le |x| \le (\beta - \beta^{\,1-t})\beta^{\,e_{\max}}$ , the largest element of $\mathbb{F}$ being $(\beta^t - 1)\beta^{\,e_{\max}-t+1} = (\beta - \beta^{1-t})\beta^{e_{\max}}$ . Below the lower bound the number underflows, above the upper bound it overflows.

The IEEE-754 binary64 format ("double precision") is $\beta = 2$ , $t = 53$ , $e_{\min} = -1022$ , $e_{\max} = 1023$ , so

\varepsilon_M = 2^{-52} = 2.220446049250313 \times 10^{-16}, \qquad u = 2^{-53} = 1.1102230246251565 \times 10^{-16},

the smallest positive normalised number is $2^{-1022} = 2.2250738585072014 \times 10^{-308}$ and the largest finite number is $(2-2^{-52})2^{1023} = 1.7976931348623157 \times 10^{308}$ . Since $53 \log_{10} 2 = 15.95$ , a double carries just under sixteen decimal digits.

Definition(The rounding error of one operation)Original

Let $x \in \mathbb{R}$ be in range. Then

\mathrm{fl}(x) = x\,(1+\delta) \qquad \text{for some } \delta \text{ with } |\delta| \le u .

If moreover $a, b \in \mathbb{F}$ , $\ast$ is one of $+,-,\times,\div$ (with $b \neq 0$ for $\div$ ), and the exact value $a \ast b$ is in range, then the IEEE-754 requirement that an operation return the correctly rounded exact result gives

\mathrm{fl}(a \ast b) = (a \ast b)(1+\delta), \qquad |\delta| \le u ,

and likewise $\mathrm{fl}(\sqrt{a}) = \sqrt{a}\,(1+\delta)$ with $|\delta| \le u$ . This pair of statements is called the standard model of floating-point arithmetic.

Idea of the proofOriginal

Everything follows from the spacing. The numbers of $\mathbb{F}$ between two consecutive powers $\beta^{e}$ and $\beta^{e+1}$ sit on a uniform grid of mesh $\beta^{e-t+1}$ . Rounding to the nearest grid point costs at most half a mesh in absolute terms; dividing that by $|x| \ge \beta^{e}$ turns it into a relative bound in which the exponent $e$ cancels. The cancellation of $e$ is the whole point: the bound is the same everywhere in range, which is why floating point is useful across three hundred orders of magnitude. The second half of the statement is then not a separate fact but a definition honoured: the hardware is required to compute $a \ast b$ exactly and then round, so the only error is the rounding already bounded.

ProofOriginal

Let $x$ be in range; by symmetry take $x > 0$ , and let $e$ be the integer with $\beta^{e} \le x < \beta^{e+1}$ . Being in range forces $e_{\min} \le e \le e_{\max}$ .

The grid. The elements of $\mathbb{F}$ with exponent $e$ are the numbers $m\beta^{\,e-t+1}$ for integers $m$ with $\beta^{t-1} \le m \le \beta^{t}-1$ . They run from $\beta^{e}$ up to $\beta^{e+1} - \beta^{\,e-t+1}$ in steps of $\beta^{\,e-t+1}$ . Two cases complete the grid. If $e < e_{\max}$ , the next element of $\mathbb{F}$ is $\beta^{e+1}$ (exponent $e+1$ , $m = \beta^{t-1}$ ), one further step along, so the whole closed interval $[\beta^{e}, \beta^{e+1}]$ is covered by grid points of spacing $\beta^{\,e-t+1}$ . If $e = e_{\max}$ , then $\beta^{e+1} \notin \mathbb{F}$ , but the largest element of $\mathbb{F}$ is $\beta^{e+1}-\beta^{\,e-t+1}$ and $x$ being in range means $x$ does not exceed it; so $x$ again lies in the gridded part.
Rounding cost. By step 1, $x$ lies between two consecutive grid points of spacing $\beta^{\,e-t+1}$ , hence within $\tfrac12\beta^{\,e-t+1}$ of the nearer one. Since $\mathrm{fl}(x)$ is the nearest element of all of $\mathbb{F}$ , it is at least as close:

|\mathrm{fl}(x) - x| \le \tfrac12 \beta^{\,e-t+1} .

Relative form. Divide by $x$ and use $x \ge \beta^{e}$ :

\frac{|\mathrm{fl}(x) - x|}{x} \le \frac{\tfrac12\beta^{\,e-t+1}}{\beta^{e}} = \tfrac12 \beta^{\,1-t} = u .

Setting $\delta = (\mathrm{fl}(x)-x)/x$ gives $\mathrm{fl}(x) = x(1+\delta)$ with $|\delta| \le u$ , which is the first claim.

Operations. IEEE-754 specifies that $+,-,\times,\div$ and $\sqrt{\ }$ return the correctly rounded value of the exact mathematical result. So the computed value of $a \ast b$ is $\mathrm{fl}(a \ast b)$ with $a \ast b$ the exact real number, and step 3 applied to that real number gives the stated bound. $\square$

Remark(Where the hypotheses are used)Original

In range. The bound $x \ge \beta^{e}$ in step 3 is the only place normalisation enters, and it is where the proof fails on underflow. Below $\beta^{e_{\min}}$ the grid stops shrinking: IEEE-754 fills that gap with subnormal numbers, multiples of $\beta^{\,e_{\min}-t+1}$ , so only the absolute bound $|\mathrm{fl}(x)-x| \le \tfrac12\beta^{\,e_{\min}-t+1}$ survives. For $x = 2^{-1080}$ in binary64 the nearest representable numbers are $0$ and $2^{-1074}$ , and the relative error of rounding is of order $1$ . Overflow is worse still: there is no bound at all, the result being $\infty$ .
Round to nearest. Step 2 uses it. With truncation toward zero the error is at most a whole step rather than half, and the bound doubles to $|\delta| \le \varepsilon_M = 2u$ ; the error also acquires a systematic sign, so it accumulates linearly over a long computation instead of partly cancelling.
One operation at a time. The model governs each operation separately. It says nothing about a sequence of them, and no amount of staring at it will tell you the error of a formula; that requires propagating the $\delta$ 's through the formula, which is what the condition number below does.
$a, b$ already in $\mathbb{F}$ . The second half assumes the operands are representable. If they are themselves rounded data the two errors compose, giving $(1+\delta_1)(1+\delta_2) = 1 + \delta_1 + \delta_2 + O(u^2)$ : relative errors add to first order.

Remark(Subtraction is innocent: cancellation reveals error, it does not create it)Original

A subtraction of nearly equal numbers is often blamed for a loss of accuracy. The blame is misplaced, and sharply so.

Write $\mathbb{F}^{\ast}$ for $\mathbb{F}$ together with the subnormals, the integer multiples of $\beta^{\,e_{\min}-t+1}$ of modulus below $\beta^{\,e_{\min}}$ , which IEEE-754 provides.

Lemma (Sterbenz). Let $\beta = 2$ and let $a, b \in \mathbb{F}$ be positive with $\tfrac12 b \le a \le 2b$ . Then $a - b \in \mathbb{F}^{\ast}$ , so $\mathrm{fl}(a-b) = a-b$ exactly.

Proof. Take $a \ge b$ (otherwise swap and negate, which changes nothing). Let $e$ be the exponent of $b$ , so $2^{e} \le b < 2^{e+1}$ and $b$ is an integer multiple of $2^{\,e-t+1}$ . Since $a \ge b$ the exponent of $a$ is some $e_a \ge e$ , so $a$ is a multiple of $2^{\,e_a-t+1}$ , hence of $2^{\,e-t+1}$ ; so is the difference. The hypothesis $a \le 2b$ gives $0 \le a-b \le b < 2^{e+1}$ . Set $v = a-b$ and suppose $v > 0$ . Let $e'$ be the integer with $2^{e'} \le v < 2^{e'+1}$ ; then $e' \le e$ , and $v$ , being a multiple of $2^{\,e-t+1} = 2^{\,e-e'}\cdot 2^{\,e'-t+1}$ , is a multiple of $2^{\,e'-t+1}$ . If $e' \ge e_{\min}$ , then $v = m\,2^{\,e'-t+1}$ with $2^{t-1} \le m < 2^{t}$ , a normalised element of $\mathbb{F}$ . If $e' < e_{\min}$ , then $v < 2^{\,e_{\min}}$ while $v$ is a multiple of $2^{\,e-t+1}$ and $e \ge e_{\min}$ , so $v$ is a multiple of $2^{\,e_{\min}-t+1}$ : a subnormal. Either way $v \in \mathbb{F}^{\ast}$ , and $v = 0$ is in $\mathbb{F}$ trivially. $\square$

So when two close numbers are subtracted, the subtraction itself is error-free. What happens is that the errors already present in the operands, which were small relative to the operands, are now large relative to the much smaller difference. The worked example below shows this happening digit by digit: the square root is rounded, the subtraction is exact, and the answer is wrong in the third decimal place.

Definition(Conditioning of a problem, stability of an algorithm)AdaptedAdapted from Strauss (2008), §1.5, the stability requirement for well-posedness and the ill-conditioned matrix with a very small eigenvalue, pp. 26-27

Regard a problem as a map $f : X \to Y$ from data to answer. Its relative condition number at $x$ measures how much a relative perturbation of the data is magnified in the answer:

\kappa_f(x) \;=\; \lim_{\epsilon \to 0} \ \sup_{0 < \|\Delta x\| \le \epsilon \|x\|} \ \frac{\|f(x+\Delta x) - f(x)\| \big/ \|f(x)\|}{\|\Delta x\| \big/ \|x\|} .

For a differentiable scalar $f$ with $f(x) \neq 0$ this is $\kappa_f(x) = |x f'(x)| / |f(x)|$ . A problem is ill-conditioned when $\kappa_f$ is large; Strauss (Strauss (2008), §1.5, p. 27) introduces the word for the matrix case, where the data are the right-hand side of $Au = b$ and a very small eigenvalue of $A$ magnifies a perturbation of $b$ (the matrix condition number $\|A\|\,\|A^{-1}\|$ is treated in Matrix norms and condition numbers). Conditioning is a property of the problem: it is there before any algorithm is chosen, and no algorithm removes it, because the data themselves are only known to within a rounding.

An algorithm is a second map $\tilde f : X \to Y$ , the one the machine actually computes. It is backward stable if for every $x$ there is a nearby $\tilde x$ with

\tilde f(x) = f(\tilde x), \qquad \frac{\|\tilde x - x\|}{\|x\|} = O(u) :

the computed answer is the exact answer to a slightly different question. Backward stability is a property of the algorithm, and two algorithms for the same problem can differ in it completely. The two notions combine in one inequality, proved by the definition of $\kappa_f$ : if $\tilde f(x) = f(\tilde x)$ with relative backward error $\eta$ , then

\frac{\|\tilde f(x) - f(x)\|}{\|f(x)\|} \ \le \ \kappa_f(x)\,\eta + o(\eta) .

Forward error $\lesssim$ condition number $\times$ backward error. A large forward error therefore has exactly two possible causes, and the inequality tells you which: either $\kappa_f$ is large, in which case the problem is to blame and a different algorithm will not help, or the algorithm is unstable and should be replaced.

Remark(Rounding error and truncation error are different things)AdaptedAdapted from Riley, Hobson & Bence (2006), Ch. 27, introduction, pp. 984-985

Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Ch. 27, pp. 984-985) separate two sources of numerical error, and the distinction is used throughout this part of the handbook. Rounding error comes from the finite precision with which a value is stored, and is what this unit is about; it is governed by $u$ and does not improve when you refine a discretisation. Truncation error (also called discretisation error) comes from the approximation built into the numerical procedure itself — replacing a derivative by a difference quotient (Finite-difference approximations), an integral by a finite sum (Numerical quadrature), a function by a polynomial (Polynomial and spline interpolation) — and it shrinks as the step $h$ shrinks. The two run in opposite directions as $h \to 0$ : truncation error falls like a power of $h$ while rounding error, divided by $h$ or $h^2$ in a difference quotient, grows. The optimal step is where they cross, and it is never zero.

Method(Diagnosing and curing a loss of accuracy)Original

Estimate the condition number of the problem at your data, from $\kappa_f = |xf'(x)|/|f(x)|$ for a scalar problem or from $\|A\|\|A^{-1}\|$ for a linear system. If $\kappa_f \approx 10^{k}$ , expect to lose $k$ of your sixteen digits no matter what you do. Stop here if $k$ is large: the fix is a reformulation of the problem, better data, or extended precision — not a better algorithm.
If the problem is well conditioned but the answer is bad, the algorithm is unstable. Look for the step at which a quantity much smaller than its inputs is formed. That step is where previously harmless relative errors become visible.
Remove the cancellation algebraically, using an identity that computes the small quantity directly instead of as a difference:

\sqrt{x+1}-\sqrt{x} = \frac{1}{\sqrt{x+1}+\sqrt{x}}, \qquad 1 - \cos x = 2\sin^{2}\tfrac{x}{2}, \qquad \ln(1+x) \ \text{by its series for small } x .

For the roots of $ax^2+bx+c$ , compute the root of larger modulus by the usual formula with the sign of $b$ chosen to add, then obtain the other from the product of the roots, $x_1 x_2 = c/a$ . 4. Reorder the arithmetic where you cannot remove the cancellation: sum a series from the small terms upward, and pivot before eliminating. Pivoting is the standard instance: Gaussian elimination with interchange (Riley, Hobson & Bence (2006), §27.3.1, pp. 995-996) reorders the equations before each elimination so that the largest remaining coefficient occupies the diagonal position, which bounds the multipliers by $1$ and stops the rounding errors growing. Same problem, same arithmetic count, different stability (Solving linear systems numerically). 5. Re-run in higher precision and compare. Agreement to $k$ digits between double and quadruple precision is evidence, not proof, that $k$ digits are right; disagreement is proof that they are not.

Example(A well-conditioned problem that a natural formula gets wrong in the third decimal place)Original

Take the quadratic

x^{2} - 2bx + 1 = 0, \qquad b = 10^{7},

whose roots are $r_{\pm} = b \pm \sqrt{b^{2}-1}$ with $r_{+}r_{-} = 1$ . We want the small root $r_{-}$ . All of what follows is IEEE binary64 arithmetic.

Step 1: the problem is perfectly conditioned. Treat $b$ as the datum. From $r_{-} = b - \sqrt{b^{2}-1}$ ,

\frac{dr_{-}}{db} = 1 - \frac{b}{\sqrt{b^{2}-1}} = -\frac{r_{-}}{\sqrt{b^{2}-1}}, \qquad\text{so}\qquad \kappa(b) = \left|\frac{b}{r_{-}}\frac{dr_{-}}{db}\right| = \frac{b}{\sqrt{b^{2}-1}} = 1.000000000000005 .

A relative perturbation of $b$ changes $r_{-}$ by the same relative amount. Nothing about this problem is delicate.

Step 2: the natural formula, step by step. $b = 10^{7}$ is in $\mathbb{F}$ . So is $b^{2} = 10^{14}$ , since $10^{14} = 2^{14}\cdot 5^{14}$ and $5^{14} = 6103515625 < 2^{53}$ ; so is $b^{2}-1 = 99999999999999 < 2^{53}$ , and both operations are exact. The square root is not exact: its true value is

\sqrt{99999999999999} = 9999999.99999995000000000125\ldots,

and the nearest double is

\mathrm{fl}(\sqrt{b^2-1}) = 9999999.99999994970858097076416015625 ,

low by $2.914 \times 10^{-10}$ . That is a relative error of $2.91 \times 10^{-17}$ , comfortably inside the guaranteed $u = 1.11\times10^{-16}$ : the square root did its job.

Step 3: the subtraction, which is exact. The operands satisfy $\tfrac12 b \le \mathrm{fl}(\sqrt{b^2-1}) \le 2b$ , so by the lemma above the subtraction commits no error at all and returns the exact difference of the two stored numbers,

\widetilde r_{-} = 10^{7} - 9999999.99999994970858097076416015625 = 5.029141902923583984375\times10^{-8},

which rounds to the double $5.029141902923584\times10^{-8}$ .

Step 4: compare. The true root, to twenty digits, is

r_{-} = \frac{1}{10^{7}+\sqrt{10^{14}-1}} = 5.0000000000000125\times10^{-8} .

So the computed value is wrong by $2.914\times10^{-10}$ — precisely the error made in step 2, carried through unchanged — and its relative error is $5.83 \times 10^{-3}$ . Two correct digits out of sixteen, from a problem with $\kappa = 1$ and an arithmetic that never erred by more than $3\times10^{-17}$ relatively.

Step 5: the stable formula. Using $r_{-} = 1/(b+\sqrt{b^{2}-1})$ there is no subtraction of near-equals. The same square root, now added to $b$ , gives the double

\widetilde r_{-}^{\,\text{stable}} = 5.000000000000013\times10^{-8},

with relative error $1.02\times10^{-16}$ , that is $0.92\,u$ : full accuracy. The two formulas are the same function of $b$ on paper and differ by a factor of $5.7\times10^{13}$ in accuracy. That factor is the difference between an unstable algorithm and a stable one, and the condition number $\kappa = 1$ is what tells us the blame lies there and not with the question.

Application(Variance, and why one-pass formulas are banned in production code)Original

Reducing a long record — a time series of temperatures, a run of detector counts, a Monte Carlo sample (Monte Carlo methods) — to a mean and a variance is the most-performed calculation in the applied sciences, and the textbook identity

s^{2} = \frac{1}{n-1}\Bigl(\sum_i x_i^{2} - \frac{1}{n}\Bigl(\sum_i x_i\Bigr)^{2}\Bigr)

is attractive because it needs one pass and two accumulators. It is also a subtraction of two nearly equal large numbers whenever the data have a large mean compared with their spread, which is the normal situation for an instrument reading around a baseline.

Take the five values $x_i = 10^{8}+i$ , $i = 1,\dots,5$ , whose mean is exactly $10^{8}+3$ and whose sample variance is exactly $\tfrac{1}{4}(4+1+0+1+4) = 2.5$ . In binary64: $\sum x_i = 500000015$ is formed exactly, but $\sum x_i^{2}$ has the exact value $50000003000000055$ , which is not representable — $5\times10^{16}$ lies in $[2^{55},2^{56})$ , where the spacing of $\mathbb{F}$ is $2^{\,55-52} = 8$ — and the machine returns the correctly rounded $50000003000000056$ . On the other side, $(\sum x_i)^2 = 250000015000000225$ lies in $[2^{57},2^{58})$ , where the spacing is $32$ , and rounds to $250000015000000224$ ; dividing by $n = 5$ gives $50000003000000044.8$ , which rounds to $50000003000000048$ . The difference of the two stored numbers is $8$ , and the formula returns

s^{2}_{\text{naive}} = \frac{8}{4} = 2.0 \qquad\text{against the true } 2.5 ,

an error of $20\%$ with data that are integers and a calculation with four correct inputs. The two-pass form $s^{2} = \frac{1}{n-1}\sum_i (x_i - \bar x)^2$ , which subtracts before squaring, returns $2.5$ exactly.

The same arithmetic lies behind three familiar rules of practice: subtract a reference value from a geophysical record before computing its statistics; compute differences of total energies in molecular dynamics by accumulating the difference rather than differencing the totals; and, in least-squares fitting, avoid forming the normal-equation matrix at all. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 5.1.3, Note, p. 267) record what happens if you do: fitting a power series on $[0,1]$ produces the Hilbert matrix $A_{ij} = 1/(i+j+1)$ , whose determinant at $n = 5$ is already $3.7\times10^{-12}$ and whose system they describe as becoming ill-conditioned and unstable. That is case 1 of the method above and not case 2 — the problem is ill-conditioned, so the cure is to change the problem, by orthogonalising the basis (General properties of orthogonal polynomials) or by solving the least-squares problem through a QR or singular-value decomposition (Singular value decomposition) instead of through the normal equations.

Method

Root finding

Adapted

Add to PDF

Bisection, Newton-Raphson and the secant method for a scalar equation, with the error recursions that give them orders 1, 2 and 1.618, the hypotheses each one needs, and the ways Newton's method fails when they are dropped.

Prerequisites: Taylor's theorem

NoteOriginal

What it means. Three algorithms, and the whole of the engineering choice between them sits in one line each. Bisection cannot fail and cannot be fast: it buys one bit per step, and you can say in advance how many steps you need. Newton is fast when it works — the number of correct digits doubles each step — and it works only near the root, needs the derivative, and fails in several instructive ways. The secant method is Newton with the derivative replaced by a difference quotient: slightly slower per step, but it costs one function evaluation instead of two, which makes it the better buy when the derivative is expensive. Production code uses all three: bracket with bisection, finish with the other two.

Definition(Roots, multiplicity, order of convergence)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.2, Eqs. (27.18)-(27.19), p. 993

Let $f$ be real-valued on an interval. A root (or zero) of $f$ is a point $\xi$ with $f(\xi) = 0$ . If $f$ is $m$ times differentiable at $\xi$ and

f(\xi) = f'(\xi) = \cdots = f^{(m-1)}(\xi) = 0, \qquad f^{(m)}(\xi) \neq 0 ,

the root has multiplicity $m$ ; a root of multiplicity $1$ is simple, and simple is exactly the condition $f'(\xi) \neq 0$ .

An iterative method produces $x_0, x_1, x_2, \dots$ ; write $e_n = x_n - \xi$ for the error. Suppose $x_n \to \xi$ with $x_n \neq \xi$ for all $n$ . The method converges with order $p \ge 1$ and asymptotic error constant $C$ if

\frac{|e_{n+1}|}{|e_n|^{p}} \longrightarrow C \in (0,\infty) .

Order $1$ with $C < 1$ is called linear convergence (a fixed number of extra correct digits per step, namely $-\log_{10} C$ of them); order $2$ is quadratic (the number of correct digits doubles). Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.2, Eqs. (27.18)-(27.19), p. 993) give the corresponding criterion for a one-point iteration $x_{n+1} = F(x_n)$ with $F$ smooth and $F(\xi) = \xi$ : the iteration is of order $N$ when $F'(\xi) = \cdots = F^{(N-1)}(\xi) = 0$ . For this to match the definition above exactly one must add, in our own voice, that $F^{(N)}(\xi) \neq 0$ , since otherwise the limiting constant $C$ is zero and the order is higher than $N$ . Either way the criterion follows by expanding $F(\xi + e_n)$ in a Taylor series (Taylor's theorem). They note that the error form of the definition is the more general one, since it does not require the iteration to be a differentiable function of a single previous iterate — and indeed neither bisection nor the secant method is of that form.

Method(The three methods and their rates)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.1.3 Eq. (27.9) p. 990 for (a), §27.1.4 Eqs. (27.10)-(27.11) pp. 990-991 for (b), and §27.2 pp. 992-994 for the order of convergence; Kreyszig (1978), §5.1, Problem 15, pp. 305-306 for the local hypothesis in (b). Part (c) is our own: neither source treats the unbracketed secant method

(a) Bisection. Let $f$ be continuous on $[a,b]$ with $f(a)f(b) < 0$ . Put $[a_1,b_1] = [a,b]$ and, for $n \ge 1$ ,

x_n = \tfrac12 (a_n + b_n), \qquad [a_{n+1}, b_{n+1}] = \begin{cases} [a_n, x_n] & \text{if } f(a_n) f(x_n) \le 0, \\ [x_n, b_n] & \text{otherwise.}\end{cases}

The test must be the non-strict one. If a midpoint lands exactly on a root, so that $f(x_n) = 0$ , the rule then keeps that point as the new right-hand endpoint and the bracket never loses it; with $<$ in its place the rule would throw the root away and the conclusion below would be false (take $f(x) = x - \tfrac12$ on $[0,1]$ , where $x_1 = \tfrac12$ is the root and the strict test sends the bracket off to $1$ ).

Then $a_n \uparrow \xi$ and $b_n \downarrow \xi$ for a limit $\xi$ with $f(\xi) = 0$ , and for every $n$

|x_n - \xi| \ \le \ \frac{b-a}{2^{\,n}} .

The bound is known before the computation starts: reaching a tolerance $\tau$ takes $n = \lceil \log_2((b-a)/\tau) \rceil$ steps and $n$ evaluations of $f$ .

(b) Newton-Raphson. Let $\xi$ be a simple root of $f$ and let $f$ be twice continuously differentiable on a closed interval $I$ containing $\xi$ in its interior. There is a $\delta > 0$ such that for every starting point $x_0$ with $|x_0 - \xi| \le \delta$ the iteration

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

is well defined, stays in $[\xi-\delta, \xi+\delta]$ and converges to $\xi$ , with

\frac{e_{n+1}}{e_n^{2}} \longrightarrow \frac{f''(\xi)}{2 f'(\xi)} .

The order is therefore $2$ whenever $f''(\xi) \neq 0$ , and better than $2$ when $f''(\xi) = 0$ .

(c) Secant method. Under the same hypotheses on $f$ and $\xi$ , there is a $\delta > 0$ such that for distinct $x_0, x_1$ in $[\xi-\delta,\xi+\delta]$ the iteration

x_{n+1} = x_n - f(x_n)\,\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}

is well defined and converges to $\xi$ , with

\frac{e_{n+1}}{e_n e_{n-1}} \longrightarrow \frac{f''(\xi)}{2 f'(\xi)} ,

and the errors satisfy $K|e_n| \le \theta^{\,F_n}$ for constants $K>0$ and $\theta<1$ , where $F_0 = F_1 = 1$ and $F_{n+1} = F_n + F_{n-1}$ are the Fibonacci numbers. If in addition $|e_{n+1}|/|e_n|^{p} \to C \in (0,\infty)$ for some $p$ , then

p = \frac{1+\sqrt5}{2} = 1.6180339887\ldots

and $C = |f''(\xi)/2f'(\xi)|^{\,p-1}$ . Each step costs one new evaluation of $f$ and none of $f'$ . Part (c) is not taken from either source — Riley, Hobson and Bence treat only the bracketed chord method of their §27.1.2, and Kreyszig's problems concern one-point iterations — and the proof below establishes it in full.

Idea of the proofOriginal

All three are one idea applied three times: replace $f$ near the current guess by something whose root you can write down, and take that root as the next guess.

Bisection replaces $f$ by its sign alone. One bit of information per evaluation, and a bracket that is never lost — which is why its error bound is a theorem rather than a hope.

Newton replaces $f$ by its tangent at $x_n$ . Taylor's theorem says the tangent misses $f$ by a term in $e_n^2$ , and dividing that miss by the slope converts it directly into the next error: one squaring per step.

The secant method replaces $f$ by the line through the last two points. That line misses $f$ by a term in $e_n e_{n-1}$ rather than $e_n^2$ , so each step multiplies two old errors instead of squaring one. Taking logarithms turns the recursion $|e_{n+1}| \sim K|e_n||e_{n-1}|$ into the Fibonacci recursion, and the golden ratio appears because that is how fast the Fibonacci numbers grow.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §27.1.3-27.1.4 and §27.2, pp. 990-994

(a) Bisection.

The bracket survives. We show inductively that $f(a_n) f(b_n) \le 0$ for every $n$ ; it holds at $n=1$ because $f(a)f(b) < 0$ . Suppose it holds at $n$ . If $f(a_n)f(x_n) \le 0$ , the rule puts $[a_{n+1},b_{n+1}] = [a_n,x_n]$ and the invariant is the test itself. Otherwise $f(a_n)f(x_n) > 0$ , and multiplying the inductive hypothesis $f(a_n)f(b_n) \le 0$ by this positive number gives $f(a_n)^2 f(x_n) f(b_n) \le 0$ ; since $f(a_n) \neq 0$ this forces $f(x_n)f(b_n) \le 0$ , which is the invariant for $[a_{n+1},b_{n+1}] = [x_n, b_n]$ . In either case the half kept is one across which $f$ does not keep a strict sign. The inequality has to be the non-strict one from here on: as soon as some $f(x_n)$ vanishes the product is $0$ rather than negative, and that is exactly the case the rule must not discard.
The lengths. By construction $b_{n+1} - a_{n+1} = \tfrac12 (b_n - a_n)$ , so $b_n - a_n = (b-a)/2^{\,n-1}$ .
A common limit. $(a_n)$ is non-decreasing, $(b_n)$ is non-increasing, $a_n < b_n$ , and the gap tends to $0$ ; so both converge to one point $\xi \in [a,b]$ , and $a_n \le \xi \le b_n$ for every $n$ .
It is a root. $f$ is continuous, so $f(a_n) \to f(\xi)$ and $f(b_n) \to f(\xi)$ . Passing to the limit in $f(a_n)f(b_n) \le 0$ gives $f(\xi)^2 \le 0$ , hence $f(\xi) = 0$ . This is the step that uses continuity, and it is the intermediate value theorem in the form we need.
The bound. $x_n$ is the midpoint of $[a_n,b_n]$ and $\xi$ lies in that interval, so $|x_n - \xi| \le \tfrac12 (b_n - a_n) = (b-a)/2^{\,n}$ .

(b) Newton-Raphson. Fix $\delta_0 > 0$ with $J = [\xi - \delta_0, \xi + \delta_0] \subset I$ . Since $f'$ is continuous and $f'(\xi) \neq 0$ , shrinking $\delta_0$ makes $f'$ non-vanishing on $J$ . Put

K = \frac{\max_{J} |f''|}{2 \min_{J} |f'|} < \infty ,

both extrema being attained because $J$ is compact and $f', f''$ continuous.

One step, exactly. Let $x_n \in J$ . Apply Taylor's theorem with the Lagrange remainder (Taylor's theorem) to $f$ about $x_n$ , evaluated at $\xi$ : there is $\eta_n$ between $x_n$ and $\xi$ with

0 = f(\xi) = f(x_n) + f'(x_n)(\xi - x_n) + \tfrac12 f''(\eta_n) (\xi - x_n)^2 .

Rearrange into the iteration. Divide by $f'(x_n) \neq 0$ and use $x_{n+1} = x_n - f(x_n)/f'(x_n)$ :

0 = (x_n - x_{n+1}) - e_n + \frac{f''(\eta_n)}{2f'(x_n)} e_n^{2}, \qquad\text{that is}\qquad e_{n+1} = \frac{f''(\eta_n)}{2 f'(x_n)}\, e_n^{2} .

This identity is exact; no term has been discarded. It is the whole proof, and everything else is bookkeeping. 3. Staying put. Choose $\delta \le \delta_0$ with $\theta := K\delta < 1$ . If $|e_n| \le \delta$ then $|e_{n+1}| \le K e_n^2 \le K\delta |e_n| = \theta |e_n| \le \delta$ , so the iterates never leave $[\xi-\delta,\xi+\delta]$ , where $f'$ does not vanish and the next step is defined. 4. Convergence, and how fast. Multiplying step 2 by $K$ gives $K|e_{n+1}| \le (K|e_n|)^2$ , so by induction $K|e_n| \le \theta^{\,2^{n}} \to 0$ : the error is squashed by one squaring per step, which is the doubling of correct digits. Finally $x_n \to \xi$ forces $\eta_n \to \xi$ (it is trapped between $x_n$ and $\xi$ ) and $f'(x_n) \to f'(\xi)$ by continuity, so the exact relation of step 2 gives $e_{n+1}/e_n^2 \to f''(\xi)/2f'(\xi)$ .

(c) Secant method. The tangent of Newton is replaced by a chord, and the Taylor remainder by the error of linear interpolation.

Lemma (error of a chord). Let $f$ be twice continuously differentiable on an interval $J$ , let $u \neq v$ in $J$ , let $p$ be the line through $(u,f(u))$ and $(v,f(v))$ , and let $z \in J$ . Then $f(z) - p(z) = \tfrac12 f''(\eta)(z-u)(z-v)$ for some $\eta \in J$ .

Proof of the lemma. For $z \in \{u,v\}$ both sides vanish. Otherwise set $K_z = (f(z)-p(z))/((z-u)(z-v))$ and $g(x) = f(x) - p(x) - K_z (x-u)(x-v)$ . Then $g$ vanishes at the three distinct points $u, v, z$ , so by Rolle's theorem $g'$ vanishes at two distinct points between them and $g''$ at one, say $\eta$ . Since $p'' \equiv 0$ , $g''(\eta) = f''(\eta) - 2K_z = 0$ . $\square$ (This is the case $n=1$ of the interpolation error formula of Polynomial and spline interpolation.)

One step, exactly. Let $p$ be the chord through $(x_{n-1}, f(x_{n-1}))$ and $(x_n, f(x_n))$ , with slope $s_n = (f(x_n)-f(x_{n-1}))/(x_n - x_{n-1})$ . The secant formula says exactly that $x_{n+1}$ is the root of $p$ , so $p(x) = s_n(x - x_{n+1})$ , and at $x = \xi$ ,

p(\xi) = s_n(\xi - x_{n+1}) = -s_n e_{n+1} .

Insert the lemma at $z = \xi$ . Since $f(\xi) = 0$ ,

-p(\xi) = f(\xi) - p(\xi) = \tfrac12 f''(\eta_n)(\xi - x_{n-1})(\xi - x_n) = \tfrac12 f''(\eta_n)\, e_{n-1} e_n ,

so, using step 1,

e_{n+1} = \frac{f''(\eta_n)}{2 s_n}\, e_{n-1} e_n .

The slope is a derivative. By the mean value theorem $s_n = f'(\zeta_n)$ for some $\zeta_n$ between $x_{n-1}$ and $x_n$ ; in particular $s_n \neq 0$ once both iterates are in a neighbourhood of $\xi$ on which $f'$ does not vanish, so the step is defined.
Staying put, and convergence. With $J$ , $K$ as in part (b), step 2 gives $|e_{n+1}| \le K |e_n||e_{n-1}|$ . Choose $\delta$ with $\theta := K\delta < 1$ . If $|e_{n-1}|, |e_n| \le \delta$ then $|e_{n+1}| \le K\delta^2 = \theta\delta \le \delta$ , so all iterates stay in $[\xi-\delta,\xi+\delta]$ . Writing $\epsilon_n = K|e_n|$ , the recursion reads $\epsilon_{n+1} \le \epsilon_n \epsilon_{n-1}$ with $\epsilon_0, \epsilon_1 \le \theta$ , so by induction $\epsilon_n \le \theta^{\,F_n}$ where $F_0 = F_1 = 1$ and $F_{n+1} = F_n + F_{n-1}$ . Since $F_n \to \infty$ , $e_n \to 0$ .
Where the golden ratio comes from. The Fibonacci numbers satisfy $F_n = (\varphi^{\,n+1} - (-\varphi)^{-n-1})/\sqrt5$ with $\varphi = (1+\sqrt5)/2$ , so $F_n$ grows like $\varphi^{\,n}$ : the number of correct digits is multiplied by $\varphi$ at each step, which is the sense in which the order is $\varphi$ . For the sharper statement, assume $|e_{n+1}| / |e_n|^{p} \to C$ . Then $|e_n| \approx C|e_{n-1}|^{p}$ and step 2 gives

C\,|e_n|^{p} \approx K_\ast |e_n| |e_{n-1}| \quad\text{with } K_\ast = \left|\frac{f''(\xi)}{2f'(\xi)}\right| .

Substituting $|e_n| = C|e_{n-1}|^{p}$ on both sides,

C^{\,p+1} |e_{n-1}|^{\,p^{2}} \approx K_\ast C\, |e_{n-1}|^{\,p+1} .

For this to hold as $|e_{n-1}| \to 0$ the exponents must agree, $p^2 = p+1$ , whose positive root is $p = \varphi$ ; the coefficients then give $C^{p} = K_\ast$ , that is $C = K_\ast^{1/p} = K_\ast^{\,p-1}$ , using $1/\varphi = \varphi - 1$ . $\square$

Remark(Where the hypotheses are used, and how Newton fails without them)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.1.4, the behaviour near a stationary point, p. 991

Bisection uses continuity and nothing else — no derivative, no smoothness, not even a formula for $f$ . What it does need is the sign change $f(a)f(b)<0$ , and that is a real restriction: $f(x) = x^2$ has a root on $[-1,1]$ and no sign change, so bisection cannot be started, while $f(x) = 1/x$ has a sign change on $[-2,1]$ and no root, so what the method converges to need not be a root if $f$ is discontinuous: the rule produces $x_1 = -\tfrac12$ , $x_2 = \tfrac14$ , $x_3 = -\tfrac18$ , $x_4 = \tfrac1{16}, \dots \to 0$ , which is the pole and not a zero. (The bracket has to be chosen so that no midpoint is the pole itself; on $[-1,1]$ the first midpoint is $0$ and $f(x_1)$ is undefined.) Step 4 of the proof is where continuity enters and where the $1/x$ example breaks down. If the bracket contains several roots, bisection finds one of them and gives no warning about the others.
Newton uses $f'(\xi) \neq 0$ to divide in step 2 of the proof. At a root of multiplicity $m \ge 2$ it still converges, but only linearly. Write $f(x) = (x-\xi)^m g(x)$ with $g(\xi) \neq 0$ ; then $f/f' = (x-\xi)g/\bigl(mg + (x-\xi)g'\bigr)$ and

e_{n+1} = e_n\left(1 - \frac{g(x_n)}{m g(x_n) + e_n g'(x_n)}\right) \longrightarrow \Bigl(1 - \frac1m\Bigr) e_n ,

so the error is multiplied by $1-1/m$ each step: at a double root it merely halves, and reaching $10^{-15}$ from $10^{-1}$ takes about $47$ steps instead of $4$ . The same algebra shows the cure when $m$ is known: iterating with $x_{n+1} = x_n - m f(x_n)/f'(x_n)$ gives $e_{n+1} = e_n^2 g'(x_n)/(mg(x_n) + e_n g'(x_n)) \to e_n^2 g'(\xi)/(m g(\xi))$ , and quadratic convergence is restored.

Newton uses $f \in C^2$ for the Lagrange remainder in step 1. Dropping it costs the rate and sometimes the convergence: $f(x) = \operatorname{sign}(x)|x|^{1/3}$ has the single root $0$ and is differentiable away from it, but $f/f' = 3x$ , so Newton reads $x_{n+1} = -2x_n$ and the iterate doubles in modulus and alternates in sign from any starting point.
Newton uses " $x_0$ close enough", which is step 3 of the proof, and this is the hypothesis that fails in practice. Three ways:
1. A small derivative throws the iterate far. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.1.4, p. 991) note that the scheme behaves badly when some $x_n$ falls near a stationary point of $f$ , where $f'(x_n) \approx 0$ , and that a nearby maximum or minimum can leave the iterates hopping across the root rather than closing on it.
2. It can cycle. For $f(x) = x^3 - 2x + 2$ and $x_0 = 0$ : $f(0) = 2$ , $f'(0) = -2$ , so $x_1 = 1$ ; then $f(1) = 1$ , $f'(1) = 1$ , so $x_2 = 0$ . The iteration repeats $0, 1, 0, 1, \dots$ for ever, and the only real root, $-1.769292354\ldots$ , is never approached.
3. It can converge to the wrong root, when $f$ has several and the tangent at $x_0$ points at a different one. Nothing in the method expresses a preference. The contraction-mapping view explains why all of this is local: Kreyszig (Kreyszig (1978), §5.1, Problem 15, pp. 305-306) sets the exercise of showing that $g(x) = x - f(x)/f'(x)$ is a contraction in some neighbourhood of a simple zero, so Banach's fixed point theorem (The contraction mapping theorem) applies there and says nothing anywhere else.
The secant method needs $f(x_n) \neq f(x_{n-1})$ to divide, and it keeps no bracket, so an iterate can leave the neighbourhood and never come back. Its bracketed cousin, in which the two points are always kept on opposite sides of the root — Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.1.2, Eq. (27.8), pp. 988-989) call it linear interpolation; it is also known as false position — is safe but slow, because, as they observe, when the curvature of $f$ keeps one sign near the root one endpoint tends to stay frozen, and the convergence degrades to linear.
Floating-point arithmetic sets a floor on all three. Near a simple root the computed $f$ carries an absolute error of roughly $u$ times the size of the terms making it up (Floating-point arithmetic and error), so $f$ cannot be distinguished from zero within a band of width $\approx \text{noise}/|f'(\xi)|$ around $\xi$ , and no method can do better. For the example below the noise in $f$ is about $1.1\times10^{-16}\times 18.4 \approx 2\times10^{-15}$ and $|f'(\xi)| = 11.16$ , giving a floor of about $1.8\times10^{-16}$ . One unit in the last place at $2.09$ is $2^{-51} = 4.4\times10^{-16}$ , so the floor sits at about $0.4$ of an ulp: full double precision is attainable here, and the last bit is the only thing in doubt. At a double root the floor is $\sqrt{2\,\text{noise}/|f''(\xi)|}$ instead, typically around $10^{-8}$ : half the digits are simply not there, whatever the algorithm.

Method(Solving

f(x)=0

in one variable)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.1, the four schemes and the closing recommendation to combine them, pp. 985-992

Bracket first. Evaluate $f$ on a coarse grid over the range of interest and find a pair $a < b$ with $f(a)f(b) < 0$ . Plot $f$ if you can; most root-finding disasters are failures to look at the function. If $f$ has a root of even multiplicity there is no sign change, and you must bracket the minimum of $|f|$ instead.
Bisect a few times to get inside the basin of the fast methods, say until $b_n - a_n$ is a per cent of the original bracket. This costs $\lceil \log_2 100 \rceil = 7$ evaluations and buys the guarantee.
Switch to a fast method.
- If $f'$ is available cheaply and analytically, iterate $x_{n+1} = x_n - f(x_n)/f'(x_n)$ .
- If not, iterate the secant formula from the current bracket endpoints. Never form $f'$ by a finite difference merely to feed Newton: that costs two evaluations per step, the same as the secant method gets two steps out of, and it introduces a truncation error the secant method does not have.
Guard the fast step. Keep the bracket $[a_n,b_n]$ up to date. If a proposed $x_{n+1}$ falls outside it, or if the step has not at least halved, $|x_{n+1}-x_n| > \tfrac12|x_n - x_{n-1}|$ , discard it and bisect instead. Both tests use only computed quantities; the error $|x_n - \xi|$ itself is of course not available, which is why the step length and the bracket width stand in for it. This hybrid inherits the guarantee of bisection and the speed of Newton; some combination of a bracketing method with Newton-Raphson at the end is what Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.1, p. 992) recommend for most situations.
Stop on the step, not on the residual. Terminate when $|x_{n+1}-x_n| \le \tau(1+|x_{n+1}|)$ , with a relative tolerance $\tau$ no smaller than a few times the unit roundoff, and cap the iteration count. A small $|f(x_n)|$ on its own means nothing: it is small whenever $f'$ is small, whether or not $x_n$ is near a root.
Choose by cost. Per function evaluation, Newton multiplies the number of correct digits by $2^{1/2} = 1.414$ (two evaluations, $f$ and $f'$ , buy one squaring) while the secant method multiplies it by $\varphi = 1.618$ (one evaluation per step). So when $f'$ costs about as much as $f$ , the secant method is the more efficient of the two; Newton wins when $f'$ is nearly free, as it is for a polynomial evaluated by Horner's rule.

Example(The three methods on

x^{3}-2x-5=0

)Original

Take $f(x) = x^{3}-2x-5$ , the equation Newton himself used. Then $f(2) = -1 < 0$ and $f(3) = 16 > 0$ , so there is a root in $(2,3)$ ; and $f'(x) = 3x^{2}-2 \ge 10 > 0$ there, so the root is unique in the interval and simple. To twenty digits,

\xi = 2.0945514815423265915 .

The asymptotic error constant predicted by the theorem is

\frac{f''(\xi)}{2f'(\xi)} = \frac{6\xi}{2(3\xi^{2}-2)} = \frac{3\xi}{3\xi^{2}-2} = \frac{6.28365444}{11.16143773} = 0.5629789458 .

Bisection on $[2,3]$ . $x_1 = 2.5$ , $x_2 = 2.25$ , $x_3 = 2.125$ , $x_4 = 2.0625$ , $x_5 = 2.09375$ , $x_6 = 2.109375$ ; the errors are $0.4054$ , $0.1554$ , $0.03045$ , $-0.03205$ , $-0.0008015$ , $0.01482$ . They are not monotone — only the bound $|x_n-\xi| \le 2^{-n}$ is — and $x_5$ happens to be far better than $x_6$ . To guarantee $10^{-10}$ requires $2^{-n} \le 10^{-10}$ , that is $n \ge \log_2 10^{10} = 33.22$ , so $34$ evaluations of $f$ .

Newton from $x_0 = 2$ . $f(2) = -1$ , $f'(2) = 10$ , so $x_1 = 2 - (-1)/10 = 2.1$ . Next, $f(2.1) = 9.261 - 4.2 - 5 = 0.061$ and $f'(2.1) = 3(4.41)-2 = 11.23$ , so $x_2 = 2.1 - 0.061/11.23 = 2.0945681211041852$ . Continuing:

$n$	$x_n$	$e_n = x_n - \xi$	$e_n / e_{n-1}^{2}$
0	$2$	$-9.4551\times10^{-2}$
1	$2.1$	$5.44852\times10^{-3}$	$0.60946$
2	$2.09456812110418522$	$1.66396\times10^{-5}$	$0.56051$
3	$2.09455148169819930$	$1.55873\times10^{-10}$	$0.56297$
4	$2.09455148154232659$	$1.36783\times10^{-20}$	$0.5629790$

The last column converges to $0.5629789458$ , the predicted constant, and $\log_{10}|e_n| = -1.02, -2.26, -4.78, -9.81, -19.86$ doubles at each step. Reaching $1.6\times10^{-10}$ at $x_3$ costs three evaluations of $f$ and three of $f'$ ; bisection needed $34$ evaluations for $10^{-10}$ .

Secant from $x_0 = 2$ , $x_1 = 3$ . The first step is $x_2 = 3 - 16 \times (3-2)/(16-(-1)) = 3 - 16/17 = 35/17 = 2.0588235294117647$ . Continuing:

$n$	$x_n$	$e_n$	$e_n /(e_{n-1}e_{n-2})$
2	$2.05882352941176471$	$-3.5728\times10^{-2}$	$0.41733$
3	$2.08126365984502285$	$-1.32878\times10^{-2}$	$0.41075$
4	$2.09482414609405235$	$2.72665\times10^{-4}$	$0.57434$
5	$2.09454943103524732$	$-2.05051\times10^{-6}$	$0.56595$
6	$2.09455148122759912$	$-3.14727\times10^{-10}$	$0.5629174$
7	$2.09455148154232695$	$3.63319\times10^{-16}$	$0.5629794$

Again the last column tends to $0.5629789458$ . The measured order, $\log|e_{n+1}|/\log|e_n|$ , takes the successive values $1.596$ , $1.671$ , $1.625$ , $1.632$ , $1.623$ , $1.622$ , $1.620$ as the iteration proceeds, against the predicted $\varphi = 1.6180$ .

The cost comparison. Reaching $x_6$ , with error $3.1\times10^{-10}$ , takes six evaluations of $f$ (two to start, one per step thereafter) and none of $f'$ . Newton reached the comparable error $1.6\times10^{-10}$ at $x_3$ , also after six evaluations — but three of them were of $f'$ . With a cubic polynomial the derivative is nearly free and the two methods tie; with an $f$ whose derivative must be assembled by hand or by automatic differentiation, the secant method does the same work for half the code and half the cost.

Application(Wien's displacement law)Original

The spectral radiance of a black body per unit wavelength is Planck's function

B_\lambda(T) = \frac{2hc^{2}}{\lambda^{5}}\,\frac{1}{e^{\,hc/\lambda k_B T}-1} .

Where does it peak? Write $x = hc/(\lambda k_B T)$ , so that $\lambda \propto 1/x$ and $B_\lambda \propto x^{5}/(e^{x}-1)$ . Setting the derivative to zero,

\frac{d}{dx}\log\frac{x^5}{e^x-1} = \frac5x - \frac{e^{x}}{e^{x}-1} = 0 \qquad\Longleftrightarrow\qquad 5\bigl(1 - e^{-x}\bigr) = x ,

so the peak is at the positive root of

g(x) = x - 5 + 5e^{-x} = 0 .

There is no closed form, which is why this calculation is a root-finding problem and has been since 1893. The root is simple: $g'(x) = 1 - 5e^{-x}$ , and at a root $5e^{-x} = 5-x$ , so $g'(x_\ast) = x_\ast - 4$ , which is not zero unless $x_\ast = 4$ — and $g(4) = 5e^{-4} - 1 = -0.9084 \neq 0$ . So Newton's theorem applies.

Newton from $x_0 = 5$ . $g(5) = 5e^{-5} = 0.03368973$ and $g'(5) = 1-5e^{-5} = 0.96631027$ , so

x_1 = 5 - \frac{0.03368973}{0.96631027} = 4.96513569583650 .

Then $g(x_1) = 2.07153\times10^{-5}$ gives $x_2 = 4.96511423175260$ , and $g(x_2) = 8.03594\times10^{-12}$ gives

x_3 = 4.96511423174428 ,

with $g(x_3) = 1.2\times10^{-24}$ : sixteen correct digits after three steps, and the errors $3.5\times10^{-2}$ , $2.1\times10^{-5}$ , $8.3\times10^{-12}$ , $1.2\times10^{-24}$ show the doubling plainly. Bisection from the bracket $[4,6]$ would need $\lceil \log_2\bigl(2/(5\times10^{-16})\bigr)\rceil = \lceil 51.83 \rceil = 52$ evaluations for the same accuracy.

The physical constant. With the SI defining values $h = 6.62607015\times10^{-34}\ \mathrm{J\,s}$ , $c = 299792458\ \mathrm{m\,s^{-1}}$ and $k_B = 1.380649\times10^{-23}\ \mathrm{J\,K^{-1}}$ , so that $hc/k_B = 1.4387768775\times10^{-2}\ \mathrm{m\,K}$ , the Wien displacement constant is

b = \lambda_{\max} T = \frac{hc}{k_B x_\ast} = \frac{1.4387768775\times10^{-2}}{4.96511423174428}\ \mathrm{m\,K} = 2.8977719552\times10^{-3}\ \mathrm{m\,K} ,

which is the value tabulated by CODATA. For the Sun, with effective temperature $T = 5772\ \mathrm{K}$ , this puts the peak of the wavelength spectrum at

\lambda_{\max} = \frac{2.8977719552\times10^{-3}}{5772}\ \mathrm{m} = 5.0204\times10^{-7}\ \mathrm{m} = 502.04\ \mathrm{nm} ,

in the green. The same equation with $5$ replaced by $3$ gives the peak of the spectrum per unit frequency, at a different wavelength — the two peaks disagree because a spectral density transforms with the Jacobian $d\nu/d\lambda$ , and this is the standard trap that the root-finding step makes concrete rather than removing. Root finding of exactly this kind appears throughout physics whenever a dispersion relation, an equation of state or a transcendental eigenvalue condition has to be inverted: the quantisation condition for a finite square well, the Kepler equation $E - e\sin E = M$ of orbital mechanics, whose root is simple for every eccentricity $e<1$ because the derivative $1 - e\cos E$ is at least $1-e>0$ , and the bisection-plus-Newton solvers inside every equation-of-state routine.

Method

Polynomial and spline interpolation

Adapted

Add to PDF

The unique polynomial through n+1 points, in Lagrange and Newton form, with the exact error formula that explains both its accuracy near the nodes and Runge's divergence away from them, and the cubic spline that replaces it when the nodes are many and equally spaced.

Prerequisites: Taylor's theorem

NoteOriginal

What it means. You have a function at a finite set of points — a table, a measurement, the output of an expensive code — and you want it in between. There is exactly one polynomial of degree $n$ through $n+1$ points, and an exact formula for what it misses. The formula carries two lessons that pull in opposite directions. Near the nodes, and for a handful of them, polynomial interpolation is extremely accurate and underwrites almost every other method in this part: quadrature rules, finite differences and finite elements are all obtained by integrating or differentiating an interpolant. Pushed to many equally spaced nodes it fails spectacularly, and the cure is not a better polynomial but a different object: a chain of low-degree pieces glued smoothly, the spline.

Definition(Interpolant, Lagrange basis, divided differences)Original

Let $x_0 < x_1 < \cdots < x_n$ be nodes, distinct points of an interval $[a,b]$ , and let $f_0, \dots, f_n$ be prescribed values; when they come from a function we write $f_i = f(x_i)$ . A polynomial $p$ interpolates these data if $p(x_i) = f_i$ for $i = 0,\dots,n$ . Write $\mathcal{P}_n$ for the real polynomials of degree at most $n$ , a vector space of dimension $n+1$ .

The Lagrange basis for the nodes is

\ell_k(x) = \prod_{\substack{j=0 \\ j \neq k}}^{n} \frac{x - x_j}{x_k - x_j}, \qquad k = 0,\dots,n ,

each $\ell_k \in \mathcal{P}_n$ with $\ell_k(x_i) = \delta_{ik}$ ; the denominators are non-zero precisely because the nodes are distinct. The Lagrange form of the interpolant is

p(x) = \sum_{k=0}^{n} f_k \,\ell_k(x) .

The divided difference $f[x_i, x_{i+1}, \dots, x_{i+k}]$ is defined as the coefficient of $x^{k}$ in the polynomial of degree at most $k$ that interpolates the data at the $k+1$ nodes $x_i, \dots, x_{i+k}$ — its leading coefficient. In particular $f[x_i] = f_i$ . The Newton form of the interpolant is

p(x) = \sum_{k=0}^{n} f[x_0,\dots,x_k] \prod_{j=0}^{k-1}(x - x_j) = f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) + \cdots,

the empty product being $1$ . Finally write

w(x) = \prod_{i=0}^{n}(x-x_i) ,

the monic node polynomial of degree $n+1$ , which governs the error.

Method(Existence, uniqueness, the Newton recursion and the error)Original

Let $x_0, \dots, x_n$ be distinct nodes in $[a,b]$ .

(i) Existence and uniqueness. For any values $f_0,\dots,f_n$ there is exactly one $p \in \mathcal{P}_n$ with $p(x_i) = f_i$ for every $i$ , and it is given by either the Lagrange form or the Newton form above; the two are the same polynomial written in two bases.

(ii) The divided-difference recursion. For $k \ge 1$ ,

f[x_i, \dots, x_{i+k}] = \frac{f[x_{i+1},\dots,x_{i+k}] - f[x_i,\dots,x_{i+k-1}]}{x_{i+k}-x_i} .

(iii) The error. Suppose in addition that $f$ is $n+1$ times continuously differentiable on $[a,b]$ and $f_i = f(x_i)$ . Then for every $x \in [a,b]$ there is a point $\xi_x$ in the open interval spanned by $x, x_0, \dots, x_n$ with

f(x) - p(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\; w(x), \qquad w(x) = \prod_{i=0}^{n}(x-x_i) .

In particular

\max_{a \le x \le b} |f(x)-p(x)| \ \le \ \frac{\max_{a\le t\le b}|f^{(n+1)}(t)|}{(n+1)!}\ \max_{a \le x \le b}|w(x)| ,

a product of one factor that depends only on $f$ and one that depends only on where the nodes are put.

Idea of the proofOriginal

Uniqueness is the engine, and it is one line: two interpolants of the same data differ by a polynomial of degree at most $n$ with $n+1$ roots, so they are equal. Existence is then free, because the Lagrange basis writes an interpolant down. The Newton form and its recursion are uniqueness used twice: the difference of the interpolants on $x_0,\dots,x_k$ and on $x_0,\dots,x_{k-1}$ vanishes at the first $k$ nodes and has degree $k$ , so it must be the leading coefficient times $(x-x_0)\cdots(x-x_{k-1})$ — which telescopes into the Newton form.

The error formula is Rolle's theorem counted carefully. Fix a point $x$ away from the nodes and define the constant that would make the formula true; then the function $f - p - (\text{that constant}) w$ has $n+2$ zeros, namely the $n+1$ nodes and $x$ . Each application of Rolle costs one zero and gains one derivative, so after $n+1$ applications something of the form $f^{(n+1)} - (n+1)!\,(\text{constant})$ has a zero, and that identifies the constant.

ProofOriginal

(i) Uniqueness. Let $p, q \in \mathcal{P}_n$ both interpolate the data. Then $p-q \in \mathcal{P}_n$ vanishes at the $n+1$ distinct nodes. A non-zero real polynomial of degree at most $n$ has at most $n$ roots, so $p - q \equiv 0$ . This is the only place the distinctness of the nodes is used, and it is used decisively.

Existence. Each $\ell_k$ lies in $\mathcal{P}_n$ , being a product of $n$ linear factors, and $\ell_k(x_i) = \delta_{ik}$ because the factor $(x_i - x_i)$ appears in the numerator when $i \neq k$ while at $x_i = x_k$ numerator and denominator coincide. Hence $\sum_k f_k \ell_k$ lies in $\mathcal{P}_n$ and takes the value $f_i$ at $x_i$ .

The Newton form. Let $p_k \in \mathcal{P}_k$ interpolate the data at $x_0,\dots,x_k$ , so $p_n = p$ and $p_0 \equiv f_0$ . The difference $d_k = p_k - p_{k-1}$ lies in $\mathcal{P}_k$ and vanishes at $x_0, \dots, x_{k-1}$ , because both interpolants agree with the data there. So $d_k$ has the $k$ roots $x_0,\dots,x_{k-1}$ and degree at most $k$ , which forces

d_k(x) = c_k \prod_{j=0}^{k-1}(x-x_j)

for a constant $c_k$ ; comparing coefficients of $x^{k}$ , $c_k$ is the leading coefficient of $p_k$ , which is by definition $f[x_0,\dots,x_k]$ . Summing $p = p_0 + \sum_{k=1}^{n} d_k$ gives the Newton form.

(ii) The recursion. Let $u \in \mathcal{P}_{k-1}$ interpolate at $x_i,\dots,x_{i+k-1}$ and $v \in \mathcal{P}_{k-1}$ interpolate at $x_{i+1},\dots,x_{i+k}$ . Define

q(x) = \frac{(x - x_i)\,v(x) - (x-x_{i+k})\,u(x)}{x_{i+k}-x_i} .

Then $q \in \mathcal{P}_k$ . At $x = x_i$ the first term vanishes and $q(x_i) = u(x_i) = f_i$ ; at $x = x_{i+k}$ the second vanishes and $q(x_{i+k}) = v(x_{i+k}) = f_{i+k}$ ; at an interior node $x_m$ with $i<m<i+k$ both $u$ and $v$ equal $f_m$ , so $q(x_m) = f_m \bigl((x_m-x_i)-(x_m-x_{i+k})\bigr)/(x_{i+k}-x_i) = f_m$ . So $q$ interpolates at all $k+1$ nodes, and by uniqueness it is the interpolant. Reading off the coefficient of $x^{k}$ on both sides gives the stated recursion, the left side being $f[x_i,\dots,x_{i+k}]$ and the right $\bigl(f[x_{i+1},\dots,x_{i+k}] - f[x_i,\dots,x_{i+k-1}]\bigr)/(x_{i+k}-x_i)$ .

(iii) The error. If $x$ is one of the nodes both sides vanish, so let $x$ be none of them; then $w(x) \neq 0$ . Put

K = \frac{f(x)-p(x)}{w(x)}, \qquad g(t) = f(t) - p(t) - K\,w(t) .

$g$ is $n+1$ times continuously differentiable on $[a,b]$ , since $f$ is and $p, w$ are polynomials.
$g$ vanishes at the $n+2$ distinct points $x_0,\dots,x_n$ (where $f = p$ and $w = 0$ ) and $x$ (by the choice of $K$ ).
Rolle's theorem applied to each of the $n+1$ intervals between consecutive zeros gives $n+1$ distinct zeros of $g'$ , strictly inside the interval $I_x$ spanned by $x, x_0,\dots,x_n$ . (This is the engine that also produces the Lagrange remainder in Taylor's theorem, and the two results coincide in the limit: running all the nodes together at one point $x_0$ turns $p$ into the Taylor polynomial of order $n$ and $w$ into $(x-x_0)^{n+1}$ , and the formula below into the Lagrange remainder.) Repeating, $g''$ has $n$ zeros there, and after $n+1$ applications $g^{(n+1)}$ has at least one zero $\xi_x$ in the interior of $I_x$ .
Differentiate $g$ exactly $n+1$ times: $p^{(n+1)} \equiv 0$ because $\deg p \le n$ , and $w^{(n+1)} \equiv (n+1)!$ because $w$ is monic of degree $n+1$ . So $g^{(n+1)}(t) = f^{(n+1)}(t) - K(n+1)!$ , and $g^{(n+1)}(\xi_x) = 0$ gives $K = f^{(n+1)}(\xi_x)/(n+1)!$ .

Substituting $K$ back into its definition gives the formula. The bound follows by taking absolute values and maximising the two factors separately. $\square$

Remark(Where the hypotheses are used, and Runge's phenomenon)AdaptedAdapted from Kreyszig (1978), §6.6, the remark that a single interpolating polynomial is good near the nodes but may depart badly away from them, pp. 357-358

Distinct nodes. Used in uniqueness and in the Lagrange denominators. If two nodes merge, the interpolation conditions become inconsistent (different values) or redundant (equal values); the right generalisation is Hermite interpolation, which prescribes derivatives at repeated nodes, and the spline construction below uses the two-node case of it.
$f \in C^{n+1}$ . Used in steps 1 and 4 of part (iii). Without it the formula says nothing: $|x|$ on $[-1,1]$ is continuous and interpolated by polynomials of every degree, and the equispaced interpolants diverge.
The formula is an identity, not an estimate, but it is not a convergence theorem. As $n$ grows, $(n+1)!$ in the denominator races against two growing quantities: $\max|f^{(n+1)}|$ and $\max|w|$ . Nothing forces the race to be won. Indeed it is lost in the standard example, due to Runge. Take

f(x) = \frac{1}{1+x^{2}} \quad\text{on } [-5,5] ,

a function that is smooth on the real line, and interpolate it at $n+1$ equally spaced nodes. Computing the maximum of $|f-p_n|$ on a grid of $20001$ points gives $1.9157$ at $n=10$ and $59.822$ at $n=20$ : not convergence but divergence, and all of it within one spacing of the endpoints, while near the middle the interpolant is excellent. Kreyszig (Kreyszig (1978), §6.6, pp. 357-358) records exactly this behaviour — good near the nodes, badly off away from them — as the reason to prefer piecewise polynomials.

It is the node polynomial that can be controlled, not $f^{(n+1)}$ . For the same $n = 20$ on $[-5,5]$ , $\max|w|$ is $1.1142\times10^{11}$ with equally spaced nodes. Choosing instead the Chebyshev nodes, the zeros of $T_{21}$ rescaled to the interval,

x_i = 5\cos\frac{(2i+1)\pi}{2(n+1)}, \qquad i = 0,\dots,n,

makes $w$ the monic polynomial of degree $21$ of smallest maximum modulus on $[-5,5]$ , with $\max|w| = 2(10/4)^{21} = 4.5475\times10^{8}$ (the minimax property): smaller by a factor of $245$ . The computed interpolation errors fall accordingly, to $0.10915$ at $n=10$ and $0.015334$ at $n=20$ — now decreasing with $n$ . The nodes crowd towards the endpoints, which is precisely where the equispaced interpolant misbehaved.

Method(Cubic spline interpolation)AdaptedAdapted from Kreyszig (1978), §6.6, the definition of Y(P_n), Theorem 6.6-1 and the minimum property (6), pp. 357-360

Let $P : a = t_0 < t_1 < \cdots < t_n = b$ be a partition, with $h_j = t_{j+1}-t_j$ . A cubic spline on $P$ is a function $y \in C^{2}[a,b]$ whose restriction to each $[t_j, t_{j+1}]$ is a polynomial of degree at most $3$ ; write $S_3(P)$ for the set of them, a vector space.

(i) Existence and uniqueness (Kreyszig, Theorem 6.6-1 (Kreyszig (1978), §6.6, Thm 6.6-1, pp. 358-359)). Given values $f_0,\dots,f_n$ and two numbers $\kappa_a, \kappa_b$ , there is exactly one $y \in S_3(P)$ satisfying the $n+3$ conditions

y(t_j) = f_j \ \ (j = 0,\dots,n), \qquad y'(a) = \kappa_a, \qquad y'(b) = \kappa_b .

It is found by solving for the nodal slopes $k_j = y'(t_j)$ , $j = 1,\dots,n-1$ , the tridiagonal system

h_j\, k_{j-1} \;+\; 2\,(h_{j-1}+h_j)\, k_j \;+\; h_{j-1}\, k_{j+1} \;=\; 3\left( \frac{h_j}{h_{j-1}}\,\Delta_j + \frac{h_{j-1}}{h_j}\,\Delta_{j+1}\right), \qquad \Delta_j = f_j - f_{j-1},

with $k_0 = \kappa_a$ and $k_n = \kappa_b$ known, after which $y$ on $[t_j,t_{j+1}]$ is the unique cubic with the values $f_j, f_{j+1}$ and slopes $k_j, k_{j+1}$ at the two ends.

(ii) Minimum curvature (Kreyszig, Eq. (6) (Kreyszig (1978), §6.6, pp. 359-360)). If $f \in C^{2}[a,b]$ and $y$ is the spline of part (i) built from $f_j = f(t_j)$ , $\kappa_a = f'(a)$ , $\kappa_b = f'(b)$ , then

\int_a^b \bigl(y''(t)\bigr)^{2}\,dt \ \le \ \int_a^b \bigl(f''(t)\bigr)^{2}\,dt ,

with equality if and only if $f = y$ . Among all twice-differentiable functions through the data with those end slopes, the spline is the one of least total squared curvature — which is why a draughtsman's flexible rod, minimising its bending energy, takes a spline's shape, and why the object is named after the rod.

Idea of the proofOriginal

A cubic on one subinterval is determined by the values and the slopes at its two ends, so prescribing one slope per node already forces the pieces to agree in value and in slope at every node. The second derivative is the only thing left to arrange, and it gives one equation per interior node — exactly as many equations as there are unknown slopes. That the resulting matrix is invertible is the one analytic fact needed, and it comes free: the diagonal entry $2(h_{j-1}+h_j)$ beats the sum $h_j + h_{j-1}$ of its neighbours by a factor of two, and a strictly diagonally dominant matrix is non-singular.

The minimum property is an integration by parts. The spline's third derivative is constant on each subinterval, so integrating $y'''$ against $f'-y'$ piece by piece produces only the differences of $f-y$ at the nodes, and those are zero because $y$ interpolates $f$ . Hence $y''$ is orthogonal in $L^{2}$ to $f''-y''$ , and Pythagoras does the rest.

ProofAdaptedAdapted from Kreyszig (1978), §6.6, the proof of Theorem 6.6-1 through the Hermite cubic and the continuity of the second derivative, and the integration by parts giving (6), pp. 358-360

(i) Work on one subinterval first. Fix $j$ , write $h = h_j$ and $s = (t-t_j)/h \in [0,1]$ .

The Hermite cubic is unique. A cubic $p$ with prescribed $p(t_j), p(t_{j+1}), p'(t_j), p'(t_{j+1})$ is unique: the difference of two such cubics has double zeros at $t_j$ and at $t_{j+1}$ , so it is divisible by $(t-t_j)^2(t-t_{j+1})^2$ , a polynomial of degree $4$ ; a cubic divisible by a quartic is zero.
And it exists, given explicitly by

p(t) = f_j H_{00}(s) + f_{j+1}H_{01}(s) + h\,k_j\,H_{10}(s) + h\,k_{j+1}H_{11}(s),

H_{00} = 2s^{3}-3s^{2}+1, \quad H_{01} = -2s^{3}+3s^{2}, \quad H_{10} = s^{3}-2s^{2}+s, \quad H_{11} = s^{3}-s^{2} .

Checking the four conditions is direct: at $s=0$ the values are $H_{00}=1$ and $H_{01}=H_{10}=H_{11}=0$ ; at $s=1$ , $H_{01}=1$ and the others vanish; and $dp/dt = (dp/ds)/h$ with $H_{00}'(0)=H_{01}'(0)=H_{11}'(0)=0$ , $H_{10}'(0)=1$ , and $H_{00}'(1)=H_{01}'(1)=H_{10}'(1)=0$ , $H_{11}'(1)=1$ . 3. Two derivatives of the pieces. Differentiating twice in $t$ ,

p''(t) = \frac{1}{h^{2}}\Bigl( f_j(12s-6) + f_{j+1}(6-12s) + h k_j (6s-4) + h k_{j+1}(6s-2)\Bigr),

so at the two ends

p''(t_j) = \frac{6(f_{j+1}-f_j)}{h^{2}} - \frac{4k_j + 2k_{j+1}}{h}, \qquad p''(t_{j+1}) = -\frac{6(f_{j+1}-f_j)}{h^{2}} + \frac{2k_j + 4k_{j+1}}{h} .

Assemble. Define $y$ on $[t_j,t_{j+1}]$ by step 2, for every $j$ , using one slope $k_j$ per node. By construction the pieces agree in value and in first derivative at every interior node, so $y \in C^{1}[a,b]$ whatever the $k_j$ . Membership of $C^{2}$ is the condition that the two expressions of step 3 agree at each interior node $t_j$ , $1 \le j \le n-1$ : the right-hand end of the piece on $[t_{j-1},t_j]$ against the left-hand end of the piece on $[t_j,t_{j+1}]$ ,

-\frac{6\Delta_j}{h_{j-1}^{2}} + \frac{2k_{j-1}+4k_j}{h_{j-1}} = \frac{6\Delta_{j+1}}{h_j^{2}} - \frac{4k_j + 2k_{j+1}}{h_j} ,

with $\Delta_j = f_j - f_{j-1}$ . Multiplying by $\tfrac12 h_{j-1}h_j$ and collecting the $k$ 's on the left gives exactly the stated tridiagonal system. 5. The system is uniquely solvable. It has $n-1$ equations in the $n-1$ unknowns $k_1,\dots,k_{n-1}$ , with $k_0$ and $k_n$ known and moved to the right-hand side. In row $j$ the diagonal entry is $2(h_{j-1}+h_j)$ and the off-diagonal entries are $h_j$ and $h_{j-1}$ , all positive, so

|{\rm diagonal}| = 2(h_{j-1}+h_j) > h_{j-1} + h_j = \sum_{\text{off-diagonal}} |{\cdot}| :

the matrix $A$ is strictly diagonally dominant. Such a matrix is non-singular: if $Az = 0$ with $z \neq 0$ , pick an index $i$ with $|z_i|$ maximal, so $|z_i| > 0$ ; then

|a_{ii}||z_i| = \Bigl| \sum_{m \neq i} a_{im}z_m \Bigr| \le \Bigl(\sum_{m\neq i}|a_{im}|\Bigr)|z_i| < |a_{ii}||z_i| ,

a contradiction. (Kreyszig (Kreyszig (1978), §6.6, p. 359) draws the same conclusion from his Theorem 5.2-1 on iterative solution under a row-sum condition; the argument above is the direct one.) So the slopes exist and are unique, and with them so is $y$ by step 1. That proves (i).

(ii) Let $f \in C^{2}[a,b]$ and let $y$ be the spline built from its values and end slopes, so that $f - y$ vanishes at every node and $f'-y'$ vanishes at $a$ and $b$ .

Integrate by parts. Since $y''$ is continuous and $y'''$ exists and is continuous on the interior of each subinterval,

\int_a^b y''\,(f''-y'')\,dt = \Bigl[\,y''(f'-y')\,\Bigr]_a^b - \int_a^b y'''\,(f'-y')\,dt .

The boundary term vanishes because $f'-y'$ is zero at both ends; this is the step that uses the prescribed end slopes. 2. The remaining integral vanishes. On $[t_j,t_{j+1}]$ the function $y$ is a cubic, so $y''' = c_j$ is constant there, and

\int_{t_j}^{t_{j+1}} y'''(f'-y')\,dt = c_j\Bigl[(f-y)(t_{j+1}) - (f-y)(t_j)\Bigr] = 0

because $y$ interpolates $f$ at the nodes. Summing over $j$ gives $\int_a^b y'''(f'-y') \, dt = 0$ , hence $\int_a^b y''(f''-y'')\,dt = 0$ . 3. Pythagoras. Expand

\int_a^b (f''-y'')^{2} = \int_a^b (f'')^{2} - 2\int_a^b f''y'' + \int_a^b (y'')^{2} = \int_a^b (f'')^{2} - \int_a^b (y'')^{2} ,

where the last step used $\int y''(f''-y'') = 0$ twice over, in the form $\int f''y'' = \int (y'')^2$ . The left side is non-negative, which is the inequality. 4. Equality. Equality forces $\int_a^b (f''-y'')^{2} = 0$ , and since $f''-y''$ is continuous, $f'' \equiv y''$ . Then $f-y$ is affine; it vanishes at $t_0$ and $t_1$ , two distinct points, so $f \equiv y$ . $\square$

Remark(Why exactly two extra conditions, and which to choose)Original

The count in Theorem (i) — $n+1$ interpolation conditions plus $2$ end conditions — is not a coincidence but the dimension of $S_3(P)$ :

\dim S_3(P) = n+3 .

To see it, build a spline left to right. Its restriction to $[t_0,t_1]$ is an arbitrary cubic: $4$ parameters. On each later subinterval the cubic must match the value, the slope and the second derivative inherited at the left endpoint, which fixes three of its four coefficients and leaves one free (equivalently, the constant $y'''$ on that piece). That is $1$ parameter for each of the remaining $n-1$ subintervals, so $4 + (n-1) = n+3$ in all. Interpolation at the $n+1$ nodes is $n+1$ independent linear conditions, so exactly two degrees of freedom remain and exactly two end conditions are needed.

Which two is a modelling choice, and the theorem above covers one of them. Clamped (or complete) splines prescribe $y'(a)$ and $y'(b)$ , as Kreyszig does; use them when the end slopes are known, as they are when the data come from a function you can differentiate. Natural splines instead set $y''(a) = y''(b) = 0$ ; they need no extra information and are the ones that minimise the bending integral over all interpolants rather than over those with given end slopes, but they lose accuracy near the ends, since a smooth $f$ rarely has vanishing curvature there. Periodic splines match $y'$ and $y''$ at the two ends and are the right choice for data on a circle. Each choice replaces the first and last rows of the same tridiagonal system.

Method(Interpolating in practice)Original

Decide how many nodes the interpolant will use at once. Few (two to six) and local: use a polynomial. Many: use a spline, or a global polynomial only if you are free to place the nodes.
If you may choose the nodes, choose Chebyshev nodes, $x_i = \tfrac{a+b}{2} + \tfrac{b-a}{2}\cos\frac{(2i+1)\pi}{2(n+1)}$ , which minimise $\max|w|$ and so the error bound. If the nodes are given and equally spaced, do not raise the degree: interpolate locally on a short window of nodes and move the window.
Build the Newton form, not the Lagrange form, when you will evaluate at many points. Form the divided-difference table by the recursion, column by column; there are $n(n+1)/2$ entries beyond the data, each costing two subtractions and a division, so about $\tfrac32 n^{2}$ operations once. Evaluate by nesting,

p(x) = f[x_0] + (x-x_0)\Bigl( f[x_0,x_1] + (x-x_1)\bigl(f[x_0,x_1,x_2] + \cdots\bigr)\Bigr),

which costs $n$ multiplications, $n$ subtractions and $n$ additions: $3n$ operations per point, against about $2n^{2}$ for the Lagrange form evaluated as written. Adding one more node costs one more column, $O(n)$ work, with nothing rebuilt. 4. For a spline, assemble and solve the tridiagonal system. Never call a general linear solver: eliminating one subdiagonal entry costs a division and two multiply-subtracts, and back substitution three operations per row, so the whole solve is about $8n$ operations, against the roughly $n^{3}/3$ that Kreyszig (Kreyszig (1978), §5.2, p. 307) quotes for the usual direct methods on a dense system (Solving linear systems numerically). Then store the four Hermite coefficients per subinterval and evaluate by locating the subinterval with a binary search, $O(\log n)$ , followed by three multiply-adds. 5. Check the result against the error formula rather than against your expectations: compute $\max|w|$ on the interval of interest, bound $|f^{(n+1)}|$ if you can, and compare with the observed difference at a point that is not a node. If the two disagree by orders of magnitude, a node is wrong or the data are not samples of a smooth function. 6. Never interpolate noisy data. The interpolant reproduces the noise exactly and its derivatives amplify it. Fit instead, by least squares (Orthogonal projection and least squares) or with a smoothing spline.

Example(Quadratic interpolation and a cubic spline, on the same three points)Original

Take $f(x) = \ln x$ at the three nodes $x_0 = 1$ , $x_1 = 2$ , $x_2 = 4$ , so

f_0 = 0, \qquad f_1 = \ln 2 = 0.6931471806, \qquad f_2 = \ln 4 = 1.3862943611 ,

and estimate $\ln 3$ . The exact value is $1.0986122887$ .

The quadratic interpolant, Newton form. The divided-difference table:

f[x_0,x_1] = \frac{0.6931471806-0}{2-1} = 0.6931471806, \qquad f[x_1,x_2] = \frac{1.3862943611-0.6931471806}{4-2} = 0.3465735903,

f[x_0,x_1,x_2] = \frac{0.3465735903-0.6931471806}{4-1} = -0.1155245301 .

Hence $p(x) = 0.6931471806\,(x-1) - 0.1155245301\,(x-1)(x-2)$ and

p(3) = 0.6931471806 \times 2 - 0.1155245301 \times 2 \times 1 = 1.3862943611 - 0.2310490602 = 1.1552453009 .

As a check, the Lagrange form at $x=3$ uses $\ell_0(3) = \frac{(3-2)(3-4)}{(1-2)(1-4)} = -\tfrac13$ , $\ell_1(3) = \frac{(3-1)(3-4)}{(2-1)(2-4)} = 1$ , $\ell_2(3) = \frac{(3-1)(3-2)}{(4-1)(4-2)} = \tfrac13$ , giving $0 \cdot(-\tfrac13) + 0.6931471806 + \tfrac13 (1.3862943611) = 1.1552453009$ : the same number, as uniqueness requires.

The error, and the point $\xi$ it names. The actual error is

f(3)-p(3) = 1.0986122887 - 1.1552453009 = -0.0566330122 .

The formula predicts $f(3)-p(3) = \tfrac{1}{3!}f'''(\xi)\,w(3)$ with $f'''(x) = 2/x^{3}$ and $w(3) = (3-1)(3-2)(3-4) = -2$ , that is $-2/(3\xi^{3})$ . Solving $-2/(3\xi^{3}) = -0.0566330122$ gives $\xi^{3} = 11.7716971$ and

\xi = 2.2748164 ,

which does lie in $(1,4)$ , as the theorem says it must. The crude bound, with $\max_{[1,4]}|f'''| = 2$ attained at $x=1$ , is $|{\rm error}| \le \tfrac{2}{6}\times 2 = 0.667$ : true but twelve times too pessimistic, which is typical, since the bound maximises the two factors at different points.

The clamped cubic spline on the same data. Take $n = 2$ with $t_0=1, t_1=2, t_2=4$ , so $h_0 = 1$ , $h_1 = 2$ , and clamp with the true end slopes $\kappa_a = f'(1) = 1$ and $\kappa_b = f'(4) = 0.25$ . Here $\Delta_1 = f_1-f_0 = 0.6931471806$ and $\Delta_2 = f_2-f_1 = 0.6931471806$ . There is one interior node, so the tridiagonal system is a single equation ( $j = 1$ ):

h_1 k_0 + 2(h_0+h_1)k_1 + h_0 k_2 = 3\left(\frac{h_1}{h_0}\Delta_1 + \frac{h_0}{h_1}\Delta_2\right),

2(1) + 6k_1 + 1(0.25) = 3\bigl(2 \times 0.6931471806 + \tfrac12 \times 0.6931471806\bigr) = 3 \times 1.7328679514 = 5.1986038543 ,

so $6k_1 = 2.9486038543$ and

k_1 = 0.4914339757 ,

against the true $f'(2) = 0.5$ : the spline recovers the slope at the interior node to $1.7\%$ without being told it.

Now evaluate on $[t_1,t_2] = [2,4]$ , where $h = 2$ and $x = 3$ means $s = \tfrac12$ . The Hermite values at $s = \tfrac12$ are $H_{00} = H_{01} = \tfrac12$ , $H_{10} = \tfrac18$ , $H_{11} = -\tfrac18$ , so

y(3) = \tfrac12(0.6931471806) + \tfrac12(1.3862943611) + 2(0.4914339757)\tfrac18 + 2(0.25)\bigl(-\tfrac18\bigr)

= 0.3465735903 + 0.6931471806 + 0.1228584939 - 0.0625 = 1.1000792648 .

The spline's error is $1.0986122887 - 1.1000792648 = -0.0014669761$ , smaller than the quadratic's by a factor of $38.6$ — on the same three data points, with no new information except the two end slopes. As a check on the construction, the second derivative from the left-hand piece at $t_1 = 2$ is

-\frac{6\Delta_1}{h_0^{2}} + \frac{2k_0+4k_1}{h_0} = -4.1588830834 + (2 + 1.9657359028) = -0.1931471806 ,

and from the right-hand piece

\frac{6\Delta_2}{h_1^{2}} - \frac{4k_1+2k_2}{h_1} = 1.0397207708 - \frac{1.9657359028 + 0.5}{2} = -0.1931471806 :

equal, as $C^{2}$ demands.

Application(How big does a table have to be?)Original

Before every physical calculation that uses a tabulated quantity — an opacity, an equation of state, a collision cross-section, a radiative transfer coefficient, the output of an expensive simulation sampled on a grid — someone must decide how finely to tabulate. The error formula answers that question in one line, and the answer decides the size of the file.

Suppose a function is to be tabulated at spacing $h$ and read back by linear interpolation between the two neighbouring entries. That is the case $n = 1$ of the error formula: on $[x_0, x_0+h]$ the node polynomial $w(x) = (x-x_0)(x-x_0-h)$ has maximum modulus $h^{2}/4$ at the midpoint, so

|f - p| \ \le \ \frac{\max|f''|}{2}\cdot\frac{h^{2}}{4} = \frac{h^{2}}{8}\max|f''| .

For $f = \sin$ on $[0,\pi]$ , where $\max|f''| = 1$ , an accuracy of $10^{-6}$ needs $h \le \sqrt{8 \times 10^{-6}} = 2.828\times10^{-3}$ , that is $\lceil \pi/h \rceil = 1111$ intervals and $1112$ stored values. Evaluating the interpolant on a fine grid near $x = \pi/2$ with $h = \pi/1111$ gives a largest error of $9.995\times10^{-7}$ — the bound is attained to four digits here, because $f''$ is nearly constant across one interval.

Now read the same table back with a cubic through the four nearest entries, evaluating only in the central interval. Then $n = 3$ , and on that central interval $w(x) = (x-x_0)(x-x_1)(x-x_2)(x-x_3)$ has maximum modulus at the midpoint, where the four factors are $\tfrac32 h, \tfrac12 h, -\tfrac12 h, -\tfrac32 h$ , giving $\max|w| = \tfrac{9}{16}h^{4}$ . Hence

|f-p| \ \le \ \frac{\max|f^{(4)}|}{4!}\cdot \frac{9h^{4}}{16} = \frac{3h^{4}}{128}\max|f^{(4)}| ,

and for $\sin$ the same $10^{-6}$ now needs only $h \le (128\times10^{-6}/3)^{1/4} = 8.082\times10^{-2}$ : $39$ intervals and $40$ stored values, with one extra entry beyond each end to supply the four-point stencil for the outermost intervals. Direct computation with $h = 0.0808$ around $x = \pi/2$ gives a largest error of $9.984\times10^{-7}$ , against the predicted bound $9.990\times10^{-7}$ . The table shrinks by a factor of $28$ , at the cost of four multiply-adds instead of one per lookup — which is why almost every production table is read with a cubic rule, and why a smooth spline is used when the reader also needs derivatives, as it does when the table feeds a Newton iteration (Root finding) or a differential-equation solver (Integrators for ODEs).

The same error formula is the common ancestor of the rest of this part. Integrating the interpolant over the nodes gives the quadrature rules and their error terms (Numerical quadrature); differentiating it gives the finite-difference formulas and their truncation errors (Finite-difference approximations); and interpolating by piecewise polynomials on a mesh, exactly as the spline does, gives the shape functions of the finite element method (The finite element method).

Method

Numerical quadrature

Adapted

Add to PDF

The trapezium and Simpson rules with their exact error terms, the expansion of the trapezium error in even powers of the spacing that makes Richardson extrapolation and the Romberg table work, and the Gauss rules, whose n nodes buy exactness to degree 2n-1 because they are the zeros of an orthogonal polynomial.

Prerequisites: Polynomial and spline interpolation, General properties of orthogonal polynomials, Taylor's theorem, Euler-Maclaurin summation, Legendre polynomials

NoteOriginal

What it means. An integral that has no closed form is still a number, and a quadrature rule is a recipe that returns it: sample the integrand at finitely many points and take a weighted sum. Every rule in use is got by integrating an interpolant (Polynomial and spline interpolation), so the interpolation error formula hands over the quadrature error formula, and the whole subject is the question of where to put the nodes. Two answers dominate. Keep the nodes equally spaced, which you must when the data are a table, and use many short panels: that is the trapezium and Simpson rules, and the trapezium rule turns out to have an error series in even powers of the spacing alone, which can be extrapolated away two orders at a time. Or choose the nodes freely, in which case they should be the zeros of an orthogonal polynomial, and $n$ of them give exactness to degree $2n-1$ instead of $n-1$ : twice the return for the same number of function evaluations.

Definition(Quadrature rule, degree of exactness, interpolatory and composite rules)AdaptedAdapted from Kreyszig (1978), §4.11, the process (3) and the uniqueness argument (8)-(10), pp. 277-279

Let $w$ be a non-negative weight function on $(a,b)$ with $\int_a^b w < \infty$ and all moments $\mu_k = \int_a^b x^k w(x)\,dx$ finite, as in General properties of orthogonal polynomials. The unweighted case is $w \equiv 1$ on a bounded interval.

A quadrature rule with $n$ nodes is a choice of distinct points $x_1, \dots, x_n$ in $[a,b]$ and real weights $W_1, \dots, W_n$ , used as

Q(f) = \sum_{i=1}^{n} W_i\, f(x_i) \ \approx\ I(f) = \int_a^b f(x)\, w(x)\,dx .

Both $I$ and $Q$ are linear functionals on $C[a,b]$ , and so is the error functional

E(f) = I(f) - Q(f) .

The rule has degree of exactness $d$ if $E$ vanishes on every polynomial of degree at most $d$ and not on some polynomial of degree $d+1$ . Kreyszig regards $Q$ as a bounded functional on $C[a,b]$ with norm $\sum_i |W_i|$ , which is the point of view his convergence theorem needs.

The rule is interpolatory if its weights are

W_i = \int_a^b \ell_i(x)\, w(x)\,dx , \qquad \ell_i(x) = \prod_{j \neq i} \frac{x-x_j}{x_i - x_j} ,

the Lagrange basis of Polynomial and spline interpolation for these nodes; equivalently, $Q(f) = I(p)$ where $p$ is the interpolant of $f$ at the nodes. For any prescribed nodes these are the unique weights with degree of exactness at least $n-1$ , because exactness on $1, x, \dots, x^{n-1}$ is a linear system whose matrix is a Vandermonde matrix in the distinct nodes, hence invertible.

Newton-Cotes rules are the interpolatory rules on equally spaced nodes with $w \equiv 1$ . A composite rule splits $[a,b]$ into $N$ panels of equal length and applies one short Newton-Cotes rule on each. Writing $h = (b-a)/N$ , $x_j = a + jh$ and $f_j = f(x_j)$ , the two in universal use are the composite trapezium rule

T(h) = h\Bigl(\tfrac12 f_0 + f_1 + f_2 + \cdots + f_{N-1} + \tfrac12 f_N \Bigr)

and, for $N$ even, the composite Simpson rule

S(h) = \frac{h}{3}\Bigl( f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 2f_{N-2} + 4f_{N-1} + f_N \Bigr) .

Method(The two Newton-Cotes rules and their errors)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.4.1, (27.36)-(27.38), pp. 1002-1003, and §27.4.2, (27.39)-(27.40) and the error statement, p. 1004

Take $w \equiv 1$ throughout this statement.

(i) One trapezium panel. If $f \in C^2[c, c+h]$ then there is a $\xi \in (c, c+h)$ with

\int_c^{c+h} f(x)\,dx \;-\; \frac{h}{2}\bigl(f(c) + f(c+h)\bigr) \;=\; -\frac{h^3}{12}\, f''(\xi) .

The rule has degree of exactness $1$ .

(ii) One Simpson panel. If $f \in C^4[c, c+2h]$ then there is a $\xi \in (c, c+2h)$ with

\int_c^{c+2h} f(x)\,dx \;-\; \frac{h}{3}\bigl(f(c) + 4f(c+h) + f(c+2h)\bigr) \;=\; -\frac{h^5}{90}\, f^{(4)}(\xi) .

The rule has degree of exactness $3$ , one more than the degree of the parabola it integrates.

(iii) Composite trapezium. If $f \in C^2[a,b]$ there is an $\eta \in (a,b)$ with

I - T(h) \;=\; -\frac{(b-a)\,h^2}{12}\, f''(\eta) .

(iv) Composite Simpson. If $f \in C^4[a,b]$ and $N$ is even there is an $\eta \in (a,b)$ with

I - S(h) \;=\; -\frac{(b-a)\,h^4}{180}\, f^{(4)}(\eta) .

Riley, Hobson and Bence state (iii) and (iv) with the mean value replaced by an average of the derivative over $[a,b]$ , in their equation (27.38) and the corresponding line for Simpson; the mean-value forms above are the same statements made exact, and are what is proved below.

Idea of the proofOriginal

Each error functional kills polynomials up to some degree, so Taylor's theorem with integral remainder, applied with exactly that many terms, writes the error as a fixed kernel integrated against a derivative of $f$ . The kernel does not depend on $f$ at all; it can be computed once from the rule. If it happens to have one sign — and for these two rules it does — the integral is squeezed between the extreme values of that derivative, and the intermediate value theorem produces the point $\xi$ . The constant in front is then read off by feeding the rule a single monomial.

The extra degree of exactness in Simpson's rule, degree $3$ from a parabola, is the symmetry of the panel: the error functional is even about the midpoint, so it annihilates the odd function $(x-c-h)^3$ for free.

The composite statements follow by summing the panels and averaging: a sum of $N$ values of a continuous function is $N$ times one of its values.

ProofOriginal

Lemma (Peano kernel). Let $E$ be a linear functional on $C^{k}[a,b]$ of the form $E(f) = \int_a^b f\,w - \sum_i W_i f(x_i)$ , and suppose $E$ vanishes on every polynomial of degree at most $k-1$ . Then for $f \in C^k[a,b]$ ,

E(f) = \int_a^b K(t)\, f^{(k)}(t)\,dt, \qquad K(t) = \frac{1}{(k-1)!}\,E_x\bigl[(x-t)_+^{\,k-1}\bigr] ,

where $u_+ = \max(u,0)$ and $E_x$ means $E$ applied in the variable $x$ with $t$ held fixed. If moreover $K$ does not change sign on $[a,b]$ , then there is a $\xi \in (a,b)$ with

E(f) = \frac{f^{(k)}(\xi)}{k!}\; E\bigl(x^{k}\bigr) .

Proof of the lemma. Taylor's theorem with integral remainder (Taylor's theorem) at the left endpoint gives, for every $x \in [a,b]$ ,

f(x) = \sum_{j=0}^{k-1} \frac{f^{(j)}(a)}{j!}(x-a)^{j} + \frac{1}{(k-1)!}\int_a^{x}(x-t)^{k-1}f^{(k)}(t)\,dt = \sum_{j=0}^{k-1}\frac{f^{(j)}(a)}{j!}(x-a)^{j} + \frac{1}{(k-1)!}\int_a^{b}(x-t)_+^{k-1}f^{(k)}(t)\,dt ,

the upper limit being extended to $b$ because the truncated power vanishes for $t > x$ . Apply $E$ . The polynomial part has degree $k-1$ and is annihilated by hypothesis. In the remaining term $E$ is a finite combination of an integral in $x$ and $n$ point evaluations in $x$ , all of which commute with the $t$ -integral, the integrand $(x,t) \mapsto (x-t)_+^{k-1}f^{(k)}(t)$ being bounded and measurable on the square and continuous off the diagonal. That gives the kernel formula.

Now suppose $K \le 0$ throughout (the case $K \ge 0$ is the same with the inequalities reversed), and put $m = \min_{[a,b]} f^{(k)}$ , $M = \max_{[a,b]} f^{(k)}$ , which exist because $f^{(k)}$ is continuous on a compact interval. Multiplying $m \le f^{(k)}(t) \le M$ by $K(t) \le 0$ reverses the inequalities, and integrating gives

M\int_a^b K \ \le\ E(f) \ \le\ m \int_a^b K .

Taking $f(x) = x^{k}$ , whose $k$ -th derivative is the constant $k!$ , in the kernel formula gives $\int_a^b K = E(x^k)/k!$ , which is therefore non-positive and, for the two rules below, non-zero. Hence $E(f)\big/\int_a^b K$ lies between $m$ and $M$ , and since $f^{(k)}$ is continuous it attains that value at some $\xi$ ; the intermediate value theorem gives $\xi$ in the open interval unless $f^{(k)}$ is constant, in which case any interior point serves. $\square$

(i) The trapezium panel. Translate so that $c = 0$ ; the rule is $Q(f) = \tfrac{h}{2}(f(0)+f(h))$ and $E$ annihilates $1$ and $x$ by direct substitution, so take $k=2$ . For $t \in (0,h)$ ,

K(t) = E_x\bigl[(x-t)_+\bigr] = \int_t^h (x-t)\,dx - \frac{h}{2}\Bigl( (0-t)_+ + (h-t)\Bigr) = \frac{(h-t)^2}{2} - \frac{h(h-t)}{2} = -\frac{t\,(h-t)}{2} ,

using $(0-t)_+ = 0$ . This is $\le 0$ on $[0,h]$ and vanishes only at the two endpoints, so the lemma applies. Finally

E(x^2) = \int_0^h x^2\,dx - \frac{h}{2}\bigl(0 + h^2\bigr) = \frac{h^3}{3} - \frac{h^3}{2} = -\frac{h^3}{6} ,

so $E(f) = f''(\xi)\,(-h^3/6)/2! = -h^3 f''(\xi)/12$ . Degree of exactness is exactly $1$ because $E(x^2) \neq 0$ .

(ii) The Simpson panel. Translate so that the panel is $[-h,h]$ with midpoint $0$ ; the rule is $Q(f) = \tfrac{h}{3}(f(-h)+4f(0)+f(h))$ . Substituting $1, x, x^2$ gives $E = 0$ in each case ( $\int_{-h}^{h}x^2 = 2h^3/3$ and $\tfrac{h}{3}(h^2+0+h^2) = 2h^3/3$ ), and $E(x^3) = 0$ because both the integral over a symmetric interval and the symmetric weighted sum vanish for an odd function. So $E$ annihilates $\mathcal{P}_3$ and we take $k=4$ . For $t \in [0,h]$ , the three evaluation points contribute only through $x = h$ , since $(-h-t)_+ = (0-t)_+ = 0$ , and

3!\,K(t) = \int_t^h (x-t)^3\,dx - \frac{h}{3}(h-t)^3 = \frac{(h-t)^4}{4} - \frac{h(h-t)^3}{3} = -\frac{(h-t)^3\,(h+3t)}{12} \ \le\ 0 .

For $t \in [-h,0]$ the same computation with the term $4(0-t)^3 = -4t^3$ retained gives $K(-t) = K(t)$ , so the kernel is even and non-positive throughout. Then

E(x^4) = \int_{-h}^{h}x^4\,dx - \frac{h}{3}\bigl(h^4 + 0 + h^4\bigr) = \frac{2h^5}{5} - \frac{2h^5}{3} = -\frac{4h^5}{15},

and the lemma gives $E(f) = f^{(4)}(\xi)(-4h^5/15)/4! = -h^5 f^{(4)}(\xi)/90$ . Since $E(x^4) \neq 0$ the degree of exactness is exactly $3$ .

(iii) and (iv) Composite rules. Sum part (i) over the $N$ panels $[x_{j}, x_{j+1}]$ : there are points $\xi_j$ with

I - T(h) = -\frac{h^3}{12}\sum_{j=0}^{N-1} f''(\xi_j) = -\frac{h^3}{12}\, N\, \overline{f''}, \qquad \overline{f''} = \frac1N\sum_j f''(\xi_j) .

The average $\overline{f''}$ lies between the smallest and the largest of the $f''(\xi_j)$ , hence between $\min f''$ and $\max f''$ , so by the intermediate value theorem $\overline{f''} = f''(\eta)$ for some $\eta \in (a,b)$ . With $Nh = b-a$ this is (iii). Part (iv) is identical, summing part (ii) over the $N/2$ double panels and using $(N/2)h^5 = (b-a)h^4/2$ , which with the $1/90$ gives $(b-a)h^4/180$ . $\square$

Method(The trapezium error runs in even powers of h)Original

Let $q \ge 1$ and $f \in C^{2q}[a,b]$ , and let $T(h)$ be the composite trapezium rule with $h = (b-a)/N$ . Then

T(h) \;=\; I \;+\; \sum_{p=1}^{q-1} c_{2p}\, h^{2p} \;+\; \mathcal{E}_q(h), \qquad c_{2p} = \frac{B_{2p}}{(2p)!}\Bigl[f^{(2p-1)}(b) - f^{(2p-1)}(a)\Bigr] ,

with the Bernoulli numbers $B_2 = \tfrac16$ , $B_4 = -\tfrac1{30}$ , $B_6 = \tfrac1{42}$ of Euler-Maclaurin summation, and with

\bigl|\mathcal{E}_q(h)\bigr| \ \le\ \frac{2\,|B_{2q}|}{(2q)!}\; h^{2q} \int_a^b \bigl|f^{(2q)}(x)\bigr|\,dx .

Two features matter and neither is obvious from part (iii) of the previous statement. Only even powers of $h$ occur. And every coefficient depends on $f$ only through the two endpoints: if $f$ is periodic with period $b-a$ , or if all its odd derivatives agree at the two ends, every $c_{2p}$ vanishes and the trapezium rule converges faster than any power of $h$ .

Method(Richardson extrapolation and the Romberg table)Original

(i) One extrapolation. Suppose a computable quantity satisfies $A(h) = A + c\,h^{p} + O(h^{r})$ as $h \to 0$ , with $r > p > 0$ and $c$ independent of $h$ . Then

A_{\text{new}}(h) = \frac{2^{p}A(h/2) - A(h)}{2^{p}-1} = A + O(h^{r}) :

the leading error term is cancelled exactly, at the cost of one more evaluation of $A$ .

(ii) Romberg. Let $h_i = (b-a)/2^{i}$ and $R_{i,0} = T(h_i)$ , and define

R_{i,k} = \frac{4^{k} R_{i,k-1} - R_{i-1,k-1}}{4^{k}-1}, \qquad 1 \le k \le i .

If $f \in C^{2k+2}[a,b]$ then $R_{i,k} - I = O(h_i^{2k+2})$ , so each column of the table gains two orders of accuracy over the one to its left, and the diagonal $R_{i,i}$ converges very fast.

(iii) The first column is Simpson. $R_{i,1}$ is exactly the composite Simpson rule with spacing $h_i$ . So Simpson's rule is not a separate idea: it is one Richardson step applied to the trapezium rule.

The table costs nothing beyond the trapezium sums, and those are nested: the nodes of $T(h_{i-1})$ are half the nodes of $T(h_i)$ , so

T(h_i) = \tfrac12 T(h_{i-1}) + h_i \sum_{j \text{ odd},\, 0<j<2^{i}} f(a + j h_i) ,

and the whole table through row $i$ costs $2^{i}+1$ evaluations of $f$ in all.

Idea of the proofOriginal

The even-power expansion is the Euler-Maclaurin formula of Euler-Maclaurin summation, which compares a sum over integers with an integral, rescaled to spacing $h$ . Each derivative in that formula picks up one factor of $h$ from the chain rule and the whole formula is multiplied by one more factor of $h$ to turn a sum into a trapezium rule; since only odd-order derivatives occur in Euler-Maclaurin, only even powers of $h$ survive. The interior values of $f$ cancel in telescoping pairs, which is why only the endpoints appear.

Given such an expansion, extrapolation is just solving two equations for two unknowns: $A(h)$ and $A(h/2)$ both estimate $A$ , with errors in the known ratio $2^{p}$ , so the combination that eliminates the error term is forced. Doing it repeatedly, each time against the next surviving power, is the Romberg table.

ProofOriginal

The even-power expansion. Put $g(s) = f(a + sh)$ for $s \in [0,N]$ , so $g \in C^{2q}[0,N]$ and $g^{(m)}(s) = h^{m}f^{(m)}(a+sh)$ . Euler-Maclaurin summation with parameter $q$ applied to $g$ on the integer interval $[0,N]$ gives

\sum_{j=0}^{N} g(j) = \int_0^{N} g(s)\,ds + \frac{g(0)+g(N)}{2} + \sum_{p=1}^{q}\frac{B_{2p}}{(2p)!}\bigl[g^{(2p-1)}(N)-g^{(2p-1)}(0)\bigr] + R_q , \qquad |R_q| \le \frac{|B_{2q}|}{(2q)!}\int_0^{N}\bigl|g^{(2q)}\bigr| .

Multiply by $h$ and read each piece back in terms of $f$ . On the left, $h\sum_{j=0}^{N}g(j) = T(h) + \tfrac{h}{2}(f_0+f_N)$ , because $T$ halves the two end values. On the right, the substitution $x = a+sh$ gives $h\int_0^N g = \int_a^b f$ , the half-sum term is $\tfrac{h}{2}(f_0+f_N)$ and cancels, and $h\,g^{(2p-1)}(s) = h^{2p}f^{(2p-1)}(a+sh)$ , which at $s = N$ and $s=0$ gives $h^{2p}f^{(2p-1)}(b)$ and $h^{2p}f^{(2p-1)}(a)$ . Likewise $h|R_q| \le h \cdot \tfrac{|B_{2q}|}{(2q)!}h^{2q}\cdot \tfrac1h\int_a^b|f^{(2q)}|$ . Hence

T(h) = I + \sum_{p=1}^{q}c_{2p}h^{2p} + hR_q, \qquad |hR_q| \le \frac{|B_{2q}|}{(2q)!}h^{2q}\int_a^b\bigl|f^{(2q)}\bigr| .

Moving the last term of the sum into the remainder, $\mathcal{E}_q(h) = c_{2q}h^{2q} + hR_q$ , and

|c_{2q}| = \frac{|B_{2q}|}{(2q)!}\Bigl|\int_a^b f^{(2q)}\Bigr| \le \frac{|B_{2q}|}{(2q)!}\int_a^b\bigl|f^{(2q)}\bigr|

by the fundamental theorem of calculus, which doubles the bound. This proves the statement, and the case $q=1$ reproduces the order $h^2$ already known from the mean-value form.

(i) Richardson. Write $A(h) = A + ch^{p} + \rho(h)$ with $|\rho(h)| \le Ch^{r}$ . Then

\frac{2^{p}A(h/2)-A(h)}{2^{p}-1} = \frac{2^{p}A + 2^{p}c (h/2)^{p} - A - ch^{p}}{2^{p}-1} + \frac{2^{p}\rho(h/2)-\rho(h)}{2^{p}-1} = A + 0 + O(h^{r}) ,

since $2^{p}(h/2)^{p} = h^{p}$ makes the two $c$ -terms cancel, and the remainder is bounded by $(2^{p}2^{-r}+1)Ch^{r}/(2^{p}-1)$ .

(ii) Romberg. Induct on $k$ . For $k=0$ the previous statement with $q = k+1 = 1$ gives $R_{i,0} - I = O(h_i^{2})$ . Suppose that for some $k \ge 1$ , and for $f \in C^{2k+2}$ ,

R_{i,k-1} = I + \sum_{p=k}^{m} d_{p}h_i^{2p} + O(h_i^{2m+2})

for coefficients $d_p$ independent of $i$ — which is true at $k=1$ by the even-power expansion with $q$ as large as the smoothness allows. Applying (i) with $p$ replaced by $2k$ (so $2^{p}$ becomes $4^{k}$ ) to the pair $(R_{i-1,k-1}, R_{i,k-1})$ , where $h_i = h_{i-1}/2$ , cancels the $h^{2k}$ term and leaves an expansion starting at $h_i^{2k+2}$ with coefficients again independent of $i$ . That is the inductive step, and it needs $f \in C^{2k+2}$ for the expansion to reach that far.

(iii) Simpson. With $h = h_i$ and $2h = h_{i-1}$ ,

R_{i,1} = \frac{4T(h)-T(2h)}{3} .

Write out both sums. $T(h) = h(\tfrac12 f_0 + f_1 + \cdots + f_{N-1} + \tfrac12 f_N)$ with $N$ even, and $T(2h) = 2h(\tfrac12 f_0 + f_2 + f_4 + \cdots + f_{N-2} + \tfrac12 f_N)$ , which uses only the even-indexed nodes. Then

\frac{4T(h)-T(2h)}{3} = \frac{h}{3}\Bigl[(2-1)f_0 + 4\!\!\sum_{j \text{ odd}}\!\! f_j + (4-2)\!\!\sum_{\substack{j \text{ even}\\ 0<j<N}}\!\! f_j + (2-1)f_N\Bigr] = \frac{h}{3}\Bigl[f_0 + 4\!\!\sum_{j \text{ odd}}\!\! f_j + 2\!\!\sum_{\substack{j \text{ even}\\ 0<j<N}}\!\! f_j + f_N\Bigr],

which is $S(h)$ . $\square$

Method(Gauss rules)AdaptedAdapted from Riley, Hobson & Bence (2006), §27.4.3, the scheme and the exactness to degree 2n-1, pp. 1005-1007

Let $w$ be a weight on $(a,b)$ as above, let $(\pi_n)$ be its monic orthogonal family (General properties of orthogonal polynomials), and let $x_1 < x_2 < \cdots < x_n$ be the zeros of $\pi_n$ , which that unit shows to be real, simple and interior to $(a,b)$ . Let $W_i = \int_a^b \ell_i w$ be the interpolatory weights for these nodes. Then:

Exactness. $E(f) = 0$ for every polynomial $f$ of degree at most $2n-1$ .
And no more. $E(\pi_n^2) = \int_a^b \pi_n^2\, w > 0$ , and $\pi_n^2$ has degree $2n$ ; so the degree of exactness is exactly $2n-1$ . No $n$ -node rule whatever can do better.
Positive weights. $W_i = \int_a^b \ell_i^2\, w > 0$ for every $i$ , and $\sum_i W_i = \mu_0 = \int_a^b w$ .
Uniqueness. It is the only $n$ -node rule with degree of exactness $2n-1$ .
Error. If $f \in C^{2n}[a,b]$ there is a $\xi \in (a,b)$ with

E(f) = \frac{f^{(2n)}(\xi)}{(2n)!}\int_a^b \pi_n(x)^2\,w(x)\,dx .

Gauss-Legendre. For $w \equiv 1$ on $[-1,1]$ the nodes are the zeros of the Legendre polynomial $P_n$ and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.4.3, p. 1006) give the weights as

W_i = \frac{2}{(1-x_i^{2})\,\bigl|P_n'(x_i)\bigr|^{2}} ,

tabulating them to ten decimals for $n$ up to $20$ in their Table 27.8, p. 1008. Since $P_n$ has leading coefficient $(2n)!/(2^{n}(n!)^{2})$ and $\int_{-1}^{1}P_n^2 = 2/(2n+1)$ (Legendre polynomials), the error constant is

\int_{-1}^{1}\pi_n^{2}\,dx = \frac{2}{2n+1}\left(\frac{2^{n}(n!)^{2}}{(2n)!}\right)^{2} = \frac{2^{2n+1}(n!)^{4}}{(2n+1)\bigl((2n)!\bigr)^{2}} ,

which is $\tfrac{2}{3}, \tfrac{8}{45}, \tfrac{8}{175}, \tfrac{128}{11025}$ for $n = 1,2,3,4$ . A general interval is reached by the affine substitution $x = \tfrac{a+b}{2}+\tfrac{b-a}{2}z$ , which multiplies the integral by $(b-a)/2$ .

Other weights absorb whatever is awkward about the integrand, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.4.3, pp. 1005, 1009) list the standard ones: $e^{-x}$ on $(0,\infty)$ for Gauss-Laguerre, $e^{-x^{2}}$ on $\mathbb{R}$ for Gauss-Hermite, and $(1-x^{2})^{\pm 1/2}$ on $(-1,1)$ for the two Gauss-Chebyshev families, whose nodes and weights they give in closed form, $x_i = \cos\frac{(i-\frac12)\pi}{n}$ with all weights equal to $\pi/n$ in the first-kind case.

Method(Convergence as the number of nodes grows)AdaptedAdapted from Kreyszig (1978), §4.11, Theorems 4.11-3 (Polya) and 4.11-4 (Steklov), pp. 279-280

Let $w \equiv 1$ on $[a,b]$ and let $Q_n$ , $n = 1,2,\dots$ , be a sequence of rules, $Q_n$ having $n$ nodes and $n$ weights $W_i^{(n)}$ , each exact on polynomials of degree at most $n-1$ .

Polya. $Q_n(f) \to \int_a^b f$ for every $f \in C[a,b]$ if and only if there is a constant $c$ with

\sum_{i=1}^{n}\bigl|W_i^{(n)}\bigr| \le c \qquad \text{for all } n .

Steklov. If every weight is non-negative, this bound is automatic, since taking $f \equiv 1$ gives $\sum_i W_i^{(n)} = b-a$ ; so such a sequence converges for every continuous $f$ .

By part 3 of the previous statement the Gauss rules have positive weights, so Steklov's theorem applies to them: Gauss quadrature converges for every continuous integrand as $n$ grows. The equally spaced Newton-Cotes rules do not have this property. Computing their weights exactly, all are positive up to $8$ nodes; at $9$ nodes negative weights appear for the first time and $\sum_i|W_i^{(n)}|/(b-a) = 6857/4725 = 1.4512$ , at $11$ nodes it is $152921/49896 = 3.0648$ , and the quantity grows without bound thereafter. Raising the order of a single Newton-Cotes rule is therefore not a convergent strategy, which is the reason composite rules exist.

Idea of the proofOriginal

Everything about Gauss rules comes from one division. A polynomial $f$ of degree $2n-1$ can be divided by $\pi_n$ to give $f = q\pi_n + r$ with $q$ and $r$ of degree below $n$ . The rule does not see the quotient term at all, because every node is a root of $\pi_n$ ; the integral does not see it either, because $\pi_n$ is orthogonal to everything of lower degree, which is exactly what $q$ is. What is left on both sides is the remainder $r$ , of degree below $n$ , on which any interpolatory rule is exact. The two sides therefore agree — and the argument plainly fails at degree $2n$ , where $q$ can be constant and $\int q \pi_n w$ need not vanish.

Positivity of the weights is the same trick run on $\ell_i^2$ , which has degree $2n-2$ and so is integrated exactly, while the rule evaluates it as $W_i$ .

The error formula replaces the interpolant at $n$ points by the Hermite interpolant that also matches $n$ derivatives: it has degree $2n-1$ , so the rule integrates it exactly, and the difference $f - H$ carries the square of the node polynomial, which has a fixed sign.

Polya's theorem is the uniform boundedness principle in disguise: the rules are bounded functionals on $C[a,b]$ , they already converge on the polynomials, which are dense, and a family of functionals that converges on a dense set converges everywhere precisely when its norms are bounded.

ProofOriginal

1. Exactness. Let $\deg f \le 2n-1$ . Polynomial division by the monic $\pi_n$ gives $f = q\,\pi_n + r$ with $\deg q \le n-1$ and $\deg r \le n-1$ . Then

I(f) = \int_a^b q\,\pi_n\, w + \int_a^b r\,w = 0 + I(r) ,

the first integral vanishing because $\pi_n$ is orthogonal to every polynomial of degree below $n$ (General properties of orthogonal polynomials) and $\deg q \le n-1$ . On the other side, $\pi_n(x_i) = 0$ for every node, so

Q(f) = \sum_i W_i\bigl(q(x_i)\pi_n(x_i) + r(x_i)\bigr) = Q(r) .

Finally $I(r) = Q(r)$ because the rule is interpolatory and $\deg r \le n-1$ . Hence $E(f) = 0$ .

2. The degree is not higher. Take $f = \pi_n^2$ , of degree $2n$ . Every node is a root, so $Q(f) = 0$ , while $I(f) = \int \pi_n^2 w > 0$ because $w \ge 0$ is a weight and $\pi_n^2 \ge 0$ is not identically zero. For an arbitrary $n$ -node rule with nodes $y_i$ , the same argument with $\prod_i (x-y_i)^2$ shows no $n$ -node rule is exact on all of $\mathcal{P}_{2n}$ .

3. Positivity. $\ell_i^2$ has degree $2n-2 \le 2n-1$ , so by part 1,

\int_a^b \ell_i^2\,w = \sum_j W_j\,\ell_i(x_j)^2 = W_i ,

since $\ell_i(x_j) = \delta_{ij}$ . The left side is strictly positive: $\ell_i^2 w \ge 0$ and $\ell_i^2$ is a non-zero polynomial, so it is bounded below by a positive constant on some subinterval where $w$ is not almost everywhere zero. Taking $f \equiv 1$ in part 1 gives $\sum_i W_i = \mu_0$ .

4. Uniqueness. Suppose nodes $y_1,\dots,y_n$ and weights $V_1,\dots,V_n$ give a rule exact on $\mathcal{P}_{2n-1}$ . Put $u(x) = \prod_i (x-y_i)$ , monic of degree $n$ . For any polynomial $q$ with $\deg q \le n-1$ the product $qu$ has degree at most $2n-1$ , so

\int_a^b q\,u\,w = \sum_i V_i\, q(y_i)\,u(y_i) = 0 ,

each $u(y_i)$ being zero. So $u$ is a monic polynomial of degree $n$ orthogonal to $\mathcal{P}_{n-1}$ , and the monic orthogonal polynomial with that property is unique, so $u = \pi_n$ and the $y_i$ are the Gauss nodes. The weights are then forced, being the unique ones exact on $\mathcal{P}_{n-1}$ .

5. The error formula. Let $f \in C^{2n}[a,b]$ and let $H$ be the Hermite interpolant of $f$ at the nodes: the unique polynomial of degree at most $2n-1$ with $H(x_i) = f(x_i)$ and $H'(x_i) = f'(x_i)$ for every $i$ . (Uniqueness: a difference of two such polynomials has a double root at each of the $n$ nodes, hence is divisible by $\pi_n^2$ , of degree $2n$ ; a polynomial of degree at most $2n-1$ divisible by one of degree $2n$ is zero. Existence then follows because the $2n$ interpolation conditions are a square linear system on the $2n$ -dimensional space $\mathcal{P}_{2n-1}$ , and a square system with trivial kernel is solvable.)

Pointwise remainder. Fix $x$ in $[a,b]$ that is not a node, and set $\lambda = \bigl(f(x)-H(x)\bigr)/\pi_n(x)^2$ and $g(s) = f(s)-H(s)-\lambda\,\pi_n(s)^2$ . Then $g$ vanishes at the $n$ nodes and at $x$ , that is at $n+1$ distinct points, so by Rolle's theorem $g'$ vanishes at $n$ interior points of the gaps; and $g'$ also vanishes at each of the $n$ nodes, because $f'=H'$ there and $\pi_n^2$ has a double root. These $2n$ zeros of $g'$ are distinct, the nodes not being among the interior points just produced. Applying Rolle $2n-1$ more times, $g^{(2n)}$ has a zero $\zeta$ in $(a,b)$ . Since $H^{(2n)} \equiv 0$ and $\pi_n^2$ is monic of degree $2n$ , so that $(\pi_n^2)^{(2n)} \equiv (2n)!$ , this gives $\lambda = f^{(2n)}(\zeta)/(2n)!$ . Hence for every $x$ in $[a,b]$ , nodes included,

f(x) - H(x) = \frac{f^{(2n)}(\zeta_x)}{(2n)!}\,\pi_n(x)^2 .

Integrate. Let $m$ and $M$ be the minimum and maximum of $f^{(2n)}/(2n)!$ on $[a,b]$ . Then $m\,\pi_n(x)^2 \le f(x)-H(x) \le M\,\pi_n(x)^2$ for every $x$ . Multiply by $w \ge 0$ and integrate:

m \int_a^b \pi_n^2 w \ \le\ \int_a^b (f-H)\,w \ \le\ M\int_a^b \pi_n^2 w .

By part 1 the rule is exact on $H$ , and $H$ agrees with $f$ at the nodes, so $\int_a^b H w = Q(H) = Q(f)$ and the middle quantity is exactly $E(f)$ . Dividing by the positive number $\int \pi_n^2 w$ puts $E(f)/\int\pi_n^2 w$ between $m$ and $M$ ; the continuous function $f^{(2n)}/(2n)!$ attains that value at some $\xi$ , interior unless it is constant. $\square$

Polya and Steklov. These are Theorems 4.11-3 and 4.11-4 of Kreyszig (Kreyszig (1978), §4.11, pp. 279-280); his proof identifies $\|Q_n\| = \sum_i|W_i^{(n)}|$ , notes that the polynomials are dense in $C[a,b]$ by the Weierstrass approximation theorem and that convergence already holds on them by exactness, and then appeals to the characterisation of weak* convergence for a sequence of functionals on a Banach space, which requires exactly a uniform bound on the norms (Uniform boundedness and open mapping theorems, Weak convergence).

Remark(Where the hypotheses are used)Original

$f \in C^{2}$ for the trapezium rule, $C^{4}$ for Simpson, $C^{2n}$ for the $n$ -point Gauss rule. These enter at exactly one place each: the Peano lemma needs $f^{(k)}$ continuous to squeeze the kernel integral and to apply the intermediate value theorem, and the Gauss proof needs $f^{(2n)}$ continuous for the same reason. Drop the smoothness and the formulas say nothing. For $f(x) = \sqrt{x}$ on $[0,1]$ , whose second derivative is unbounded at the origin, the composite trapezium error is not $O(h^2)$ but $O(h^{3/2})$ : at $N = 4, 8, 16, 32, 64, 128, 256$ panels the computed errors are $2.3384\times10^{-2}$ , $8.5364\times10^{-3}$ , $3.0855\times10^{-3}$ , $1.1077\times10^{-3}$ , $3.9586\times10^{-4}$ , $1.4101\times10^{-4}$ , $5.0118\times10^{-5}$ , falling by factors $0.365, 0.361, 0.359, 0.357, 0.356, 0.355$ , which are creeping down to $2^{-3/2} = 0.3536$ and not to $2^{-2} = 0.25$ . Romberg extrapolation on such an integrand is worse than useless, because it extrapolates a power that is not there.
The kernel must not change sign. That is what converts a bound into an equality with a definite point $\xi$ . It holds for the trapezium and Simpson panels, as computed, and for every Gauss rule, because there the role of the kernel is played by $\pi_n^2 \ge 0$ . It fails for Newton-Cotes rules of high order, which is another face of their negative weights.
Distinct nodes. Used for the Lagrange basis and so for the very existence of interpolatory weights; the Gauss nodes are distinct because the zeros of an orthogonal polynomial are simple.
$w \ge 0$ . Used three times: for $\int\pi_n^2 w > 0$ , for the positivity of the weights, and in the last step of the error proof, where multiplying an inequality by $w$ must not reverse it. A signed weight destroys all three.
Exactness on $\mathcal{P}_{n-1}$ in Polya's theorem is what makes the rules converge on the dense set of polynomials; the uniform bound on $\sum|W_i|$ is what transfers that to all of $C[a,b]$ . Neither alone suffices, and Kreyszig's statement is an equivalence, so the bound is not merely convenient.
Periodicity is a hypothesis too. The final remark of the even-power expansion is not an approximation but an identity: for a smooth periodic integrand over a full period every $c_{2p}$ vanishes, and the plain trapezium rule is spectrally accurate. This is why Fourier and spectral codes (Spectral methods) integrate with the trapezium rule and nothing fancier, and why The discrete Fourier transform and the FFT computes Fourier coefficients by an equally weighted sum.

Method(Choosing and running a quadrature)Original

Decide whether you may choose the nodes. If the integrand is a table, you may not: use a composite rule on the nodes you have, Simpson if their number is odd (an even number of panels), the trapezium rule otherwise. If you can evaluate $f$ anywhere, use a Gauss rule or Romberg.
Check the smoothness first, not last. Count how many continuous derivatives the integrand has on the closed interval. An endpoint singularity (an inverse square root, a logarithm) invalidates every error formula here; remove it by substitution, by splitting off an analytically integrable piece, or by choosing a Gauss weight that contains it, as Gauss-Chebyshev contains $(1-x^2)^{-1/2}$ .
For a smooth integrand on a finite interval, use Gauss-Legendre. Map to $[-1,1]$ by $x = \tfrac{a+b}{2}+\tfrac{b-a}{2}z$ , take the nodes and weights from a table such as Table 27.8 (Riley, Hobson & Bence (2006), p. 1008), and multiply the result by $(b-a)/2$ . Estimate the error by comparing the $n$ -point and $2n$ -point rules; the $n$ -point error is well approximated by their difference, since the $2n$ -point value is far more accurate.
For a tabulated or repeatedly refined integrand, build the Romberg table. Compute $T(h_0)$ with one panel; then double the number of panels, reusing all previous values through $T(h_i) = \tfrac12T(h_{i-1}) + h_i\sum_{j\text{ odd}}f(a+jh_i)$ ; then fill the row by $R_{i,k} = (4^{k}R_{i,k-1}-R_{i-1,k-1})/(4^{k}-1)$ . Stop when $|R_{i,i}-R_{i-1,i-1}|$ is below the tolerance, and report $R_{i,i}$ .
Never extrapolate further than the smoothness allows. Column $k$ of the Romberg table presupposes $f \in C^{2k+2}$ . If the diagonal stops improving, or starts getting worse, the expansion has run out: stop at the last column that was still improving.
Watch the rounding floor. Each column of the Romberg table amplifies rounding error by roughly $(4^{k}+1)/(4^{k}-1) \le 5/3$ , so the floor is only a few times the floor of $T$ itself, about $n\varepsilon\max|f|$ for $n$ evaluations with unit roundoff $\varepsilon$ (Floating-point arithmetic and error). Refining beyond that point buys nothing: the error stops falling and begins to rise.
In several dimensions, do not use a product rule beyond three or four dimensions. A product of $n$ -point rules costs $n^{d}$ evaluations. Beyond that use Monte Carlo, whose error is $O(M^{-1/2})$ independently of the dimension (Monte Carlo methods).

Example(One integral, four rules, with the predicted error constants checked)Original

Take

I = \int_0^1 e^{x}\,dx = e - 1 = 1.718281828459045 ,

an integral whose exact value is known, so that every error below is an error and not an estimate. All derivatives of $f(x)=e^{x}$ equal $e^{x}$ , which makes the predicted constants explicit.

Trapezium rule. With $N$ panels and $h = 1/N$ , the predicted error is $I - T(h) = -h^2 e^{\eta}/12$ for some $\eta \in (0,1)$ , that is $T(h) - I = h^{2}e^{\eta}/12$ . The computed values are

$N$	$T(h)$	$T(h)-I$	$12\,(T(h)-I)/h^{2}$	$\eta$
1	1.8591409142	$1.4086\times10^{-1}$	1.690309	0.524911
2	1.7539310925	$3.5649\times10^{-2}$	1.711165	0.537174
4	1.7272219046	$8.9401\times10^{-3}$	1.716495	0.540284
8	1.7205185922	$2.2368\times10^{-3}$	1.717835	0.541065

Every $\eta$ lies in $(0,1)$ , as the theorem requires, and the fourth column settles on $e-1 = 1.718282$ , the value forced by the even-power expansion: its leading coefficient is $c_2 = \tfrac{1}{12}[f'(1)-f'(0)] = (e-1)/12$ , so $12(T-I)/h^2 \to e-1$ and $\eta \to \ln(e-1) = 0.541325$ .

Simpson's rule. The prediction is $S(h)-I = h^{4}e^{\eta}/180$ :

$N$	$S(h)$	$S(h)-I$	$180\,(S(h)-I)/h^{4}$
2	1.7188611519	$5.7932\times10^{-4}$	1.668451
4	1.7183188419	$3.7013\times10^{-5}$	1.705580
8	1.7182841547	$2.3262\times10^{-6}$	1.715091

again approaching $e-1$ . Each halving of $h$ divides the error by about $16$ , as an $h^4$ rule must.

Romberg. Build the table from $T(1), T(\tfrac12), \dots, T(\tfrac1{16})$ , which together use $17$ evaluations of $f$ because the nodes are nested. The entries, and in brackets their errors:

$N$	$R_{i,0}$	$R_{i,1}$	$R_{i,2}$	$R_{i,3}$	$R_{i,4}$
1	1.859140914
2	1.753931092	1.718861152
4	1.727221905	1.718318842	1.718282688
8	1.720518592	1.718284155	1.718281842	1.718281829
16	1.718841129	1.718281974	1.718281829	1.7182818285	1.7182818285

with diagonal errors

1.41\times10^{-1},\quad 5.79\times10^{-4},\quad 8.60\times10^{-7},\quad 3.36\times10^{-10},\quad 3.31\times10^{-14} .

Two checks on the construction. The second column reproduces the Simpson values above exactly ( $1.718861152$ and $1.718318842$ are $S(\tfrac12)$ and $S(\tfrac14)$ ), as part (iii) of the Richardson statement says it must. And the errors in each column fall by the predicted factors: column $0$ by $4$ , column $1$ by $16$ , column $2$ by $64$ . Seventeen evaluations have produced thirteen correct digits, where the plain trapezium rule with the same $17$ evaluations gives four.

Gauss-Legendre. Substituting $x = (z+1)/2$ turns $I$ into $\tfrac12\int_{-1}^{1}e^{(z+1)/2}dz$ . The two-point rule has $z = \pm 1/\sqrt3$ , both weights $1$ , so the nodes in $x$ are $0.211324865$ and $0.788675135$ with weights $\tfrac12$ each:

G_2 = \tfrac12\bigl(1.235313601 + 2.200479155\bigr) = 1.717896378, \qquad I - G_2 = 3.8545\times10^{-4} .

The three-point rule has $z = 0, \pm\sqrt{3/5}$ with weights $\tfrac89, \tfrac59$ , so $x = 0.112701665,\, 0.5,\, 0.887298335$ with weights $\tfrac59\cdot\tfrac12 = 0.277777778$ , $\tfrac89\cdot\tfrac12 = 0.444444444$ , $0.277777778$ :

G_3 = 0.277777778(1.119297958) + 0.444444444(1.648721271) + 0.277777778(2.428559625) = 1.718281004,

I - G_3 = 8.2409\times10^{-7} .

Check these against the error formula. On $[0,1]$ the constant is $\tfrac12 \cdot 2^{-2n}\int_{-1}^{1}\pi_n^2\,dz\,/(2n)!$ , the $2^{-2n}$ coming from differentiating $e^{(z+1)/2}$ in $z$ and the $\tfrac12$ from the Jacobian. For $n=2$ , with $\int\pi_2^2 = 8/45$ , this is $\tfrac12\cdot\tfrac1{16}\cdot\tfrac{8/45}{24} = 2.31481\times10^{-4}$ , and $(I-G_2)$ divided by it is $1.665146$ , so $\xi = \ln 1.665146 = 0.509913 \in (0,1)$ . For $n=3$ , with $\int\pi_3^2 = 8/175$ , the constant is $\tfrac12\cdot\tfrac1{64}\cdot\tfrac{8/175}{720} = 4.96032\times10^{-7}$ , giving $1.661358$ and $\xi = 0.507636 \in (0,1)$ . Both error formulas are confirmed to six figures.

The scoreboard. Counting evaluations of $f$ : Gauss with $3$ points reaches $8.2\times10^{-7}$ ; Simpson needs $9$ points to reach $2.3\times10^{-6}$ , and the trapezium rule would need $h^{2}(e-1)/12 \le 8.24\times10^{-7}$ , that is $h \le 2.399\times10^{-3}$ and about $418$ points, to match the three-point Gauss rule. Romberg, which pays the trapezium rule's evaluations but extrapolates them, reaches $3\times10^{-14}$ with $17$ .

Application(Matrix elements, cross-sections and the finite element method)Original

Almost every integral a physicist actually computes is a quadrature, and the choice of rule is dictated by the structure of the integrand rather than by taste.

Angular integrals in scattering and radiative transfer. A phase function or a scattering amplitude is integrated over $\mu = \cos\theta \in [-1,1]$ against a weight of $1$ . The integrand is usually well approximated by a polynomial of modest degree in $\mu$ — that is precisely what a truncated Legendre expansion says — so Gauss-Legendre with the number of nodes set by the retained order is exact, not approximate. A radiative-transfer code that keeps $L$ Legendre moments uses an $n$ -point Gauss-Legendre angular grid with $2n-1 \ge 2L$ , and the discrete ordinates of such codes are nothing but Gauss nodes. If instead the integrand carries $(1-\mu^{2})^{-1/2}$ , as grazing-incidence problems do, Gauss-Chebyshev puts the singularity in the weight and the quadrature never evaluates it.

Matrix elements in quantum mechanics and quantum chemistry. An expectation value $\int \psi^{*}\hat{O}\psi$ over an infinite range with a Gaussian or exponential envelope is exactly what Gauss-Hermite and Gauss-Laguerre are built for: the envelope is the weight, so it is handled analytically, and the rule needs only to resolve the polynomial prefactor. A harmonic-oscillator matrix element between states $n$ and $m$ has a prefactor of degree $n+m$ in the oscillator coordinate, so $\lceil (n+m+1)/2 \rceil$ Gauss-Hermite nodes evaluate it exactly, in arithmetic with no truncation error at all.

Element integrals in the finite element method. Assembling a stiffness matrix (The finite element method) means evaluating one small integral per element, millions of times. The integrands are products of polynomial shape functions and a Jacobian, of known degree, so a Gauss rule of just sufficient order is used and the quadrature is exact to rounding. This is why Gauss rules are hard-wired into finite element codes: the alternative, refining a composite rule until it converges, would multiply the cost of assembly by a large factor and give a worse answer.

Periodic integrands. Orbit averages, time averages over a cycle, and the Fourier coefficients of a periodic signal are integrals over a full period of a smooth periodic function. There the even-power expansion collapses entirely and the equally spaced trapezium rule is the optimal rule, converging faster than any power of the spacing. Using Simpson's rule instead — as is often done by reflex — is strictly worse, because its unequal weights break the symmetry that produced the collapse.

A warning from the error constants. Every statement above about accuracy multiplies a power of $h$ by a high derivative of the integrand. A resonance of width $\Gamma$ in a cross-section has $f^{(2n)} \sim \Gamma^{-2n}$ near the peak, so the error constant is astronomically large until $h \lesssim \Gamma$ . No amount of extrapolation repairs this; the integral must be split at the resonance and the narrow piece resolved on its own scale. The error formula is useful exactly because it tells you this before the computation rather than after.

Method

Solving linear systems numerically

Adapted

Add to PDF

Gaussian elimination as an LU factorization with row interchanges, its cubic cost and why the interchanges are not optional; the Jacobi and Gauss-Seidel iterations and the spectral radius that decides whether they converge; and the conjugate gradient method, which minimises the energy norm of the error over a Krylov subspace and converges at a rate set by the square root of the condition number.

Prerequisites: Linear systems and rank, Matrix norms and condition numbers, Floating-point arithmetic and error, Jordan normal form, The contraction mapping theorem, Chebyshev polynomials

NoteOriginal

What it means. Solving $Ax = b$ is the inner loop of almost every computation in applied science, and there are exactly two strategies. A direct method does a fixed amount of arithmetic — about $n^3/3$ multiply-add pairs — and returns the answer; its only enemy is rounding, and the whole content of the subject is the row interchanges that keep rounding under control. An iterative method produces a sequence converging to the answer; it never touches a matrix entry that is zero, so it is the only option when $n$ is a million and $A$ is sparse, and its content is the question of how fast the sequence converges. This unit gives the direct method, the two classical iterations, and the one iteration that is not classical and not stationary — the conjugate gradient method — which in a precise sense extracts the best possible answer from the work it has done.

Definition(LU factorization, partial pivoting, splitting, residual)AdaptedAdapted from Riley, Hobson & Bence (2006), §8.18.2, Eq. (8.125)-(8.126), pp. 296-297, and §27.3.1, p. 995; Kreyszig (1978), §5.2, the splitting A = B - G and the iteration (9), p. 309

Let $A$ be a real $n \times n$ matrix and $b \in \mathbb{R}^n$ .

An LU factorization of $A$ is $A = LU$ with $L$ unit lower triangular ( $\ell_{ii} = 1$ , $\ell_{ij} = 0$ for $j > i$ ) and $U$ upper triangular. Given one, $Ax = b$ is solved by forward substitution $Ly = b$ followed by back substitution $Ux = y$ , each a triangular system solved one entry at a time.

A permutation matrix $P$ is the identity with its rows permuted; $PA$ is $A$ with its rows permuted the same way. Gaussian elimination with partial pivoting is the algorithm that, at step $k$ , exchanges row $k$ with the row below it carrying the largest entry of column $k$ in modulus, then eliminates below the diagonal. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.3.1, p. 995) call this "elimination with interchange" and give the reason: it stops a small number arriving on the diagonal and being divided by.

A splitting of $A$ is a decomposition $A = M - N$ with $M$ invertible and cheap to invert. It defines the stationary iteration

M x^{m+1} = N x^{m} + b, \qquad \text{equivalently} \qquad x^{m+1} = C x^{m} + d, \quad C = M^{-1}N = I - M^{-1}A, \quad d = M^{-1}b ,

with iteration matrix $C$ . Writing $A = D - L - U$ with $D$ the diagonal, $-L$ the strictly lower part and $-U$ the strictly upper part, the two classical choices are

\textbf{Jacobi: } M = D, \qquad\qquad \textbf{Gauss-Seidel: } M = D - L .

Stakgold and Holst write the same thing as $x^{m+1} = x^m + B(b - Ax^m)$ with $B = M^{-1}$ an approximate inverse of $A$ , so that $C = I - BA$ is the error propagator; Jacobi is $B = D^{-1}$ and Gauss-Seidel is $B = (D-L)^{-1}$ .

For an approximation $\hat{x}$ the residual is $r = b - A\hat{x}$ and the error is $x - \hat{x} = A^{-1}r$ . They are different things, and $\kappa(A) = \|A\|\,\|A^{-1}\|$ of Matrix norms and condition numbers is exactly the factor by which they may differ in relative terms.

Method(Gaussian elimination: existence, pivots, and cost)Original

Let $A$ be $n \times n$ .

(i) Without pivoting. $A$ has an LU factorization with $L$ unit lower triangular if and only if the elimination never meets a zero pivot; a sufficient condition is that every leading principal submatrix $A_{1:k,1:k}$ is nonsingular, $k = 1,\dots,n$ . When it exists and $A$ is nonsingular — equivalently, when $U$ has no zero on its diagonal — the factorization is unique. (For singular $A$ it need not be: $\bigl(\begin{smallmatrix}0&0\\0&0\end{smallmatrix}\bigr) = LU$ for every unit lower triangular $L$ and $U = 0$ .)

(ii) With pivoting. For every nonsingular $A$ there is a permutation matrix $P$ with

PA = LU ,

$L$ unit lower triangular and $U$ upper triangular and nonsingular. Partial pivoting produces such a $P$ , and the multipliers it produces satisfy $|\ell_{ij}| \le 1$ for all $i > j$ .

(iii) Cost. Factorization costs

\sum_{k=1}^{n-1}(n-k)^2 = \frac{n^3}{3}-\frac{n^2}{2}+\frac{n}{6}

multiply-subtract pairs together with $n(n-1)/2$ divisions, so $\tfrac23 n^3 + O(n^2)$ individual floating-point operations. Each subsequent right-hand side costs $n(n-1)$ further multiply-subtract pairs and $n$ divisions, that is $O(n^2)$ : solving $k$ systems with the same matrix costs $n^3/3 + k n^2$ , not $k n^3/3$ . Kreyszig (Kreyszig (1978), §5.2, p. 307) quotes the figure as "about $n^3/3$ arithmetical operations", counting a multiply and the subtraction that follows it as one; the two counts agree up to that convention.

(iv) Determinant. $\det A = \det(P)^{-1}\prod_{i} u_{ii} = (-1)^{s}\prod_i u_{ii}$ , where $s$ is the number of row interchanges. This is the only practical way to compute a determinant of even moderate size; the cofactor expansion of Determinants costs $n!$ operations.

(v) Growth. With partial pivoting, every intermediate entry satisfies

\max_{i,j}\bigl|u_{ij}\bigr| \ \le\ 2^{\,n-1}\max_{i,j}\bigl|a_{ij}\bigr| ,

and the bound is attained. The quantity on the left over the quantity on the right is the growth factor, and it is what controls how much accuracy the elimination can lose.

Idea of the proofOriginal

Elimination is a factorization in disguise. Subtracting $\ell_{ik}$ times row $k$ from row $i$ is multiplication on the left by a matrix that differs from the identity in one entry; the product of all of them is unit lower triangular, and so is its inverse, whose entries are precisely the multipliers with no further arithmetic. So the triangular array of multipliers that elimination throws away is $L$ , and the final upper triangular array is $U$ . That is why one factorization serves any number of right-hand sides.

Pivoting enters because the algorithm divides by the pivot. If the pivot is small relative to the entries below it, the multipliers are large, the updated entries are large, and the original data are swamped: the computation loses the information it was given. Choosing the largest available pivot bounds every multiplier by $1$ , and the growth statement is the worst that can then happen, each elimination step at most doubling the largest entry.

The cost is one sum: at step $k$ a square block of side $n-k$ is updated, one multiply-subtract per entry.

ProofOriginal

(i) and (ii). Run the elimination. Suppose that after $k-1$ steps the matrix has the form $A^{(k)}$ , zero below the diagonal in the first $k-1$ columns. If every entry of column $k$ on or below the diagonal is zero, then the first $k$ columns of $A^{(k)}$ are linearly dependent; since $A^{(k)} = G A$ for an invertible $G$ (a product of the elementary operations so far), that would make $A$ singular. So for nonsingular $A$ some entry in that part of column $k$ is non-zero, and partial pivoting brings the largest in modulus, say $a^{(k)}_{pk}$ , onto the diagonal by exchanging rows $k$ and $p$ . The multipliers are then

\ell_{ik} = \frac{a^{(k)}_{ik}}{a^{(k)}_{kk}}, \qquad i > k, \qquad\text{so}\qquad |\ell_{ik}| \le 1

by the choice of pivot. Writing $E_k$ for the unit lower triangular matrix that performs step $k$ 's eliminations and $P_k$ for its row exchange, $n-1$ steps give $E_{n-1}P_{n-1}\cdots E_1 P_1 A = U$ . Pushing all the permutations to the right — each $P_j$ past an $E_i$ with $i<j$ only permutes the two rows of the multiplier column of $E_i$ , producing another unit lower triangular matrix with the same entries reordered — turns this into $\tilde E\, P A = U$ with $P = P_{n-1}\cdots P_1$ and $\tilde E$ unit lower triangular. Then $L = \tilde E^{-1}$ is unit lower triangular, and $PA = LU$ . If no interchange is ever needed, $P = I$ and (i) holds; a non-zero pivot at step $k$ is available without interchange exactly when $\det A_{1:k,1:k} \ne 0$ , since that determinant equals $u_{11}\cdots u_{kk}$ by the same factorization restricted to the leading block.

Uniqueness. If $L_1U_1 = L_2U_2$ with both $L_i$ unit lower triangular and both $U_i$ nonsingular upper triangular, then $L_2^{-1}L_1 = U_2U_1^{-1}$ . The left side is unit lower triangular and the right side is upper triangular, and a matrix that is both is the identity. Hence $L_1 = L_2$ and $U_1 = U_2$ .

(iii) Cost. At step $k$ the multipliers for rows $k+1,\dots,n$ cost $n-k$ divisions. Each of those $n-k$ rows is then updated in columns $k+1,\dots,n$ , which is $n-k$ entries, each one multiply and one subtract: $(n-k)^2$ pairs. Summing over $k = 1,\dots,n-1$ and substituting $j = n-k$ ,

\sum_{j=1}^{n-1}j^2 = \frac{(n-1)n(2n-1)}{6} = \frac{n^3}{3}-\frac{n^2}{2}+\frac{n}{6}, \qquad \sum_{j=1}^{n-1}j = \frac{n(n-1)}{2}.

Forward substitution computes $y_i = b_i - \sum_{j<i}\ell_{ij}y_j$ , costing $i-1$ pairs, so $n(n-1)/2$ in all and no divisions (the diagonal of $L$ is $1$ ); back substitution computes $x_i = (y_i - \sum_{j>i}u_{ij}x_j)/u_{ii}$ , costing $n-i$ pairs and one division, so $n(n-1)/2$ pairs and $n$ divisions.

(iv) Taking determinants in $PA = LU$ and using $\det L = 1$ (unit triangular) and $\det P = (-1)^{s}$ gives the formula.

(v) Growth. At step $k$ the updated entries are $a^{(k+1)}_{ij} = a^{(k)}_{ij} - \ell_{ik}a^{(k)}_{kj}$ . Since $|\ell_{ik}| \le 1$ ,

\bigl|a^{(k+1)}_{ij}\bigr| \le \bigl|a^{(k)}_{ij}\bigr| + \bigl|a^{(k)}_{kj}\bigr| \le 2\max_{i,j}\bigl|a^{(k)}_{ij}\bigr| .

So the largest modulus at most doubles at each of the $n-1$ steps, giving the factor $2^{n-1}$ . The bound is attained by the matrix with $1$ on the diagonal, $-1$ strictly below it and $1$ throughout the last column: computing with it for $n = 3,4,5,6$ gives growth factors exactly $4, 8, 16, 32$ , that is $2^{n-1}$ . $\square$

Method(Stationary iterations: when they converge)AdaptedAdapted from Stakgold & Holst (2011), §10.6, Algorithm 10.6.1, the error recursion (10.6.4) and Theorem 10.6.1, p. 773, and (10.6.27)-(10.6.30) and Theorem 10.6.4, pp. 799-800; Kreyszig (1978), §5.2, Theorem 5.2-1 with (5) and (7), p. 309, and 5.2-2, 5.2-3, pp. 310-312

Let $A = M - N$ be a splitting and $C = M^{-1}N = I - M^{-1}A$ the iteration matrix, so that the error $e^m = x - x^m$ obeys $e^{m+1} = C e^{m}$ and hence $e^m = C^m e^0$ .

(i) The exact criterion. The iteration converges to $x$ for every starting vector $x^0$ if and only if

\rho(C) < 1 ,

$\rho$ being the spectral radius. This is Theorem 10.6.1 of Stakgold and Holst (Stakgold & Holst (2011), §10.6, p. 773), and the asymptotic error reduction per step is $\rho(C)$ : the number of steps to gain one decimal digit is $-1/\log_{10}\rho(C)$ .

(ii) A checkable sufficient condition, with error bounds. If

\alpha := \|C\|_\infty = \max_{j}\sum_{k=1}^{n}|c_{jk}| < 1

then the iteration converges, and Kreyszig's Theorem 5.2-1 (Kreyszig (1978), §5.2, Thm 5.2-1 and Eq. (7), p. 309) gives the bounds, in the max metric $d(x,y) = \max_j|x_j - y_j|$ ,

d(x^{(m)}, x) \ \le\ \frac{\alpha}{1-\alpha}\,d\bigl(x^{(m-1)},x^{(m)}\bigr) \ \le\ \frac{\alpha^{m}}{1-\alpha}\,d\bigl(x^{(0)},x^{(1)}\bigr) .

The first is computable from two successive iterates and is the one to use as a stopping test; the second is computable in advance.

(iii) Jacobi. $C_J = D^{-1}(L+U) = I - D^{-1}A$ , that is $x^{m+1}_j = \bigl(b_j - \sum_{k\ne j}a_{jk}x^m_k\bigr)/a_{jj}$ , which needs $a_{jj} \ne 0$ . Condition (ii) becomes strict row diagonal dominance,

\sum_{k \ne j}|a_{jk}| < |a_{jj}| \qquad (j = 1,\dots,n) ,

which is Kreyszig's row sum criterion (12*).

(iv) Gauss-Seidel. $C_{GS} = (D-L)^{-1}U$ , that is $x^{m+1}_j = \bigl(b_j - \sum_{k<j}a_{jk}x^{m+1}_k - \sum_{k>j}a_{jk}x^{m}_k\bigr)/a_{jj}$ : the components already updated in the current sweep are used immediately. Kreyszig states, without proof, that strict row diagonal dominance is sufficient for this iteration too.

(v) The two are not comparable. Kreyszig's Problems 5 and 6 (p. 313) exhibit a $3\times3$ system on which Jacobi diverges and Gauss-Seidel converges, and another on which Jacobi converges and Gauss-Seidel diverges. Neither method dominates the other in general. What is true in the important special case is recorded by Stakgold and Holst (Stakgold & Holst (2011), §10.6, p. 800): for consistently ordered matrices with Property A — which the standard finite-difference Laplacian is — $\rho(C_{GS}) = \rho(C_J)^2$ , so Gauss-Seidel does in one sweep what Jacobi does in two.

(vi) The bad news for discretized PDEs. For the box-method discretization of Poisson's equation on a uniform mesh of spacing $h$ , Stakgold and Holst give $\rho(C_J) = \cos \pi h$ and $\rho(C_{GS}) = \cos^2\pi h$ , both of the form $1 - O(h^2)$ , while the optimally relaxed SOR iteration achieves $1 - O(h)$ . As $h \to 0$ all of them approach $1$ : the classical iterations slow down without limit as the mesh is refined. This is the reason multigrid and preconditioned Krylov methods exist (Relaxation methods for elliptic problems).

Idea of the proofOriginal

Everything follows from $e^m = C^m e^0$ . A matrix power tends to zero exactly when every eigenvalue has modulus below one: the Jordan form turns $C^m$ into blocks $\lambda^m$ times a polynomial in $m$ , and an exponential beats a polynomial in either direction. The sufficient condition is cruder and cheaper: an induced norm dominates the spectral radius, so a norm below one forces the spectral radius below one, and it does so with a constant one can compute from the matrix entries without finding a single eigenvalue. Kreyszig gets it as a contraction on $\mathbb{R}^n$ with the max metric, and the error bounds are then nothing but the standard bounds of the Banach fixed point theorem.

The difference between Jacobi and Gauss-Seidel is that Gauss-Seidel overwrites in place. It costs the same arithmetic and less storage, and on the matrices that come from discretizing an elliptic equation it halves the number of sweeps — but that is a property of those matrices, not a theorem about the two methods.

ProofOriginal

(i) Subtracting $x = Cx + d$ from $x^{m+1} = Cx^m + d$ gives $e^{m+1} = Ce^m$ , hence $e^m = C^m e^0$ , so the question is whether $C^m \to 0$ .

If $\rho(C) \ge 1$ : take an eigenvalue $\lambda$ with $|\lambda| \ge 1$ and an eigenvector $v \ne 0$ , and start at $x^0 = x - v$ , so $e^0 = v$ . Then $e^m = \lambda^m v$ does not tend to $0$ .

If $\rho(C) < 1$ : put $C = S J S^{-1}$ in Jordan form (Jordan normal form). It suffices that $J^m \to 0$ , and $J$ is block diagonal, so it suffices that $J_k^m \to 0$ for one Jordan block $J_k = \lambda I + N$ of size $s$ , where $N$ is nilpotent with $N^{s}=0$ . Since $\lambda I$ and $N$ commute, the binomial theorem gives

J_k^{m} = \sum_{t=0}^{s-1}\binom{m}{t}\lambda^{m-t}N^{t},

a sum of $s$ terms each bounded by $m^{t}|\lambda|^{m-t}\|N\|^t/t!$ . With $|\lambda| < 1$ each $m^{t}|\lambda|^{m}$ tends to $0$ , because $m^t |\lambda|^m = \exp(t\log m - m\log(1/|\lambda|)) \to 0$ . Hence $J_k^m \to 0$ .

(ii) The map $T x = Cx + d$ satisfies, in the max norm,

\|Tx - Ty\|_\infty = \|C(x-y)\|_\infty = \max_j\Bigl|\sum_k c_{jk}(x_k-y_k)\Bigr| \le \Bigl(\max_j \sum_k |c_{jk}|\Bigr)\|x-y\|_\infty = \alpha\,\|x-y\|_\infty ,

so $T$ is a contraction with constant $\alpha$ on the complete metric space $(\mathbb{R}^n, d)$ . The Banach fixed point theorem (The contraction mapping theorem) gives a unique fixed point, convergence from any start, and exactly the two error bounds quoted. (Consistency with (i): $\rho(C) \le \|C\|_\infty = \alpha < 1$ , since an induced norm bounds the spectral radius.)

(iii) $M = D$ is invertible precisely when no diagonal entry vanishes, and then $C_J = I - D^{-1}A$ has entries $c_{jk} = -a_{jk}/a_{jj}$ for $k\ne j$ and $0$ for $k = j$ . So $\alpha = \max_j \sum_{k\ne j}|a_{jk}|/|a_{jj}|$ , and $\alpha < 1$ is the stated dominance.

(iv) $M = D - L$ is lower triangular with the $a_{jj}$ on the diagonal, hence invertible under the same condition, and solving $Mx^{m+1} = Ux^m + b$ one row at a time from the top gives the component formula, since row $j$ involves $x^{m+1}_k$ only for $k < j$ . The claim about diagonal dominance is Kreyszig's and is reported, not proved, here: he states at the end of §5.2 that condition (12) is sufficient for the Gauss-Seidel iteration as well, referring the reader elsewhere for better conditions.

(v) and (vi) are quoted from the two sources as stated; no proof is offered here. $\square$

Method(The conjugate gradient method)AdaptedAdapted from Stakgold & Holst (2011), §10.6, Algorithm 10.6.2 with B = I, p. 778, the Krylov minimisation and minimax reduction, p. 779, and the contraction bound (10.6.9), p. 780

Let $A$ be symmetric positive definite, $b$ given, $x$ the exact solution. Write $(u,v)_A = u^{\mathsf T}Av$ and $\|u\|_A = (u,u)_A^{1/2}$ for the energy inner product and norm, and

K_i = \operatorname{span}\bigl\{r^0,\, Ar^0,\, \dots,\, A^{i-1}r^0\bigr\}

for the $i$ -th Krylov subspace. The method of Hestenes and Stiefel is Algorithm 10.6.2 of Stakgold and Holst (Stakgold & Holst (2011), §10.6, p. 778) with the preconditioner taken to be the identity: starting from any $x^0$ , with $r^0 = b - Ax^0$ and $p^0 = r^0$ ,

\alpha_i = \frac{(r^i, r^i)}{(p^i, Ap^i)},\quad x^{i+1} = x^i + \alpha_i p^i,\quad r^{i+1} = r^i - \alpha_i Ap^i,\quad \beta_{i+1} = \frac{(r^{i+1},r^{i+1})}{(r^i,r^i)},\quad p^{i+1} = r^{i+1}+\beta_{i+1}p^i .

Each step costs one matrix-vector product, two inner products and three vector updates; $A$ is never modified and need not even be stored, only applied.

Then, in exact arithmetic and as long as $r^i \ne 0$ :

Orthogonality. $(r^i, r^j) = 0$ and $(p^i, Ap^j) = 0$ for $j < i$ , and $K_i = \operatorname{span}\{r^0,\dots,r^{i-1}\} = \operatorname{span}\{p^0,\dots,p^{i-1}\}$ .
Optimality. $x^i$ is the unique minimiser of $\|x - u\|_A$ over the affine subspace $u \in x^0 + K_i$ . No method that uses only the vectors $A^{j}r^0$ can do better in this norm.
Finite termination. $x^{i} = x$ for some $i \le n$ ; more sharply, for some $i \le$ the number of distinct eigenvalues of $A$ .
Rate. With $\kappa = \kappa_2(A) = \lambda_{\max}/\lambda_{\min}$ ,

\|x - x^{i}\|_A \ \le\ 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{\!i}\,\|x-x^0\|_A ,

which is Stakgold and Holst's bound (10.6.9) with $B = I$ . Compare a stationary iteration, whose rate involves $\kappa$ itself rather than $\sqrt{\kappa}$ : the number of conjugate gradient steps needed for a fixed accuracy grows like $\sqrt{\kappa}$ , not like $\kappa$ .

Positive definiteness is used throughout and cannot be dropped: $(p,Ap)$ appears in a denominator, and $\|\cdot\|_A$ is a norm only for positive definite $A$ .

Idea of the proofOriginal

For symmetric positive definite $A$ , solving $Ax=b$ is the same as minimising the quadratic $\phi(u) = \tfrac12 u^{\mathsf T}Au - b^{\mathsf T}u$ , since $\phi(u) = \tfrac12\|x-u\|_A^2 - \tfrac12 x^{\mathsf T}Ax$ . Steepest descent takes the negative gradient, which is the residual, and immediately wastes its work: the next step undoes part of the previous one. The conjugate gradient fix is to search instead along directions that are orthogonal in the energy inner product, so that minimising along a new direction never spoils the minimisation already achieved along the old ones. The iterate after $i$ steps is then the exact minimiser over the whole $i$ -dimensional space searched so far, not merely over the last line.

The miracle — and it is specific to symmetric $A$ — is that building such a direction needs only the previous one. The space searched after $i$ steps is the Krylov space, and the single number $\beta$ suffices to keep the new direction energy-orthogonal to all its predecessors, because a symmetric matrix applied to a vector of one Krylov space lands in the next, and everything earlier has already been orthogonalised away.

Optimality over Krylov spaces then converts the convergence question into a question about polynomials: which polynomial of degree $i$ equal to $1$ at the origin is smallest on the spectrum of $A$ ? If the spectrum is clustered the answer is tiny, which is why preconditioning works; if the spectrum fills an interval $[\lambda_{\min},\lambda_{\max}]$ , the answer is given by Chebyshev polynomials and produces the $\sqrt{\kappa}$ rate.

ProofOriginal

1. Orthogonality. We prove by induction on $i$ the statement $H_i$ : the vectors $r^0,\dots,r^i$ are mutually orthogonal, $p^0,\dots,p^i$ are mutually $A$ -orthogonal, and

\operatorname{span}\{r^0,\dots,r^i\} = \operatorname{span}\{p^0,\dots,p^i\} = K_{i+1}.

$H_0$ holds because $p^0 = r^0$ and $K_1 = \operatorname{span}\{r^0\}$ . Assume $H_i$ and $r^{i+1}\ne 0$ .

Step 1: $r^{i+1} \perp p^j$ for $j \le i$ . For $j = i$ ,

(r^{i+1},p^i) = (r^i,p^i) - \alpha_i(Ap^i,p^i) .

Now $p^i = r^i + \beta_i p^{i-1}$ and $r^i \perp p^{i-1}$ (that is $H_i$ : $p^{i-1} \in K_i = \operatorname{span}\{r^0,\dots,r^{i-1}\}$ and $r^i$ is orthogonal to each of those), so $(r^i,p^i) = (r^i,r^i)$ , and the definition of $\alpha_i$ makes the right side vanish. For $j < i$ : $(r^i,p^j) = 0$ because $p^j \in K_i$ , and $(Ap^i,p^j) = 0$ by the $A$ -orthogonality in $H_i$ . Hence $(r^{i+1},p^j)=0$ .

Since $\operatorname{span}\{p^0,\dots,p^i\} = \operatorname{span}\{r^0,\dots,r^i\}$ by $H_i$ , step 1 also gives $r^{i+1}\perp r^j$ for $j \le i$ .

Step 2: $p^{i+1}$ is $A$ -orthogonal to $p^0,\dots,p^i$ . From $r^{j+1} = r^j - \alpha_j Ap^j$ we get $Ap^j = (r^j - r^{j+1})/\alpha_j$ whenever $\alpha_j \ne 0$ , which holds because $(p^j,Ap^j)>0$ and $(r^j,r^j)>0$ . For $j < i$ ,

(Ap^{i+1},p^j) = (Ar^{i+1},p^j) + \beta_{i+1}(Ap^i,p^j) = \bigl(r^{i+1}, Ap^j\bigr) + 0 = \frac{1}{\alpha_j}\bigl(r^{i+1}, r^j - r^{j+1}\bigr) = 0 ,

using the symmetry of $A$ and step 1, since $j+1 \le i$ . For $j = i$ ,

(Ap^{i+1},p^i) = (r^{i+1},Ap^i) + \beta_{i+1}(p^i,Ap^i) = \frac{1}{\alpha_i}\bigl(r^{i+1}, r^i - r^{i+1}\bigr) + \beta_{i+1}\frac{(r^i,r^i)}{\alpha_i} = \frac{1}{\alpha_i}\Bigl(-(r^{i+1},r^{i+1}) + \beta_{i+1}(r^i,r^i)\Bigr) = 0

by the definition of $\beta_{i+1}$ . This is where the formula for $\beta$ comes from: it is the unique choice that makes the new direction $A$ -orthogonal to the last one, and steps 1 and 2 show that this one condition secures orthogonality to all the others for free.

Step 3: the spans. $r^{i+1} = r^i - \alpha_i Ap^i$ with $r^i \in K_{i+1}$ and $p^i \in K_{i+1}$ , so $Ap^i \in K_{i+2}$ and $r^{i+1} \in K_{i+2}$ ; likewise $p^{i+1} = r^{i+1}+\beta_{i+1}p^i \in K_{i+2}$ . Conversely $\operatorname{span}\{r^0,\dots,r^{i+1}\}$ has dimension $i+2$ (the $r^j$ are non-zero and mutually orthogonal), is contained in $K_{i+2}$ , and $\dim K_{i+2}\le i+2$ ; so the two are equal, and the same argument applies to the $p^j$ , which are non-zero and mutually $A$ -orthogonal. That is $H_{i+1}$ .

2. Optimality. By construction $x^i \in x^0 + K_i$ . Let $e^i = x - x^i$ ; then $Ae^i = b - Ax^i = r^i$ , so for every $v \in K_i$ ,

(e^i, v)_A = (Ae^i, v) = (r^i, v) = 0

by part 1. So $e^i$ is $A$ -orthogonal to $K_i$ , that is, $x^i$ is the $A$ -orthogonal projection of $x$ onto the affine subspace $x^0+K_i$ . For any other $u = x^i + w$ with $w\in K_i$ ,

\|x-u\|_A^2 = \|e^i - w\|_A^2 = \|e^i\|_A^2 - 2(e^i,w)_A + \|w\|_A^2 = \|e^i\|_A^2 + \|w\|_A^2 > \|e^i\|_A^2

unless $w = 0$ (Orthogonal projection and least squares).

3. Finite termination. The non-zero residuals are mutually orthogonal in $\mathbb{R}^n$ , so there are at most $n$ of them and $r^i = 0$ , hence $x^i = x$ , for some $i \le n$ . For the sharper count, let $\lambda_1,\dots,\lambda_m$ be the distinct eigenvalues of $A$ and put $q(\lambda) = \prod_{j=1}^{m}(1-\lambda/\lambda_j)$ , a polynomial of degree $m$ with $q(0)=1$ . Since $x^i - x^0 \in K_i$ , we may write $x^i = x^0 + s(A)r^0$ with $\deg s \le i-1$ , and then

e^i = e^0 - s(A)Ae^0 = q_i(A)e^0, \qquad q_i(\lambda) = 1-\lambda s(\lambda), \quad q_i(0)=1, \ \deg q_i \le i ,

and conversely every such $q_i$ arises from some $u \in x^0+K_i$ . By part 2, $\|e^i\|_A = \min\{\|q(A)e^0\|_A : \deg q \le i,\ q(0)=1\}$ . Taking $i = m$ and the $q$ above, $q(A) = 0$ because $q$ vanishes at every eigenvalue and $A$ is diagonalisable (Spectral theorem for Hermitian matrices); hence $\|e^m\|_A = 0$ .

\frac{\|e^i\|_A}{\|e^0\|_A} \ \le\ \min_{\substack{\deg q\le i\\ q(0)=1}}\ \max_{\lambda\in\sigma(A)}\bigl|q(\lambda)\bigr| \ \le\ \min_{\substack{\deg q\le i\\ q(0)=1}}\ \max_{\lambda_{\min}\le\lambda\le\lambda_{\max}}\bigl|q(\lambda)\bigr| .

The last minimax problem is solved by the Chebyshev polynomial $T_i$ rescaled to $[\lambda_{\min},\lambda_{\max}]$ and normalised at the origin (Chebyshev polynomials); Stakgold and Holst (Stakgold & Holst (2011), §10.6, Eq. (10.6.9) and Exercises 10.5.2-10.5.3, pp. 779-780) carry out that step and record the resulting bound, which is the one quoted. $\square$

Remark(Where the hypotheses are used)Original

Pivoting is not a refinement; without it the algorithm can return nonsense on a perfectly conditioned problem. Take

\varepsilon x_1 + x_2 = 1, \qquad x_1 + x_2 = 2, \qquad \varepsilon = 10^{-17},

whose exact solution is $x_1 = 1/(1-\varepsilon)$ , $x_2 = (1-2\varepsilon)/(1-\varepsilon)$ , both equal to $1$ to sixteen digits. Its condition number is $\kappa_2 = 2.618$ : the problem is as well conditioned as a problem can be. Eliminating without pivoting uses the multiplier $1/\varepsilon = 10^{17}$ ; in IEEE double precision $1 - 10^{17}$ and $2 - 10^{17}$ both round to $-10^{17}$ , the $1$ and the $2$ being lost entirely, so $x_2 = 1$ and then $x_1 = (1 - x_2)/\varepsilon = 0$ . The computed answer is $(0,1)$ against the true $(1,1)$ : one hundred per cent error. Exchanging the two rows first gives multiplier $\varepsilon$ and returns $(1,1)$ to full precision. The information destroyed is the coefficient $1$ in the second equation, and what destroys it is the subtraction of a huge number from a small one — cancellation, in the sense of Floating-point arithmetic and error.

The growth factor is why "with partial pivoting" is enough in practice but not in theory. The bound $2^{n-1}$ of part (v) is attained, so there exist matrices on which elimination with pivoting loses $n-1$ bits. They are rare; on matrices arising from applications the growth factor is almost always below $n$ . The theorem is what it is, and a code that must be certain checks the residual afterwards.
$\rho(C)<1$ is necessary and sufficient; $\|C\|<1$ is only sufficient. The gap is real. Diagonal dominance is a sufficient condition for Jacobi, and Kreyszig's Problems 5 and 6 show it is not necessary and that the two classical methods do not dominate one another.
$a_{jj}\neq 0$ for both iterations. A zero on the diagonal makes $M$ singular and the iteration undefined, even though the system may be perfectly solvable; a row exchange fixes it, so the hypothesis is about the arrangement of the equations, not about the problem.
Symmetry and positive definiteness for the conjugate gradient method. Symmetry is used in step 2 of the orthogonality proof, where $(Ar,p) = (r,Ap)$ , and it is what makes the short three-term recurrence possible. Positive definiteness is used for $\|\cdot\|_A$ to be a norm, for $\alpha_i$ to have a non-zero denominator, and in part 4 where $\lambda_k > 0$ . Applied to a non-symmetric or indefinite matrix the recurrence can break down with a zero denominator; the family of methods for that case (GMRES, BiCGSTAB) must store or regenerate the whole Krylov basis.
Exact arithmetic for finite termination. Part 3 is a statement about exact arithmetic only. In floating point the orthogonality of the residuals is lost gradually and the method does not terminate at step $n$ . This does not matter, because in practice it is stopped long before $n$ on the strength of part 4 — which is a statement about the rate, and survives rounding.
The residual is not the error. If $\hat x$ is any approximation with residual $r = b - A\hat x$ , then $x - \hat x = A^{-1}r$ and $\|b\| = \|Ax\|\le\|A\|\|x\|$ , so

\frac{\|x-\hat x\|}{\|x\|}\ \le\ \|A^{-1}\|\,\|r\|\cdot\frac{\|A\|}{\|b\|} = \kappa(A)\,\frac{\|r\|}{\|b\|} ,

and Matrix norms and condition numbers shows this is attained. A tiny residual certifies a tiny error only when $\kappa(A)$ is modest. Every stopping test on $\|r\|$ is therefore a statement about the error only after multiplication by an estimate of $\kappa$ .

Method(Choosing a solver)Original

Is the matrix sparse and large? If $n \lesssim 10^4$ and $A$ is dense, factorize: the cost $n^3/3$ is predictable and the answer is backward stable. If $A$ is sparse with $O(n)$ non-zeros and $n$ is large, a dense factorization is impossible on both counts — fill-in destroys the sparsity and the cost is cubic — so iterate.
Factorize once, solve many times. Never form $A^{-1}$ . Compute $PA = LU$ once at cost $n^3/3$ and reuse it for every right-hand side at $n^2$ each. Forming $A^{-1}$ costs the same as $n$ solves, needs $n^2$ storage, and then each application still costs $n^2$ — with a worse error.
Exploit structure before anything else. A tridiagonal system costs $O(n)$ , not $O(n^3)$ (Polynomial and spline interpolation uses this for splines, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §27.3.3, p. 998) give the recurrence); a banded system of bandwidth $m$ costs $O(nm^2)$ ; a symmetric positive definite system uses the Cholesky factorization $A = R^{\mathsf T}R$ at half the cost of LU and needs no pivoting.
For an iterative solve, check the splitting first. Verify strict diagonal dominance; if it holds, Jacobi and Gauss-Seidel both converge and the row-sum constant $\alpha$ gives a usable a priori step count. If it fails, compute or estimate $\rho(C)$ before spending time: an iteration with $\rho$ just below $1$ is worse than useless.
For symmetric positive definite $A$ , use conjugate gradients and precondition it. The step count scales as $\sqrt{\kappa}$ , and preconditioning — solving $M^{-1}Ax = M^{-1}b$ with $M$ an easily invertible approximation to $A$ — is the whole game, because it is the condition number of the preconditioned matrix that enters. Incomplete Cholesky and multigrid (Relaxation methods for elliptic problems) are the standard choices.
Stop on the right quantity. Use $\|r\|/\|b\| \le \text{tol}$ only with an estimate of $\kappa(A)$ in hand; otherwise use the a posteriori bound of the stationary-iteration statement, which involves only successive iterates.
Always check afterwards. Compute the residual of the returned $\hat x$ in higher precision if you can. One step of iterative refinement — solve $A\delta = r$ with the factorization you already have, set $\hat x \leftarrow \hat x + \delta$ — costs $n^2$ and typically recovers the digits lost to a large growth factor.

Example(The same three ideas on explicit numbers)Original

Part 1: elimination with interchange. Solve $Ax=b$ with

A = \begin{pmatrix} 1 & 4 & 2\\ 3 & 1 & -1\\ 2 & -5 & 5\end{pmatrix}, \qquad b = \begin{pmatrix}-1\\-2\\27\end{pmatrix}.

Step 1. The largest entry of column $1$ in modulus is $3$ , in row $2$ : exchange rows $1$ and $2$ . The multipliers are $\ell = 1/3$ for the (new) row $2$ and $\ell = 2/3$ for row $3$ , both at most $1$ in modulus as the theory promises. Subtracting:

\begin{pmatrix} 3 & 1 & -1 & \big| & -2\\ 0 & \tfrac{11}{3} & \tfrac{7}{3} & \big| & -\tfrac13\\ 0 & -\tfrac{17}{3} & \tfrac{17}{3} & \big| & \tfrac{85}{3}\end{pmatrix}.

Step 2. $|{-17/3}| > |11/3|$ , so exchange rows $2$ and $3$ — a second interchange, which plain elimination would not have made. The multiplier is $(11/3)/(-17/3) = -11/17$ , and subtracting it gives the third row $6x_3 = 18$ .

Back substitution. $x_3 = 3$ ; then $-\tfrac{17}{3}x_2 + \tfrac{17}{3}(3) = \tfrac{85}{3}$ gives $x_2 = -2$ ; then $3x_1 + (-2) - 3 = -2$ gives $x_1 = 1$ . So

x = (1,\,-2,\,3)^{\mathsf T},

which substituted back into the three original equations gives $1-8+6=-1$ , $3-2-3=-2$ and $2+10+15=27$ : the data exactly.

The factorization. The permutation sends the original rows to the order $(2,3,1)$ , and

L = \begin{pmatrix}1&0&0\\ \tfrac23&1&0\\ \tfrac13&-\tfrac{11}{17}&1\end{pmatrix}, \qquad U = \begin{pmatrix}3&1&-1\\0&-\tfrac{17}{3}&\tfrac{17}{3}\\0&0&6\end{pmatrix}.

Multiplying out, row $3$ of $LU$ is $\tfrac13(3,1,-1) - \tfrac{11}{17}(0,-\tfrac{17}{3},\tfrac{17}{3}) + (0,0,6) = (1,4,2)$ , the original first row, and the other two rows check the same way: $PA = LU$ exactly. The determinant is $\det A = (+1)\cdot 3\cdot(-\tfrac{17}{3})\cdot 6 = -102$ , the sign being $+$ because two interchanges were made; expanding the $3\times3$ determinant by cofactors gives $0 - 68 - 34 = -102$ as well. The operation count predicted by part (iii) is $n^3/3 - n^2/2 + n/6 = 9 - 4.5 + 0.5 = 5$ multiply-subtract pairs and $n(n-1)/2 = 3$ divisions, which is exactly what the two steps above used: a $2\times2$ block updated at step $1$ and a $1\times1$ at step $2$ , with three multipliers.

Part 2: Jacobi and Gauss-Seidel. Take the symmetric positive definite

A = \begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}, \qquad b = \begin{pmatrix}2\\4\\10\end{pmatrix}, \qquad x = \begin{pmatrix}1\\2\\3\end{pmatrix},

the solution being exact as substitution confirms ( $4-2=2$ , $-1+8-3=4$ , $-2+12=10$ ). The rows of $|C_J|$ sum to $0.25$ , $0.5$ , $0.25$ , so $\alpha = 0.5 < 1$ : strictly diagonally dominant, and both iterations converge. The eigenvalues of $A$ are $4-\sqrt2, 4, 4+\sqrt2$ , so

\rho(C_J) = \frac{2\cos(\pi/4)}{4} = \frac{\sqrt2}{4} = 0.3535534, \qquad \rho(C_{GS}) = 0.125 = \rho(C_J)^2 ,

the second computed from the eigenvalues of $(D-L)^{-1}U$ and agreeing exactly with the squaring relation Stakgold and Holst record for this class of matrix. Starting from $x^0 = 0$ the maximum-norm errors are

sweep $m$	Jacobi	Gauss-Seidel
1	$1.0000$	$8.7500\times10^{-1}$
2	$2.5000\times10^{-1}$	$2.1875\times10^{-1}$
3	$1.2500\times10^{-1}$	$2.7344\times10^{-2}$
4	$3.1250\times10^{-2}$	$3.4180\times10^{-3}$
6	$3.9062\times10^{-3}$	$5.3406\times10^{-5}$
8	$4.8828\times10^{-4}$	$8.3447\times10^{-7}$
12	$7.6294\times10^{-6}$

The Gauss-Seidel errors fall by exactly $1/8$ per sweep from $m=2$ on, which is $\rho(C_{GS})$ . The Jacobi errors alternate between factors $1/4$ and $1/2$ — the max-norm picks out different components on alternate sweeps — but over two sweeps they fall by $1/8 = \rho(C_J)^2$ , as the spectral radius requires. Gauss-Seidel reaches $10^{-6}$ in $8$ sweeps ( $8.34\times10^{-7}$ ) where Jacobi needs $14$ ( $9.54\times10^{-7}$ at sweep $14$ , against $3.81\times10^{-6}$ at sweep $13$ ). Kreyszig's a priori bound for Jacobi, with $\alpha = 0.5$ and $d(x^{(0)},x^{(1)}) = \|D^{-1}b\|_\infty = 2.5$ , predicts an error at most $\alpha^{m}(1-\alpha)^{-1}(2.5) = 5\cdot 2^{-m}$ , which at $m=12$ is $1.22\times10^{-3}$ : a true bound, pessimistic here by a factor of $160$ , as a bound built from a norm rather than a spectral radius must be.

Part 3: conjugate gradients on the same system. $\kappa_2(A) = (4+\sqrt2)/(4-\sqrt2) = 2.093836$ , so the predicted contraction factor is $(\sqrt\kappa-1)/(\sqrt\kappa+1) = 0.182676$ . From $x^0 = 0$ , so $r^0 = p^0 = b$ and $\|x - x^0\|_A = \sqrt{x^{\mathsf T}\!Ax} = \sqrt{x^{\mathsf T}b} = \sqrt{40} = 6.3246$ :

$i$	$\alpha_{i-1}$	$x^{i}$	$\\|r^i\\|_2$	$\\|x-x^i\\|_A$
1	$0.3125000$	$(0.6250000,\ 1.2500000,\ 3.1250000)$	$3.1125$	$1.5811$
2	$0.2484140$	$(0.8514190,\ 2.0133556,\ 3.0150250)$	$0.63751$	$0.30576$
3	$0.2300307$	$(1.0000000,\ 2.0000000,\ 3.0000000)$	$4.4\times10^{-17}$	$0$

Three steps for a $3\times3$ matrix with three distinct eigenvalues, exactly as part 3 of the statement requires; the residual at step $3$ is at the rounding level and the energy error is zero in our run. Check the bound of part 4 as well. At $i=2$ it reads

\|x-x^2\|_A \;\le\; 2\rho^{2}\,\|x-x^0\|_A = 2(0.182676)^2(6.3246) = 0.42211 ,

against the achieved $0.30576$ : the bound holds, with about $28\%$ to spare. In ratio form the achieved contraction is $0.30576/6.3246 = 0.048345$ , comfortably under the allowed $2\rho^2 = 0.066741$ . At $i=1$ it reads $2\rho\,(6.3246) = 0.36535 \times 6.3246 = 2.3107$ against the achieved $1.5811$ : again true, again loose. The reason the bound has room is that it is a minimax value over the whole interval $[\lambda_{\min},\lambda_{\max}]$ , whereas conjugate gradients only has to be small on the three points $\{4-\sqrt2,\,4,\,4+\sqrt2\}$ actually in the spectrum; a finite spectrum is easier than a full interval, which is both why the bound has slack here and why clustering the spectrum by preconditioning pays.

Application(Where the cubic cost and the condition number decide the science)Original

Structural and finite element analysis. Assembling a finite element model (The finite element method) produces a symmetric positive definite stiffness matrix that is sparse and very large: a million degrees of freedom is routine, a billion is not unusual. A dense factorization would need $10^{18}$ operations and $10^{12}$ words of storage for the million-unknown case; a preconditioned conjugate gradient solve needs a few hundred matrix-vector products, each $O(n)$ . The condition number of a finite element stiffness matrix grows like $h^{-2}$ as the mesh is refined, so $\sqrt{\kappa} \sim h^{-1}$ , and the unpreconditioned step count doubles when the mesh is halved — which is exactly why industrial codes are built around multigrid preconditioners that make $\kappa$ independent of $h$ .

Least squares and data fitting. Fitting $m$ data with $n$ parameters leads to the normal equations $A^{\mathsf T}A c = A^{\mathsf T}y$ (Orthogonal projection and least squares), which are symmetric positive definite — but $\kappa(A^{\mathsf T}A) = \kappa(A)^2$ . Fitting a degree- $12$ polynomial in the monomial basis on $[0,1]$ uses $13$ coefficients, so the normal matrix is the $13\times13$ Hilbert matrix, whose spectral condition number we compute to be $5.6\times10^{17}$ — the $12\times12$ one, belonging to the degree- $11$ fit, is already $1.7\times10^{16}$ — and no algorithm whatever can return a meaningful answer in double precision: by the rule of thumb in Matrix norms and condition numbers, every digit is gone. The cure is not a better solver but a better basis — orthogonal polynomials (General properties of orthogonal polynomials) — or a QR factorization of $A$ that never forms $A^{\mathsf T}A$ and so works with $\kappa(A)$ rather than its square.

Implicit time stepping. Every implicit integrator for a stiff system (Stability and stiff equations) and every implicit scheme for a parabolic PDE requires the solution of a linear system at every time step, with the same matrix if the step size is held fixed. That is the canonical case for factorizing once and reusing: the $n^3/3$ is paid at the first step and $n^2$ at every step thereafter. When the step size is adapted, the trade-off between refactorizing and iterating is the central implementation decision in the solver.

Circuit simulation and network flow. Modified nodal analysis produces a sparse, highly unsymmetric matrix whose structure follows the circuit topology, and here pivoting is not an academic concern: the entries span many orders of magnitude, from gate capacitances to supply resistances, and a pivot chosen by position rather than by magnitude will produce the $(0,1)$ of the remark above in a circuit with a billion transistors. Production simulators use a pivoting strategy that balances the magnitude criterion against keeping the factorization sparse, which is the practical form of the growth-factor bound.

Method

Computing eigenvalues

Adapted

Add to PDF

The three iterations every practical eigenvalue computation is built from - power iteration, whose error falls by the ratio of the two largest eigenvalues at each step; inverse iteration with a shift, which turns that ratio into a quantity the user controls and becomes cubically convergent when the shift is the Rayleigh quotient; and the QR algorithm, which is power iteration run on a whole nested family of subspaces at once - with their convergence rates, their costs and a residual bracket that says when to stop.

Prerequisites: Spectral theorem for Hermitian matrices, Eigenvalues and eigenvectors, Orthonormal bases and Gram-Schmidt, Quadratic forms and the Rayleigh quotient, Matrix norms and condition numbers, Solving linear systems numerically, Floating-point arithmetic and error

Note(What it means)AdaptedAdapted from Hassani (2013), §6.2, Eqs. (6.5)-(6.7) and Thm 6.2.5, pp. 173-174; Arfken, Weber & Harris (2013), §6.4, Finding a Diagonalizing Transformation, p. 314

The eigenvalues of an $n\times n$ matrix are the roots of its characteristic polynomial, a polynomial of degree $n$ (Hassani (2013), Eq. (6.7), p. 173). From the fifth degree onwards there is no formula in radicals for the roots of a general polynomial, so there is no finite sequence of arithmetic operations and root extractions that produces the eigenvalues of a general matrix. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §6.4, p. 314) draw the practical conclusion: numerical diagonalization proceeds by successive approximation, and it has to.

Every method in use is therefore an iteration, and essentially all of them are built from one idea: multiply a vector by a matrix repeatedly and the direction of largest growth takes over. The three sections below are that idea, applied to $A$ (power iteration), to $(A-\sigma I)^{-1}$ (inverse iteration, which lets you aim at any part of the spectrum you like), and to a whole nested family of subspaces at once (the QR algorithm, which is how a dense eigenvalue problem is actually solved). Forming the characteristic polynomial and sending it to a root finder (Root finding) is not on the list: its coefficients are catastrophically ill-conditioned functions of the matrix entries, and the information they carry about the small eigenvalues is destroyed before the root finder is called.

Definition(The iterations, the Rayleigh quotient and the residual)AdaptedAdapted from Kreyszig (1978), §9.2, Eq. (1) and Thm 9.2-1, p. 465, and Problems 8-9, p. 469; Stakgold & Holst (2011), §6.5, (6.5.17) and (6.5.19)-(6.5.20), pp. 401-402

Throughout, $A \in M_n(\mathbb{C})$ , vectors are columns, $(x,y) = y^{*}x$ (the handbook's convention, linear in its first argument), and $\lVert x\rVert = (x,x)^{1/2}$ .

The Rayleigh quotient of a non-zero vector $q$ is

\rho(q) \;=\; \frac{(Aq,q)}{(q,q)} ,

the name used by Kreyszig (Kreyszig (1978), §9.2, Problem 8, p. 469), who records that for a bounded self-adjoint operator the whole spectrum lies between the infimum and the supremum of $\rho$ ; for a Hermitian matrix with eigenvalues $\lambda_1 \ge \cdots \ge \lambda_n$ the extremes are attained, $\lambda_1 = \max_q \rho(q)$ and $\lambda_n = \min_q \rho(q)$ (Quadratic forms and the Rayleigh quotient). The residual of a unit vector $q$ is

\eta(q) \;=\; \lVert Aq - \rho(q)\,q\rVert ,

which is $\lVert Kz - mz\rVert$ in the notation of Stakgold and Holst (Stakgold & Holst (2011), §6.5, (6.5.19), p. 401), who read it as the standard deviation of a mass distribution on the eigenvalue axis.

Four iterations are defined below. In each, $q_0$ is a unit vector chosen by the user.

Power iteration. $\displaystyle q_{k+1} = \frac{Aq_k}{\lVert Aq_k\rVert}$ .
Inverse iteration with shift $\sigma$ . Solve $(A - \sigma I)z_{k+1} = q_k$ and set $q_{k+1} = z_{k+1}/\lVert z_{k+1}\rVert$ . The shift $\sigma$ is fixed, and $A - \sigma I$ is factorized once (Solving linear systems numerically).
Rayleigh quotient iteration. The same, but with the shift recomputed at every step: $\sigma_k = \rho(q_k)$ , then $(A - \sigma_k I)z_{k+1} = q_k$ and $q_{k+1} = z_{k+1}/\lVert z_{k+1}\rVert$ .
The QR iteration. Put $A_0 = A$ . For $k \ge 1$ , factor $A_{k-1} = Q_kR_k$ with $Q_k$ unitary and $R_k$ upper triangular (Orthonormal bases and Gram-Schmidt) and multiply the factors back in the other order:

A_{k-1} = Q_kR_k, \qquad A_k = R_kQ_k .

Write $\bar Q_k = Q_1Q_2\cdots Q_k$ and $\bar R_k = R_kR_{k-1}\cdots R_1$ for the accumulated factors.

The angle $\theta(q, S) \in [0,\pi/2]$ between a unit vector $q$ and a subspace $S$ is defined by $\cos\theta = \lVert P_Sq\rVert$ , with $P_S$ the orthogonal projector onto $S$ ; so $\sin\theta = \lVert (I-P_S)q\rVert$ is the distance from $q$ to $S$ . Convergence "to an eigenvector" always means $\theta \to 0$ , never $q_k \to v$ : an iterate is determined only up to a scalar of modulus one, and for a negative eigenvalue the sign genuinely alternates.

Method(Power iteration)Original

Let $A$ be Hermitian with eigenvalues $\lambda_1,\dots,\lambda_n$ listed with multiplicity and ordered so that

\lvert\lambda_1\rvert > \lvert\lambda_2\rvert \ge \cdots \ge \lvert\lambda_n\rvert ,

and let $v_1,\dots,v_n$ be a corresponding orthonormal eigenbasis (Spectral theorem for Hermitian matrices). The strict first inequality forces $\lambda_1$ to be a simple eigenvalue. Let $q_0$ be a unit vector with $(q_0,v_1) \neq 0$ , let $\theta_k$ be the angle between $q_k$ and $\operatorname{span}\{v_1\}$ , and put $D = \max_i\lvert\lambda_1-\lambda_i\rvert$ . Then for every $k \ge 0$ :

$\displaystyle \tan\theta_k \;\le\; \Bigl\lvert\frac{\lambda_2}{\lambda_1}\Bigr\rvert^{\,k}\tan\theta_0$ , so $\theta_k \to 0$ geometrically with ratio $\lvert\lambda_2/\lambda_1\rvert$ ;
$\displaystyle \bigl\lvert\lambda_1 - \rho(q_k)\bigr\rvert \;\le\; D\,\sin^{2}\theta_k \;\le\; D\,\tan^{2}\theta_0\,\Bigl\lvert\frac{\lambda_2}{\lambda_1}\Bigr\rvert^{\,2k} ,$

so the eigenvector converges at the rate $\lvert\lambda_2/\lambda_1\rvert$ and the eigenvalue estimate at the square of that rate. If in addition $\lambda_1$ is the algebraically largest eigenvalue, then $\rho(q_k) \le \lambda_1$ for every $k$ , so the estimates approach $\lambda_1$ from below.

Idea of the proofOriginal

In the eigenbasis, multiplying by $A$ multiplies the $i$ th coordinate by $\lambda_i$ . After $k$ steps the coordinate along $v_1$ has been multiplied by $\lambda_1^k$ and every other coordinate by something no larger than $\lambda_2^k$ in modulus, so the ratio of what we do not want to what we do want is multiplied by $\lvert\lambda_2/\lambda_1\rvert$ each time. Normalising at every step changes nothing in that ratio; it only keeps the numbers inside the range of the floating-point format. The Rayleigh quotient does better than the vector because it is stationary at an eigenvector: the error in $\rho$ is a weighted average of the gaps $\lambda_1-\lambda_i$ with weights equal to the squares of the unwanted coordinates.

ProofOriginal

Normalising at each step only rescales, so $q_k = A^kq_0/\lVert A^kq_0\rVert$ by induction. Expand $q_0 = \sum_i c_iv_i$ in the orthonormal eigenbasis, so $\sum_i\lvert c_i\rvert^2 = 1$ and $c_1 = (q_0,v_1) \neq 0$ by hypothesis.

The angle. Since $A^kq_0 = \sum_i c_i\lambda_i^{\,k}v_i$ and the $v_i$ are orthonormal, the component of $A^kq_0$ along $v_1$ has modulus $\lvert c_1\rvert\,\lvert\lambda_1\rvert^k$ and the component orthogonal to $v_1$ has norm $\bigl(\sum_{i\ge2}\lvert c_i\rvert^2\lvert\lambda_i\rvert^{2k}\bigr)^{1/2}$ . Dividing the second by the first — both are positive, the first because $c_1\neq0$ and $\lambda_1\neq0$ — gives the tangent of the angle exactly:

\tan\theta_k = \frac{\bigl(\sum_{i\ge2}\lvert c_i\rvert^2\lvert\lambda_i\rvert^{2k}\bigr)^{1/2}}{\lvert c_1\rvert\,\lvert\lambda_1\rvert^{k}} .

By the ordering, $\lvert\lambda_i\rvert \le \lvert\lambda_2\rvert$ for every $i\ge2$ , so the numerator is at most $\lvert\lambda_2\rvert^{k}\bigl(\sum_{i\ge2}\lvert c_i\rvert^2\bigr)^{1/2}$ . Dividing, and recognising $\bigl(\sum_{i\ge2}\lvert c_i\rvert^2\bigr)^{1/2}/\lvert c_1\rvert$ as $\tan\theta_0$ , gives claim 1.

The Rayleigh quotient. Write $q_k = \sum_i d_iv_i$ with $\sum_i\lvert d_i\rvert^2 = 1$ . Hermiticity is used here: because the $v_i$ are orthonormal and $Av_i = \lambda_iv_i$ ,

\rho(q_k) = (Aq_k,q_k) = \sum_i \lambda_i\lvert d_i\rvert^2 .

Subtracting this from $\lambda_1 = \lambda_1\sum_i\lvert d_i\rvert^2$ and noting that the $i=1$ term cancels,

\lambda_1 - \rho(q_k) = \sum_{i\ge2}(\lambda_1-\lambda_i)\lvert d_i\rvert^2 .

Every $\lambda_i$ is real, again by Hermiticity, so this is a real number of modulus at most $D\sum_{i\ge2}\lvert d_i\rvert^2$ . That sum is exactly $\sin^2\theta_k$ , since $\lvert d_1\rvert = \cos\theta_k$ . The second inequality in claim 2 follows from claim 1 and $\sin\theta\le\tan\theta$ .

One-sidedness. If $\lambda_1 \ge \lambda_i$ for all $i$ , every term of the last display is non-negative, so $\lambda_1-\rho(q_k)\ge0$ . $\square$

Remark(A residual bracket: when to stop)AdaptedAdapted from Stakgold & Holst (2011), §6.5, Thm 6.5.5 with (6.5.17) and (6.5.19), pp. 401-402

The bounds above contain $\theta_0$ and $\lambda_2/\lambda_1$ , neither of which is known while the computation is running. The following bound uses only quantities already on hand, and it is where the residual $\eta(q)$ earns its place.

Lemma. Let $A$ be Hermitian and $q$ a unit vector. Then the closed interval $[\rho(q)-\eta(q),\ \rho(q)+\eta(q)]$ contains an eigenvalue of $A$ .

Proof. Write $q = \sum_id_iv_i$ in an orthonormal eigenbasis, $\sum_i\lvert d_i\rvert^2=1$ , and $\rho = \rho(q)$ . Then $Aq-\rho q = \sum_i(\lambda_i-\rho)d_iv_i$ , so by orthonormality $\eta(q)^2 = \sum_i(\lambda_i-\rho)^2\lvert d_i\rvert^2$ . If no eigenvalue lay within $\eta(q)$ of $\rho$ , every factor $(\lambda_i-\rho)^2$ with $d_i\neq0$ would exceed $\eta(q)^2$ , and the sum would exceed $\eta(q)^2\sum_i\lvert d_i\rvert^2 = \eta(q)^2$ , a contradiction. $\square$

Stakgold and Holst (Stakgold & Holst (2011), §6.5, Thm 6.5.5, p. 402) state and prove this for a positive compact symmetric operator on a Hilbert space, where $\rho$ and $\eta$ are the mean and the standard deviation of the mass distribution that the trial element $z$ puts on the eigenvalue axis, and they note that the result is Chebyshev's inequality in disguise. The proof above is that argument in finite dimensions. Two cautions carry over. First, the bracket names an eigenvalue, not the one you are chasing: to conclude that $\lambda_1 \in [\rho-\eta,\rho+\eta]$ one needs to know in addition that no other eigenvalue is that close, which is the point of the source's condition $m - \sigma > \lambda_2$ . Second, the bracket is honest but pessimistic by roughly one factor of the residual, because by part 2 of the statement above the actual error in $\rho$ is of order $\eta^2/\mathrm{gap}$ once $q$ is close.

Remark(Where the hypotheses are used, and the operator form of the same iteration)AdaptedAdapted from Stakgold & Holst (2011), §6.5, Exercise 6.5.2, p. 406

$\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ strictly. Used in step 1 to make the ratio in the numerator contract. It fails for $A = \operatorname{diag}(1,-1)$ and $q_0 = (1,1)^{T}/\sqrt2$ : then $q_k = (1,(-1)^k)^{T}/\sqrt2$ for every $k$ , the iteration oscillates for ever between two vectors, and $\rho(q_k) = 0$ , which is not an eigenvalue at all. For a real non-symmetric matrix the same failure occurs whenever the dominant eigenvalues are a complex conjugate pair, as for a rotation matrix. When $\lambda_1$ is a repeated eigenvalue the statement survives unchanged if $\lambda_2$ is read as the largest modulus strictly below $\lvert\lambda_1\rvert$ and $v_1$ is replaced by the whole $\lambda_1$ -eigenspace: the proof uses only that the unwanted coordinates are multiplied by something smaller.
$(q_0,v_1) \neq 0$ . Used to divide by $\lvert c_1\rvert$ . In exact arithmetic an iteration started orthogonal to $v_1$ stays orthogonal to it and converges to the dominant eigenvector of the rest of the spectrum. In floating point it does not: rounding injects a component along $v_1$ of relative size about the unit roundoff $u$ (Floating-point arithmetic and error), and from then on that component grows by $\lvert\lambda_1/\lambda_2\rvert$ per step, so the iteration turns towards $v_1$ after roughly $\ln(1/u)/\ln\lvert\lambda_1/\lambda_2\rvert$ additional steps. A random $q_0$ makes the event $c_1 = 0$ one of probability zero, which is why starting vectors are chosen at random in practice.
Hermiticity. Used twice in step 2 and nowhere in step 1. Claim 1 holds for any diagonalizable $A$ with the constant $\tan\theta_0$ replaced by a constant involving the condition number of the eigenvector matrix. Claim 2 does not: without an orthonormal eigenbasis the expansion $\rho = \sum_i\lambda_i\lvert d_i\rvert^2$ is false, and the Rayleigh quotient converges only at the first-order rate $\lvert\lambda_2/\lambda_1\rvert^k$ . Doubling the exponent is a reward for symmetry.
The normalisation. Nothing in the mathematics needs it; the floating-point format does. Without it the iterates grow or shrink like $\lvert\lambda_1\rvert^k$ and overflow or underflow within a few hundred steps.
The same iteration in a Hilbert space. Stakgold and Holst (Stakgold & Holst (2011), §6.5, Exercise 6.5.2, p. 406) set out the power iteration for a nonnegative symmetric Hilbert-Schmidt operator $K$ in a form that avoids normalising at all: from a trial function $f_0$ with $Kf_0 \neq 0$ they form $f_n = K^nf_0$ , define the moments $a_k = (f_i, f_{k-i})$ — a quantity that depends only on $k$ , not on the split — and state that the quotients $a_{k+1}/a_k$ increase monotonically and are bounded above by the largest eigenvalue $\lambda_1$ , with convergence to $\lambda_1$ whenever $f_0$ is not orthogonal to the first eigenfunction. Their hypothesis that $f_0$ is not orthogonal to that eigenfunction is our $c_1\neq0$ ; their monotonicity is a sharper statement than our part 2 and is bought with the extra hypothesis that $K$ is nonnegative. The example below exhibits both in the same computation.

Method(Inverse iteration with a shift, and the Rayleigh quotient iteration)Original

Let $A$ be Hermitian with eigenvalues $\lambda_1,\dots,\lambda_n$ and orthonormal eigenbasis $v_1,\dots,v_n$ , and let $\sigma\in\mathbb{R}$ not be an eigenvalue of $A$ .

The shifted inverse. $(A-\sigma I)^{-1}$ is Hermitian, has the same eigenvectors as $A$ , and has eigenvalues $(\lambda_i-\sigma)^{-1}$ . Consequently its dominant eigenvalue corresponds to the eigenvalue of $A$ nearest $\sigma$ .
Inverse iteration. Suppose some $\lambda_j$ is strictly nearest $\sigma$ , that is $\lvert\lambda_j-\sigma\rvert < \delta := \min_{i\neq j}\lvert\lambda_i-\sigma\rvert$ , and let $q_0$ be a unit vector with $(q_0,v_j)\neq0$ . Let $\theta_k$ be the angle between the inverse iterate $q_k$ and $\operatorname{span}\{v_j\}$ . Then

\tan\theta_k \;\le\; \Bigl(\frac{\lvert\lambda_j-\sigma\rvert}{\delta}\Bigr)^{k}\tan\theta_0 , \qquad \bigl\lvert\lambda_j - \rho(q_k)\bigr\rvert \;\le\; D\sin^2\theta_k ,

with $D = \max_i\lvert\lambda_j-\lambda_i\rvert$ . The convergence factor $\lvert\lambda_j-\sigma\rvert/\delta$ is at the user's disposal: it is small exactly when $\sigma$ is a good guess. 3. Cubic convergence of the Rayleigh quotient iteration. Let $q_k$ be a unit vector, $\sigma_k=\rho(q_k)$ , $\theta_k$ the angle from $q_k$ to $\operatorname{span}\{v_j\}$ and $\delta_k = \min_{i\neq j}\lvert\lambda_i-\sigma_k\rvert$ . If $\sigma_k$ is not an eigenvalue, then one step of inverse iteration with the shift $\sigma_k$ produces $q_{k+1}$ with

\tan\theta_{k+1} \;\le\; \frac{\lvert\lambda_j-\sigma_k\rvert}{\delta_k}\tan\theta_k \;\le\; \frac{D}{\delta_k}\,\sin^{2}\theta_k\,\tan\theta_k .

Once $\theta_k$ is small enough that $\delta_k$ is bounded below by half the gap between $\lambda_j$ and the rest of the spectrum, the error is therefore cubed at every step.

Idea of the proofOriginal

Everything here is the previous theorem applied to a different matrix. Inverting turns the eigenvalue nearest $\sigma$ — the one the power method is worst at finding — into the one of largest modulus, which the power method is best at finding (Riley, Hobson & Bence (2006), §8.13, p. 273). The cost is one linear solve per step instead of one matrix-vector product, and since the matrix $A-\sigma I$ is the same at every step it is factorized once. The Rayleigh quotient iteration then notices that a shift is only as good as the information behind it, and feeds the current eigenvalue estimate back in as the next shift. Each step gains twice over: the shift is better by the square of the angle, which is the Rayleigh quotient's own accuracy, and the step then reduces the angle in proportion to the shift error. Two times the error squared, once more linearly: the exponent is three.

ProofOriginal

Set $B = (A-\sigma I)^{-1}$ . Since $\sigma$ is real and $A$ is Hermitian, $A-\sigma I$ is Hermitian, and it is invertible because $\sigma$ is not an eigenvalue; the inverse of an invertible Hermitian matrix is Hermitian, as $(B^{*}) = ((A-\sigma I)^{*})^{-1} = B$ . From $Av_i = \lambda_iv_i$ we get $(A-\sigma I)v_i = (\lambda_i-\sigma)v_i$ , and applying $B$ to both sides and dividing by the non-zero scalar $\lambda_i-\sigma$ gives $Bv_i = (\lambda_i-\sigma)^{-1}v_i$ . This is the familiar statement that the inverse keeps the eigenvectors and reciprocates the eigenvalues (Riley, Hobson & Bence (2006), §8.13, p. 273) (Arfken, Weber & Harris (2013), Exercise 6.4.5, p. 318), applied to the shifted matrix. The map $t\mapsto(t-\sigma)^{-1}$ is largest in modulus at the $\lambda_i$ closest to $\sigma$ , which proves the last clause.
The inverse iteration is the power iteration for $B$ , since $q_{k+1}$ is $Bq_k$ normalised. $B$ is Hermitian with the same eigenvectors $v_i$ , and by hypothesis its eigenvalue $(\lambda_j-\sigma)^{-1}$ is strictly largest in modulus, with the next largest of modulus $1/\delta$ . The power iteration statement applied to $B$ gives

\tan\theta_k \le \Bigl(\frac{1/\delta}{1/\lvert\lambda_j-\sigma\rvert}\Bigr)^{k}\tan\theta_0 = \Bigl(\frac{\lvert\lambda_j-\sigma\rvert}{\delta}\Bigr)^{k}\tan\theta_0 .

The Rayleigh quotient bound is part 2 of that statement applied to $A$ itself — not to $B$ — with $v_j$ in the role of $v_1$ : its proof used only that $v_j$ is a unit eigenvector of a Hermitian matrix and that $\sin\theta_k$ is the norm of the part of $q_k$ orthogonal to it, both of which hold here.

The first inequality is a direct one-step estimate, and needs no assumption beyond the one made: that $\sigma_k$ is not an eigenvalue. Expand $q_k = \sum_id_iv_i$ in the orthonormal eigenbasis. One step of inverse iteration with the shift $\sigma_k$ replaces $q_k$ by the normalisation of $(A-\sigma_kI)^{-1}q_k = \sum_i(\lambda_i-\sigma_k)^{-1}d_iv_i$ , which is defined because no $\lambda_i$ equals $\sigma_k$ . A tangent is unchanged by scaling, so, taking $d_j\neq0$ (otherwise $\tan\theta_k=\infty$ and there is nothing to prove),

\tan\theta_{k+1} = \frac{\bigl(\sum_{i\neq j}\lvert d_i\rvert^{2}/\lvert\lambda_i-\sigma_k\rvert^{2}\bigr)^{1/2}}{\lvert d_j\rvert/\lvert\lambda_j-\sigma_k\rvert} \;\le\; \frac{\lvert\lambda_j-\sigma_k\rvert}{\delta_k}\,\frac{\bigl(\sum_{i\neq j}\lvert d_i\rvert^{2}\bigr)^{1/2}}{\lvert d_j\rvert} \;=\; \frac{\lvert\lambda_j-\sigma_k\rvert}{\delta_k}\tan\theta_k ,

where the inequality is $\lvert\lambda_i-\sigma_k\rvert\ge\delta_k$ for every $i\neq j$ . (This is the computation behind part 2, run for a single step; part 2 itself may not be invoked here, because its hypothesis that $\lambda_j$ is strictly nearest the shift is exactly what part 3 declines to assume. Without that hypothesis the estimate still holds, the factor merely failing to be smaller than $1$ .) For the second inequality, part 2 of the power iteration statement, applied to $A$ and $v_j$ , bounds the shift error itself: $\lvert\lambda_j-\sigma_k\rvert = \lvert\lambda_j-\rho(q_k)\rvert \le D\sin^2\theta_k$ . Substituting gives the stated bound. If $\theta_k$ is small enough that $\delta_k \ge \tfrac12\min_{i\neq j}\lvert\lambda_i-\lambda_j\rvert =: \tfrac12 g$ , then $\tan\theta_{k+1} \le (2D/g)\sin^2\theta_k\tan\theta_k \le (2D/g)\tan^3\theta_k$ , a cubic recursion, which contracts as soon as $\tan^2\theta_k < g/(2D)$ . $\square$

Method(The QR algorithm)Original

Let $A\in M_n(\mathbb{C})$ and let $A_k$ , $\bar Q_k$ , $\bar R_k$ be as in the definition, with every $R_k$ taken to have non-negative diagonal. Write $E_j = \operatorname{span}\{e_1,\dots,e_j\}$ for the span of the first $j$ standard basis vectors.

Similarity. $A_k = \bar Q_k^{*}A\,\bar Q_k$ for every $k$ . Every iterate therefore has the same eigenvalues, trace and determinant as $A$ (Arfken, Weber & Harris (2013), Exercise 6.4.1, p. 317), and if $A$ is Hermitian then so is every $A_k$ .
It is power iteration on a flag. $A^k = \bar Q_k\bar R_k$ . If $A$ is invertible, the first $j$ columns of $\bar Q_k$ are an orthonormal basis of $A^kE_j$ , simultaneously for $j=1,\dots,n$ .
Convergence. Suppose $A = X\Lambda X^{-1}$ is diagonalizable with $\Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_n)$ and

\lvert\lambda_1\rvert > \lvert\lambda_2\rvert > \cdots > \lvert\lambda_n\rvert > 0 ,

and suppose every leading principal submatrix of $X^{-1}$ is invertible, equivalently that $X^{-1} = LU$ with $L$ unit lower triangular and $U$ upper triangular (Solving linear systems numerically). Then there is a constant $C$ , depending on $A$ but not on $k$ , such that for each $j = 1,\dots,n-1$ ,

\lvert (A_k)_{mi}\rvert \;\le\; C\Bigl\lvert\frac{\lambda_{j+1}}{\lambda_j}\Bigr\rvert^{k} \quad\text{for all } i\le j < m , \qquad\text{and}\qquad (A_k)_{jj} \longrightarrow \lambda_j \ \text{ as } k\to\infty ,

the last also for $j=n$ . In particular the subdiagonal entries go to zero, the iterates become upper triangular in the limit, and the eigenvalues appear on the diagonal in order of decreasing modulus.

Idea of the proofOriginal

Run the power iteration on $e_1$ and you get the dominant eigenvector. Run it on $e_1$ and $e_2$ at once and both columns swing towards $v_1$ , so the second carries no new information — unless you orthogonalize against the first at every step, which is exactly what taking a QR factorization does. What survives in column $j$ is then the part of the dominant behaviour that the first $j-1$ columns have not already claimed, so the span of the first $j$ columns converges to the invariant subspace belonging to the $j$ largest eigenvalues, at the rate $\lvert\lambda_{j+1}/\lambda_j\rvert$ set by the gap at that position in the list. The two lines of the algorithm, "factor and multiply back in the other order", do this for all $j$ at once while keeping the matrix similar to $A$ . Once the span of the first $j$ columns is invariant, the block of $A_k$ below the leading $j\times j$ corner vanishes; that happening for every $j$ simultaneously is triangularity, and the diagonal of a triangular matrix is its spectrum.

ProofOriginal

1. Similarity. From $A_{k-1}=Q_kR_k$ and the unitarity of $Q_k$ we get $R_k = Q_k^{*}A_{k-1}$ , hence $A_k = R_kQ_k = Q_k^{*}A_{k-1}Q_k$ . Induction on $k$ , with $A_0=A$ and $\bar Q_0 = I$ , gives $A_k = \bar Q_k^{*}A\bar Q_k$ . A similarity transformation changes neither the characteristic polynomial nor therefore the eigenvalues, trace or determinant (Arfken, Weber & Harris (2013), Exercise 6.4.1, p. 317); and if $A^{*}=A$ then $A_k^{*} = \bar Q_k^{*}A^{*}\bar Q_k = A_k$ .

2. The power-iteration identity. For $k=1$ the claim $A^1 = \bar Q_1\bar R_1$ is the factorization $A = Q_1R_1$ . Assume $A^k = \bar Q_k\bar R_k$ . By part 1, $A\bar Q_k = \bar Q_kA_k$ , and $A_k = Q_{k+1}R_{k+1}$ by definition, so

A^{k+1} = A\,\bar Q_k\bar R_k = \bar Q_kA_k\bar R_k = \bar Q_kQ_{k+1}R_{k+1}\bar R_k = \bar Q_{k+1}\bar R_{k+1} .

For the second sentence, let $\bar q_1,\dots,\bar q_n$ be the columns of $\bar Q_k$ . Reading column $i$ of $A^k = \bar Q_k\bar R_k$ and using that $\bar R_k$ is upper triangular gives $A^ke_i = \sum_{m\le i}(\bar R_k)_{mi}\bar q_m$ , so $A^kE_j \subseteq \operatorname{span}\{\bar q_1,\dots,\bar q_j\}$ . If $A$ is invertible then so is $A^k$ , so $A^kE_j$ has dimension $j$ and the inclusion is an equality. The $\bar q_m$ are orthonormal because $\bar Q_k$ is unitary.

3. Convergence. Fix $j$ with $1\le j\le n-1$ and write $r_j = \lvert\lambda_{j+1}/\lambda_j\rvert < 1$ . Let $x_1,\dots,x_n$ be the columns of $X$ , so $Ax_i = \lambda_ix_i$ , and let $S_j = \operatorname{span}\{x_1,\dots,x_j\}$ , which is $A$ -invariant. Note $A$ is invertible, since $\lvert\lambda_n\rvert>0$ .

Step (a): an explicit basis of $A^kE_j$ . Using $A^k = X\Lambda^kX^{-1}$ and the hypothesis $X^{-1}=LU$ ,

A^k = X\Lambda^kLU = X\bigl(\Lambda^kL\Lambda^{-k}\bigr)\Lambda^kU .

Both $U$ and $\Lambda^k$ are invertible and upper triangular, so $\Lambda^kU$ maps $E_j$ onto $E_j$ and

A^kE_j = XN_kE_j, \qquad N_k := \Lambda^kL\Lambda^{-k}, \quad (N_k)_{im} = L_{im}\Bigl(\frac{\lambda_i}{\lambda_m}\Bigr)^{k} .

$N_k$ is unit lower triangular. Split it into blocks along the $j$ th row and column, $N_k = \begin{pmatrix}N^{11}&0\\N^{21}&N^{22}\end{pmatrix}$ with $N^{11}$ of size $j\times j$ . Then $N_kE_j$ is spanned by the columns of $\begin{pmatrix}N^{11}\\N^{21}\end{pmatrix}$ , and multiplying on the right by the invertible $(N^{11})^{-1}$ does not change that span, so

N_kE_j = \text{column span of } \begin{pmatrix}I_j\\ \Delta_k\end{pmatrix}, \qquad \Delta_k := N^{21}(N^{11})^{-1} .

Step (b): $\Delta_k$ is small. The entries of $N^{21}$ are $L_{im}(\lambda_i/\lambda_m)^k$ with $i>j\ge m$ ; since the moduli decrease strictly, $\lvert\lambda_i/\lambda_m\rvert \le \lvert\lambda_{j+1}/\lambda_j\rvert = r_j$ for all such $i,m$ , so every entry is at most $\ell\,r_j^{\,k}$ in modulus, where $\ell = \max_{i>m}\lvert L_{im}\rvert$ . The matrix $N^{11}$ is $I_j$ plus a strictly lower triangular matrix whose entries are bounded by $\ell$ uniformly in $k$ (each is $L_{im}$ times a number of modulus below one); a strictly lower triangular $j\times j$ matrix is nilpotent of index $j$ , so $(N^{11})^{-1} = \sum_{p<j}(I_j-N^{11})^{p}$ is bounded in norm uniformly in $k$ . Hence $\lVert\Delta_k\rVert \le c_1r_j^{\,k}$ for a constant $c_1$ independent of $k$ .

Step (c): $A^kE_j$ is close to $S_j$ , in both directions. By steps (a) and (b), $A^kE_j$ is spanned by the vectors $y_i = x_i + \sum_{m>j}(\Delta_k)_{mi}x_m$ for $i\le j$ . Let $u = \sum_{i\le j}c_iy_i$ be any element of $A^kE_j$ and put $u_0 = \sum_{i\le j}c_ix_i \in S_j$ . Then $u = X\binom{c}{\Delta_kc}$ and $u-u_0 = X\binom{0}{\Delta_kc}$ , so

\lVert u-u_0\rVert \le \lVert X\rVert\,\lVert\Delta_k\rVert\,\lvert c\rvert \quad\text{and}\quad \lVert u\rVert \ge \sigma_{\min}(X)\,\lvert c\rvert ,

whence $\lVert u-u_0\rVert \le \varepsilon_k\lVert u\rVert$ with $\varepsilon_k := \kappa_2(X)\lVert\Delta_k\rVert \le C_0r_j^{\,k}$ (Matrix norms and condition numbers). Exchanging the roles of $u$ and $u_0$ in the same two estimates — given $w = \sum_{i\le j}c_ix_i\in S_j$ , compare it with $w' = \sum_{i\le j}c_iy_i\in A^kE_j$ — gives $\operatorname{dist}(w, A^kE_j) \le \varepsilon_k\lVert w\rVert$ as well. So, writing $P_k$ and $P$ for the orthogonal projectors onto $A^kE_j$ and $S_j$ , both $\lVert (I-P)P_k\rVert$ and $\lVert (I-P_k)P\rVert$ are at most $\varepsilon_k$ , and from the identity $P_k-P = P_k(I-P) - (I-P_k)P$ together with $\lVert P_k(I-P)\rVert = \lVert (I-P)P_k\rVert$ (adjoints of each other, both projectors being self-adjoint),

\lVert P_k - P\rVert \le 2\varepsilon_k .

Step (d): the block below the corner. Let $i\le j<m$ . By part 2, $\bar q_i\in A^kE_j$ and $\bar q_m\perp A^kE_j$ , so $(A_k)_{mi} = (A\bar q_i,\bar q_m)$ is a component of $(I-P_k)A\bar q_i$ and is bounded by its norm. Write $\bar q_i = u_0+r$ with $u_0\in S_j$ and $\lVert r\rVert\le\varepsilon_k$ , as step (c) permits. Then $Au_0\in S_j$ because $S_j$ is invariant, so $\lVert (I-P_k)Au_0\rVert \le \varepsilon_k\lVert Au_0\rVert \le \varepsilon_k\lVert A\rVert(1+\varepsilon_k)$ by the second estimate of step (c), while $\lVert (I-P_k)Ar\rVert\le\lVert A\rVert\varepsilon_k$ . Adding, and assuming $\varepsilon_k\le1$ (true for large $k$ ),

\lvert (A_k)_{mi}\rvert \le \lVert (I-P_k)A\bar q_i\rVert \le 3\lVert A\rVert\varepsilon_k \le C r_j^{\,k} .

Enlarging $C$ to cover the finitely many small $k$ with $\varepsilon_k>1$ gives the stated bound for all $k$ .

Step (e): the diagonal. The trace of the leading $j\times j$ block of $A_k$ is $\sum_{i\le j}(A\bar q_i,\bar q_i) = \operatorname{tr}(P_kAP_k) = \operatorname{tr}(AP_k)$ . Since $\lvert\operatorname{tr}M\rvert\le n\lVert M\rVert$ and $\lVert A(P_k-P)\rVert\le2\varepsilon_k\lVert A\rVert$ by step (c),

\bigl\lvert \operatorname{tr}(AP_k) - \operatorname{tr}(AP)\bigr\rvert \le 2n\lVert A\rVert\varepsilon_k \xrightarrow[k\to\infty]{} 0 .

In an orthonormal basis whose first $j$ vectors span $S_j$ , invariance makes $A$ block upper triangular with leading block $B$ the matrix of $A$ restricted to $S_j$ , and $PAP$ is $B$ in the corner and zero elsewhere; since $S_j$ is spanned by eigenvectors for $\lambda_1,\dots,\lambda_j$ , $\operatorname{tr}(AP) = \operatorname{tr}B = \lambda_1+\cdots+\lambda_j$ . So the trace of the leading $j\times j$ block of $A_k$ tends to $\lambda_1+\cdots+\lambda_j$ for $j=1,\dots,n-1$ , and for $j=n$ it equals $\operatorname{tr}A = \lambda_1+\cdots+\lambda_n$ exactly, by part 1. Subtracting consecutive cases,

(A_k)_{jj} = \Bigl(\text{trace of leading } j\times j\Bigr) - \Bigl(\text{trace of leading } (j-1)\times(j-1)\Bigr) \longrightarrow \lambda_j

for every $j$ , the case $j=1$ needing only the first bracket. $\square$

Remark(Where the hypotheses are used)Original

Distinct moduli, in part 3. Step (b) needs $\lvert\lambda_i/\lambda_m\rvert\le r_j<1$ for $i>j\ge m$ , which is exactly the strict ordering. Without it the algorithm need not converge at all: for $A = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ the QR factorization is $Q=A$ , $R=I$ , so $A_1 = RQ = A$ and every iterate equals $A$ , whose off-diagonal entries never decay. Here $\lambda = \pm1$ have equal modulus. The same obstruction appears for every real matrix with a complex conjugate pair of eigenvalues, since a real iteration cannot converge to a real triangular matrix with non-real diagonal; the practical remedy is to stop at the real Schur form, in which such a pair survives as a $2\times2$ diagonal block, and to read its eigenvalues off directly.
The $LU$ factorization of $X^{-1}$ . Used only in step (a), to arrange that $A^kE_j$ is a graph over $S_j$ . Without it the diagonal can converge to the eigenvalues in the wrong order. Take $A = \operatorname{diag}(1,2)$ with the eigenvalues listed as $\lambda_1 = 2$ , $\lambda_2 = 1$ , so that $X = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ and $X^{-1}=X$ has a vanishing leading entry and no $LU$ factorization. The moduli are distinct, but $A$ is already upper triangular with $Q_1=I$ , so every iterate is $\operatorname{diag}(1,2)$ and $(A_k)_{11} = 1 \neq \lambda_1$ . Nothing is lost numerically — the eigenvalues are still on the diagonal — but the clean ordering is.
$\lvert\lambda_n\rvert>0$ . Used in part 2 to know that $A^kE_j$ has dimension $j$ . A singular $A$ is handled by shifting: the eigenvalues of $A+cI$ are $\lambda_i+c$ , and the iteration is run on the shifted matrix.
Hermiticity is not assumed in any part of this statement, unlike the two earlier ones. What Hermiticity buys is part 1's last clause: every iterate stays Hermitian, so only half the matrix need be stored and touched, and the tridiagonal form of the method below is available.
Non-negative diagonal in $R_k$ . A convenience, not a necessity: it makes the factorization unique for invertible $A_{k-1}$ (Orthonormal bases and Gram-Schmidt) and so makes "the" QR iteration well defined. Any other choice of phases changes each $A_k$ by a diagonal unitary conjugation, which alters no modulus in the statement.

Method(Computing eigenvalues in practice)Original

Decide first which eigenvalues you need; that decides the algorithm.

Locate the spectrum cheaply. Gershgorin's theorem — every eigenvalue $\lambda$ of $C=(c_{jk})$ satisfies $\lvert c_{jj}-\lambda\rvert \le \sum_{k\neq j}\lvert c_{jk}\rvert$ for at least one $j$ (Kreyszig (1978), §5.2, Problem 4, p. 313) — costs $O(n^2)$ and brackets every eigenvalue in a union of discs. Use it to choose shifts and to check an answer for gross error.
The dominant eigenvalue only, with $A$ large and sparse: power iteration. One matrix-vector product per step; cost per step is the cost of that product. Stop when $\eta(q_k)$ is below the tolerance, and report the bracket $[\rho-\eta,\rho+\eta]$ of the lemma above. Expect $\ln(1/\text{tol})/\ln\lvert\lambda_1/\lambda_2\rvert$ steps, so stop early and change method if the ratio is near one.
One eigenvalue near a known $\sigma$ , or an eigenvector for an eigenvalue you already have: inverse iteration. Factor $A-\sigma I$ once with pivoted $LU$ (Solving linear systems numerically), at cost $\tfrac23n^3$ for a dense matrix and far less for a banded one, then reuse the factors each step at $O(n^2)$ . Convergence factor $\lvert\lambda_j-\sigma\rvert/\delta$ : a shift good to three digits usually gives the eigenvector in three or four steps.
Refining to full accuracy: switch to the Rayleigh quotient iteration once $\theta$ is small. It cubes the error, but $A-\sigma_kI$ changes every step, so each step costs a fresh factorization. Two or three steps are enough; do not start from a poor guess, since the iteration is only locally convergent and can converge to an eigenvalue you did not want.
All the eigenvalues of a dense matrix: never iterate on the full matrix. First reduce $A$ by unitary similarities to upper Hessenberg form (zero below the first subdiagonal), or to tridiagonal form if $A$ is Hermitian. This costs $O(n^3)$ once and is exact in the sense that it terminates: counting a multiply and an add as one operation each, applying a Householder reflection to an $m\times p$ block costs $4mp$ , so the reduction costs about $\sum_k\bigl[4(n-k)^2+4n(n-k)\bigr] \approx \tfrac{10}{3}n^3$ in general and about $\tfrac43n^3$ when symmetry lets the two-sided update be done at once.
Then run shifted QR with deflation. The Hessenberg form is preserved by a QR step, and this is what makes the algorithm affordable: with Givens rotations a Hessenberg QR step costs about $6n^2$ operations — $n-1$ rotations acting on two rows of length $n-k+1$ at $6$ operations per column, about $3n^2$ in all, and the same again applying them from the right — against about $\tfrac43n^3$ for a QR factorization of a full matrix. A tridiagonal step costs $O(n)$ . Use the shift $\sigma_k = (A_k)_{nn}$ , or Wilkinson's shift (the eigenvalue of the trailing $2\times2$ block nearer to $(A_k)_{nn}$ ) when a real matrix may have complex eigenvalues. When $\lvert (A_k)_{n,n-1}\rvert$ falls below $u(\lvert(A_k)_{n-1,n-1}\rvert+\lvert(A_k)_{nn}\rvert)$ , accept $(A_k)_{nn}$ as an eigenvalue, deflate — delete the last row and column — and continue on the leading $(n-1)\times(n-1)$ block. With shifts this takes two or three steps per eigenvalue in practice, so the whole iteration costs $O(n^2)$ for a tridiagonal matrix and the $\tfrac43n^3$ reduction dominates.
Eigenvectors. Accumulate $\bar Q_k$ if you want them all, which multiplies the cost by several; for a few, run one or two steps of inverse iteration with the computed eigenvalue as the shift, which is cheaper and more accurate.
Check the answer. Compute $\eta$ for each computed pair and use the bracket. Compare $\sum_i\lambda_i$ with $\operatorname{tr}A$ and $\prod_i\lambda_i$ with $\det A$ , both invariant under the similarities used (Arfken, Weber & Harris (2013), Exercise 6.4.1, p. 317).

The older Jacobi method — repeated plane rotations, each chosen to annihilate one off-diagonal entry, with $\tan2\theta = 2a_{ij}/(a_{jj}-a_{ii})$ , each rotation reducing the sum of squares of the off-diagonal entries by $2a_{ij}^2$ (Arfken, Weber & Harris (2013), §6.4 and Exercise 6.4.9, pp. 314, 318) — is still the most transparent way to see that a symmetric matrix can be diagonalized by iteration, and is still used where accuracy on tiny eigenvalues matters. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §6.4, p. 314) describe it as superseded by more efficient though less transparent methods; step 6 is what superseded it.

Example(All three iterations on one

4\times4

matrix)Original

Take the real symmetric matrix

A = \begin{pmatrix} 6&2&1&0\\ 2&5&1&1\\ 1&1&4&2\\ 0&1&2&3\end{pmatrix} , \qquad \det(A-\lambda I) = (\lambda-5)\,(\lambda^3-13\lambda^2+43\lambda-36) .

The factor $\lambda-5$ is visible, and the cubic has three real roots, so the eigenvalues in decreasing order are

\lambda_1 = 8.382699939755810,\quad \lambda_2 = 5,\quad \lambda_3 = 3.326149857358553,\quad \lambda_4 = 1.291150202885637 ,

the irrational ones computed to the digits shown and checked against $\sum\lambda_i = \operatorname{tr}A = 18$ and $\prod\lambda_i = \det A = 180$ . $A$ is positive definite, and the ratios that govern the three algorithms are

\frac{\lambda_2}{\lambda_1} = 0.5964665,\qquad \frac{\lambda_3}{\lambda_2} = 0.6652300,\qquad \frac{\lambda_4}{\lambda_3} = 0.3881816 .

(a) Power iteration, by hand to begin with. Start at $q_0 = e_1$ ; the exact dominant eigenvector is $v_1 = (0.6579120, 0.5834795, 0.4006479, 0.2572641)^{T}$ , so $(q_0,v_1) = 0.658 \neq 0$ and the hypothesis holds. The unnormalised iterates $A^ke_1$ are integer vectors,

(6,2,1,0),\quad (41,23,12,4),\quad (304,213,120,59),\quad (2370,1852,1115,630),

and their first entries are the moments $a_k = (A^ke_1,e_1)$ of the Stakgold-Holst exercise: $a_0,\dots,a_6 = 1, 6, 41, 304, 2370, 19039, 155666$ . The quotients

\frac{a_1}{a_0}=6,\quad \frac{41}{6}=6.8333,\quad \frac{304}{41}=7.41463,\quad \frac{2370}{304}=7.79605,\quad \frac{19039}{2370}=8.03333,\quad \frac{155666}{19039}=8.17616

increase monotonically and stay below $\lambda_1 = 8.38270$ , exactly as that exercise asserts for a nonnegative symmetric operator.

Normalising at each step and recording the Rayleigh quotient, the angle to $v_1$ and the residual:

$k$	$\rho(q_k)$	$\lambda_1-\rho(q_k)$	$\sin\theta_k$	$\eta(q_k)$	bracket
5	8.367694690	$1.50\times10^{-2}$	$6.644\times10^{-2}$	$2.256\times10^{-1}$	$[8.14207,\ 8.59332]$
10	8.382615289	$8.47\times10^{-5}$	$5.002\times10^{-3}$	$1.692\times10^{-2}$	$[8.36569,\ 8.39954]$
15	8.382699457	$4.82\times10^{-7}$	$3.776\times10^{-4}$	$1.277\times10^{-3}$	$[8.38142,\ 8.38398]$
20	8.382699937	$2.75\times10^{-9}$	$2.851\times10^{-5}$	$9.644\times10^{-5}$	$[8.38260,\ 8.38280]$

Every bracket contains $\lambda_1 = 8.3826999$ , as the lemma guarantees, and every one is pessimistic by about a factor of $\eta$ , as predicted. The ratio of successive $\sin\theta_k$ at $k=20$ is $0.5964665$ and the ratio of successive Rayleigh errors is $0.3557723$ ; the predicted values are $\lambda_2/\lambda_1 = 0.5964665$ and $(\lambda_2/\lambda_1)^2 = 0.3557723$ . The explicit bounds also hold with room to spare: $\tan\theta_0 = 1.1447$ , so part 1 allows $\tan\theta_{20}\le 1.1447\times0.5964665^{20} = 3.72\times10^{-5}$ against the observed $2.85\times10^{-5}$ , and with $D = \lambda_1-\lambda_4 = 7.0915$ part 2 allows $5.76\times10^{-9}$ against the observed $2.75\times10^{-9}$ .

(b) Inverse iteration at an interior eigenvalue. Take $\sigma = 3.5$ , a guess at $\lambda_3$ . The distances are $\lvert\lambda_i-\sigma\rvert = 4.883,\ 1.5,\ 0.1739,\ 2.209$ , so $\lambda_3$ is strictly nearest, $\delta = 1.5$ , and the predicted convergence factor is $0.17385/1.5 = 0.11590$ . From $q_0 = e_1$ :

$k$	$\rho(q_k)$	$\rho(q_k)-\lambda_3$	$\sin\theta_k$
1	3.374449339	$4.83\times10^{-2}$	$1.607\times10^{-1}$
3	3.326156621	$6.76\times10^{-6}$	$2.032\times10^{-3}$
5	3.326149859	$1.23\times10^{-9}$	$2.716\times10^{-5}$

The observed ratio $\sin\theta_5/\sin\theta_4 = 2.716\times10^{-5}/2.347\times10^{-4} = 0.1157$ matches the predicted $0.1159$ . Five steps against the twenty of part (a), for an eigenvalue the power method cannot reach at all.

(c) Rayleigh quotient iteration. Start from $q_0 = (0,0,1,1)^{T}/\sqrt2$ , for which $\rho(q_0) = 5.5$ . The successive Rayleigh quotients and angles are

\begin{aligned} \rho_1 &= 5.018153235, &\sin\theta_1 &= 1.330\times10^{-1},\\ \rho_2 &= 4.999999414, &\sin\theta_2 &= 1.094\times10^{-3},\\ \rho_3 &= 5 - 1.6\times10^{-19}, &\sin\theta_3 &= 3.456\times10^{-10},\\ \rho_4 &= 5 - 1.6\times10^{-57}, &\sin\theta_4 &= 3.156\times10^{-29} . \end{aligned}

The iteration has found $\lambda_2 = 5$ . Cubing each angle predicts the next: $(1.094\times10^{-3})^3 = 1.31\times10^{-9}$ against $3.46\times10^{-10}$ , and $(3.456\times10^{-10})^3 = 4.13\times10^{-29}$ against $3.156\times10^{-29}$ . In double precision this would have reached full accuracy at step three; the computation above was done in 60-digit arithmetic so that the cubic law could be seen twice.

(d) The unshifted QR algorithm. The hypotheses of part 3 hold: the four moduli are distinct and non-zero, and the leading principal minors of $X^{-1} = X^{T}$ (the eigenvector matrix is orthogonal here) are $-0.6579$ , $-0.3369$ , $0.7739$ , $1$ , all non-zero. The subdiagonal entries of $A_k$ :

$k$	$\lvert(A_k)_{21}\rvert$	$\lvert(A_k)_{32}\rvert$	$\lvert(A_k)_{43}\rvert$
10	$1.692\times10^{-2}$	$6.605\times10^{-2}$	$8.256\times10^{-6}$
20	$9.644\times10^{-5}$	$1.123\times10^{-3}$	$6.307\times10^{-10}$
30	$5.497\times10^{-7}$	$1.905\times10^{-5}$	$4.900\times10^{-14}$

The ratios of successive values at $k=30$ are $0.5964665$ , $0.6652300$ and $0.3881816$ : each subdiagonal entry decays at exactly the ratio $\lvert\lambda_{j+1}/\lambda_j\rvert$ belonging to its own position, confirming the rates in part 3. The diagonal at $k=30$ is $(8.3826999,\ 5.0000000,\ 3.3261499,\ 1.2911502)$ , the eigenvalues in order of decreasing modulus.

Part 2 of the statement can be checked on the same run. The first column of $\bar Q_k$ is the power iterate from $e_1$ , so column one of $A_k$ below the diagonal should reproduce the power iteration of part (a): at $k=20$ its norm is $9.64399\times10^{-5}$ and $(A_{20})_{11} = 8.382699937$ , which are exactly the residual $\eta(q_{20})$ and the Rayleigh quotient $\rho(q_{20})$ in the table of part (a).

(e) Shifted QR. Repeat with the shift $\sigma_k = (A_k)_{44}$ at every step, that is, factor $A_k-\sigma_kI = Q_{k+1}R_{k+1}$ and set $A_{k+1} = R_{k+1}Q_{k+1}+\sigma_kI$ . Starting from $\lvert(A_0)_{43}\rvert = 2$ :

\lvert(A_k)_{43}\rvert = 7.51\times10^{-1},\ 6.58\times10^{-1},\ 2.30\times10^{-1},\ 3.54\times10^{-3},\ 1.12\times10^{-8},\ 3.41\times10^{-25} \quad (k=1,\dots,6),

with $(A_6)_{44} = 1.291150202885637$ , matching $\lambda_4$ to all sixteen digits shown. The last three entries show the cubic law again, $(3.54\times10^{-3})^3 = 4.4\times10^{-8}$ and $(1.12\times10^{-8})^3 = 1.4\times10^{-24}$ ; this is the Rayleigh quotient iteration of part (c), hidden in the corner of the matrix. Six steps and one deflation, against the thirty of part (d) for a result still only accurate to seven digits in the first subdiagonal.

Application(Normal modes, and why nobody computes the whole spectrum)Original

Discretise $-u'' = \lambda u$ on $(0,1)$ with $u(0)=u(1)=0$ by second differences on a grid of $N = 200$ interior points, $h = 1/201$ (Finite-difference approximations). The matrix is $h^{-2}$ times the tridiagonal matrix with $2$ on the diagonal and $-1$ beside it, and its eigenvalues are known exactly, $\lambda_k = 4h^{-2}\sin^2(k\pi h/2)$ , which makes it a fair test. The smallest is $\lambda_1 = 9.869403$ against the continuum value $\pi^2 = 9.869604$ , a relative error of $-2.0\times10^{-5}$ of the expected size $O(h^2)$ ; the largest is $\lambda_{200} = 1.6159\times10^{5}$ .

A physicist wants the low modes: they are the ones the continuum problem has, and the high ones are discretisation artefacts. The numbers show why the choice of algorithm is not a detail.

Power iteration is useless here. It converges to the largest eigenvalue, and does so at the ratio $\lambda_{199}/\lambda_{200} = 0.9998168$ . Reaching ten digits would take $\ln(10^{-10})/\ln(0.9998168) \approx 1.26\times10^{5}$ matrix-vector products, for an eigenvalue of no physical interest.
Inverse iteration is immediate. With $\sigma = 0$ the convergence factor is $\lambda_1/\lambda_2 = 0.2500153$ , so ten digits take $\ln(10^{-10})/\ln(0.25) \approx 17$ steps. Each step is one solve with a tridiagonal matrix, $O(N)$ work after a single $O(N)$ factorization (Solving linear systems numerically): the whole computation costs less than two steps of the power iteration cost on a dense matrix of the same size.
Shift-invert reaches the interior. For the tenth mode, put $\sigma = 985$ . The nearest eigenvalue is $\lambda_{10} = 984.952862$ at distance $0.0471$ and the next nearest is $\lambda_9$ at distance $186.9$ , so the convergence factor is $2.522\times10^{-4}$ : three steps shrink the eigenvector angle by $(2.522\times10^{-4})^{3} = 1.6\times10^{-11}$ , eleven digits in the eigenvector, and, by the squaring in part 2 of the power iteration statement, far more than that in the eigenvalue. ( $\lambda_{10}$ differs from the continuum $100\pi^2 = 986.9604$ by $-2.0\times10^{-3}$ in relative terms, a hundred times the error in $\lambda_1$ — the discretisation error grows with the mode number, which is the real limit on how far up the spectrum this grid can be trusted.)

The same three-way choice recurs wherever a matrix eigenvalue problem comes out of physics: the normal modes of a structure or a molecule (Generalized eigenproblems and normal modes), the band energies of a tight-binding Hamiltonian, the stationary distribution of a Markov chain (the classic use of the power method, since the distribution is the dominant eigenvector), and the stability of a fixed point, where only the sign of the largest real part is wanted. For a dense matrix of moderate size, the reduction plus shifted QR of the method above returns the whole spectrum for a small multiple of $n^3$ operations and is what every library routine does. For the very large sparse problems of quantum chemistry it is not available at any price: a dense $n^3$ is unthinkable and even storing $\bar Q_k$ is impossible, so one computes a few eigenvalues by Krylov subspace methods, built on the same matrix-vector product as the power iteration but extracting far more from it, exactly as the conjugate gradient method does for linear systems (Solving linear systems numerically). Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §6.4, p. 314 and footnote) give the scale that this makes possible: a configuration-interaction calculation that obtained eigenvalues and eigenvectors of a matrix of dimension beyond $10^9$ .