Numerical methods · The Lumenary

Lens

MotivationH

Numbers that are almost right

Nearly four thousand years ago a scribe in Mesopotamia pressed a small clay tablet — now called YBC 7289 and kept at Yale — with a square, its diagonals, and the number $1;24,51,10$ in base $60$ . That is $1 + \frac{24}{60} + \frac{51}{60^2} + \frac{10}{60^3} \approx 1.4142130$ , the ratio of a square's diagonal to its side, correct to about six decimal places. The true value $\sqrt 2 = 1.41421356\ldots$ is irrational: no finite expansion in any base is exact. Every practical computation with it is an approximation.

Modern computers face the same fact on a vast scale. A weather model, a bridge design or a vaccine-trial analysis performs billions of arithmetic operations on numbers stored with about sixteen significant digits, so almost every operation rounds. Usually that is harmless. Sometimes it is not: in 1982 the Vancouver Stock Exchange started a new index at $1000$ and recomputed it thousands of times a day, truncating to three decimals each time; within two years it read about $520$ while the correctly computed value was above $1000$ .

Numerical analysis is the mathematics of computing with approximations: designing methods that find roots, derivatives and integrals from finitely many arithmetic operations, and proving how large their errors can be. This lesson treats the three classical sources of error — rounding, truncation and the failure of hypotheses — through floating-point arithmetic, root finding, interpolation and differentiation, and numerical integration.

It assumes calculus of one variable (the mean value theorem and Taylor's theorem with remainder) and some programming.

IntuitionH

Three kinds of error, and how fast they shrink

Computing $\int_0^1 e^{-x^2}dx$ on a laptop involves three different approximations.

Modelling error: the integral itself may be only a model of something physical. Numerical analysis takes the mathematical problem as given.
Truncation (discretisation) error: we replace the integral by a finite sum, say the trapezoidal rule with $n$ strips. Even with exact arithmetic the sum differs from the integral, typically by about $C/n^2$ .
Rounding error: every operation in the sum is done in floating point and is rounded to about $16$ significant digits.

Truncation error shrinks as we work harder; rounding error tends to grow with the number of operations. Good methods balance the two, and good analysis predicts the balance. The key measure of a method is how fast its error shrinks: an error that is halved at each step (bisection) gains one correct binary digit per step, while an error that is squared at each step (Newton's method) doubles the number of correct digits per step:

10^{-1} \to 10^{-2} \to 10^{-4} \to 10^{-8} \to 10^{-16}.

DefinitionsH

Definitions

Definition— floating-point number

A (binary, normalised) floating-point number with precision $p$ is a number $\pm m \cdot 2^{e}$ where the significand $m$ is an integer with $2^{p-1} \le m < 2^p$ and the exponent $e$ is an integer in a fixed range. The IEEE 754 binary64 format ("double precision", Python's float) has $p = 53$ and exponents giving magnitudes from about $10^{-308}$ to $10^{308}$ , plus zero, subnormal numbers, $\pm\infty$ and NaN ("not a number").

Definition— machine epsilon and unit roundoff

Machine epsilon $\varepsilon$ is the gap between $1$ and the next larger floating-point number: $\varepsilon = 2^{1-p}$ , which is $2^{-52} \approx 2.2 \times 10^{-16}$ for binary64. The unit roundoff is $u = \varepsilon/2 = 2^{-53} \approx 1.1 \times 10^{-16}$ , the largest relative error made by rounding to nearest.

Definition— absolute and relative error

If $\hat x$ approximates $x$ , its absolute error is $|\hat x - x|$ and, for $x \ne 0$ , its relative error is $|\hat x - x|/|x|$ . A relative error of about $10^{-k}$ means about $k$ correct significant digits.

Definition— order of convergence

A sequence $x_n \to r$ with errors $e_n = |x_n - r|$ converges with order $p \ge 1$ if $e_{n+1} \le C e_n^{\,p}$ for some constant $C$ and all large $n$ (with $C < 1$ when $p = 1$ ). Order $1$ is linear, order $2$ quadratic. In practice one estimates $p \approx \log(e_{n+1}/e_n)/\log(e_n/e_{n-1})$ .

Definition— condition number

For a differentiable function $f$ with $f(x) \ne 0$ , the (relative) condition number $\kappa_f(x) = \left|\frac{x f'(x)}{f(x)}\right|$ measures how much a small relative change in $x$ is magnified in $f(x)$ . A problem with a large condition number is ill-conditioned: no method can solve it accurately from inexact data.

ProofH

Floating point and rounding error

Theorem— Proposition 1 (rounding)

If $x$ is a real number in the normal range of binary64 and $\mathrm{fl}(x)$ is the nearest binary64 number, then $\mathrm{fl}(x) = x(1 + \delta)$ with $|\delta| \le u = 2^{-53}$ .

Proof

Take $x > 0$ and let $2^{E} \le x < 2^{E+1}$ . In this interval the binary64 numbers are $m \cdot 2^{E - 52}$ for $2^{52} \le m < 2^{53}$ , spaced $2^{E-52}$ apart (and $2^{E+1}$ is also one). So the nearest one is within half the spacing: $|\mathrm{fl}(x) - x| \le 2^{E - 53} \le 2^{-53} x$ .

∎

IEEE 754 requires $+$ , $-$ , $\times$ , $\div$ and $\sqrt{\ }$ to be correctly rounded: the computed result is $\mathrm{fl}$ of the exact result. This gives the standard model $\mathrm{fl}(a \circ b) = (a \circ b)(1 + \delta)$ , $|\delta| \le u$ , on which error analysis is built. So floating-point arithmetic is not random or vague; it is exact arithmetic followed by one well-defined rounding. But it is not ordinary arithmetic either: $0.1$ is stored as $\frac{3602879701896397}{2^{55}}$ , so $0.1 + 0.2 \ne 0.3$ in floating point, and addition is not associative: $(1 + 2^{-53}) + 2^{-53} = 1$ but $1 + (2^{-53} + 2^{-53}) = 1 + 2^{-52}$ .

Theorem— Theorem 2 (cancellation)

Suppose $\hat x = x(1 + \delta_1)$ and $\hat y = y(1 + \delta_2)$ with $|\delta_1|, |\delta_2| \le u$ , and $x \ne y$ . Then

\big|(\hat x - \hat y) - (x - y)\big| \le u\,(|x| + |y|), \qquad \text{so} \qquad \frac{|(\hat x - \hat y) - (x - y)|}{|x - y|} \le u\,\frac{|x| + |y|}{|x - y|}.

Proof

$(\hat x - \hat y) - (x - y) = x\delta_1 - y\delta_2$ , and $|x\delta_1 - y\delta_2| \le |x|\,u + |y|\,u$ .

∎

The subtraction itself is exact or nearly so; the danger is that it exposes errors already present in $\hat x$ and $\hat y$ . If $x$ and $y$ agree in their first $k$ digits, the factor $(|x| + |y|)/|x - y|$ is about $10^{k}$ , and $k$ digits are lost. This is catastrophic cancellation, and $(|x| + |y|)/|x - y|$ is exactly the condition number of subtraction.

Example— Two cures for cancellation

For small $x$ , $\cos x$ is close to $1$ , and $(1 - \cos x)/x^2$ computed directly loses digits: at $x = 10^{-8}$ the computer finds $\cos x$ rounds to exactly $1$ and returns $0$ instead of $\tfrac12$ . The identity $1 - \cos x = 2\sin^2(x/2)$ removes the subtraction, and $2\sin^2(x/2)/x^2$ is accurate to full precision.
The roots of $x^2 - 10^8 x + 1 = 0$ are about $10^8$ and $10^{-8}$ . The formula $\frac{10^8 - \sqrt{10^{16} - 4}}{2}$ for the small root subtracts two nearly equal numbers and is wrong by $25\%$ . Computing the large root first and using the product of the roots, $x_1 x_2 = 1$ , gives the small root to full precision.

A computed answer is backward stable if it is the exact answer to a slightly perturbed problem (perturbations of relative size about $u$ ). Then its accuracy is limited only by the problem's condition number: error $\lesssim$ condition number $\times$ backward error. This way of thinking, due largely to James Wilkinson, separates what is the method's fault from what is the problem's.

ProofH

Finding roots: bisection and Newton's method

Bisection. Given a continuous $f$ on $[a, b]$ with $f(a) \le 0 \le f(b)$ , look at the midpoint $m = (a + b)/2$ . If $f(m) \le 0$ , replace $a$ by $m$ ; otherwise replace $b$ by $m$ . Repeat.

Theorem— Theorem 3 (bisection)

After $n$ steps the interval $[a_n, b_n]$ has length $(b - a)/2^n$ , still satisfies $f(a_n) \le 0 \le f(b_n)$ , and contains a root $c$ of $f$ . The midpoint $m_n$ satisfies $|m_n - c| \le (b - a)/2^{n+1}$ . So bisection always converges, linearly with rate $\tfrac12$ .

Proof

The length is halved at each step, and the sign condition is preserved by construction: we keep $[m, b_n]$ only when $f(m) \le 0$ , and $[a_n, m]$ only when $f(m) > 0$ . By the intermediate value theorem $f$ has a root $c \in [a_n, b_n]$ , and the midpoint of an interval is within half its length of every point in it.

∎

Bisection is slow — about $3.3$ steps per decimal digit — but it needs only continuity and a sign change, and it cannot fail. Its proof rests on the completeness of $\R$ : the nested intervals shrink to a point, which is a root.

Newton's method. Replace $f$ near the current guess $x_n$ by its tangent line, and take the tangent's zero as the next guess:

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

For $f(x) = x^2 - a$ this is $x_{n+1} = \frac12\left(x_n + \frac{a}{x_n}\right)$ — the method Heron of Alexandria described in the first century CE. (It is one plausible way the Babylonian value could have been found, but how the scribes actually computed it is not recorded.)

Theorem— Theorem 4 (local quadratic convergence of Newton's method)

Let $f$ be twice continuously differentiable on an open interval containing $r$ , with $f(r) = 0$ and $f'(r) \ne 0$ . Then there is $\delta > 0$ such that for every starting value with $|x_0 - r| \le \delta$ , the Newton iterates are defined, stay within $\delta$ of $r$ , converge to $r$ , and satisfy

|x_{n+1} - r| \le M\,|x_n - r|^2, \qquad M = \frac{\max_{|x - r| \le \delta}|f''(x)|}{2\min_{|x - r| \le \delta}|f'(x)|}.

Moreover $\frac{x_{n+1} - r}{(x_n - r)^2} \to \frac{f''(r)}{2f'(r)}$ .

Proof

Since $f'$ is continuous and $f'(r) \ne 0$ , there is $\delta_0 > 0$ with $f' \ne 0$ on $I_0 = [r - \delta_0, r + \delta_0]$ . Let $M_0$ be the quotient $M$ computed on $I_0$ (finite, since $f''$ and $1/f'$ are continuous on the closed interval), and put $\delta = \min(\delta_0, 1/(2M_0))$ (or $\delta = \delta_0$ if $M_0 = 0$ ). Then $M \le M_0$ on the smaller interval.

Suppose $|x_n - r| \le \delta$ . Taylor's theorem with the Lagrange remainder at $x_n$ gives a $\xi_n$ between $x_n$ and $r$ with

0 = f(r) = f(x_n) + f'(x_n)(r - x_n) + \tfrac12 f''(\xi_n)(r - x_n)^2.

Dividing by $f'(x_n) \ne 0$ and rearranging,

x_{n+1} - r = x_n - r - \frac{f(x_n)}{f'(x_n)} = \frac{f''(\xi_n)}{2f'(x_n)}\,(x_n - r)^2 .

Hence $|x_{n+1} - r| \le M|x_n - r|^2 \le M\delta\,|x_n - r| \le \tfrac12|x_n - r|$ . By induction every iterate stays within $\delta$ of $r$ , $|x_n - r| \le \delta/2^n \to 0$ , and the quadratic bound holds at every step. Finally $\xi_n \to r$ and $x_n \to r$ , so by continuity the ratio tends to $f''(r)/(2f'(r))$ .

∎

The hypotheses are all needed, and Newton's method can fail badly when they are not met or when the start is too far away.

Example— When Newton's method fails or slows down

No convergence from a bad start. For $f(x) = x^3 - 2x + 2$ and $x_0 = 0$ : $x_1 = 0 - \frac{2}{-2} = 1$ and $x_2 = 1 - \frac{1}{1} = 0$ . The iterates cycle $0, 1, 0, 1, \dots$ forever, although $f$ has a real root near $-1.77$ . For $f(x) = \arctan x$ , starts with $|x_0| > 1.39$ (roughly) make the iterates oscillate with growing size.
Not differentiable at the root. For $f(x) = \sqrt[3]{x}$ (real cube root), $f'(0)$ does not exist and $x_{n+1} = x_n - \frac{x_n^{1/3}}{\frac13 x_n^{-2/3}} = -2x_n$ : every start except $0$ diverges.
A double root. For $f(x) = (x - 1)^2$ , $f'(1) = 0$ and $x_{n+1} - 1 = \frac12(x_n - 1)$ : convergence, but only linear. (For a root of multiplicity $m$ , the modified step $x - m f(x)/f'(x)$ restores quadratic convergence.)

The secant method replaces $f'(x_n)$ by the slope $\frac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}}$ through the last two points. It needs no derivative, and under the hypotheses of Theorem 4 it converges with order $\varphi = (1 + \sqrt5)/2 \approx 1.618$ : its errors satisfy $e_{n+1} \approx C e_n e_{n-1}$ , and the exponents of $e_n$ follow the Fibonacci recurrence (see recurrences). Practical root finders, such as Brent's method, combine bisection's safety with the speed of secant-type steps.

ProofH

Interpolation and numerical differentiation

Given $n + 1$ distinct nodes $x_0, \dots, x_n$ and values $f(x_i)$ , there is exactly one polynomial $p$ of degree at most $n$ with $p(x_i) = f(x_i)$ for all $i$ : the Lagrange form is $p(x) = \sum_i f(x_i)\prod_{j \ne i}\frac{x - x_j}{x_i - x_j}$ , and uniqueness holds because the difference of two such polynomials has $n + 1$ roots and degree at most $n$ .

Theorem— Theorem 5 (interpolation error)

If $f$ has $n + 1$ continuous derivatives on $[a, b]$ , which contains the nodes, then for every $x \in [a, b]$ there is $\xi \in (a, b)$ with

f(x) - p(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x - x_i).

Proof

If $x$ is a node, both sides are $0$ . Otherwise let $w(t) = \prod_i (t - x_i)$ and choose the constant $K$ so that $g(t) = f(t) - p(t) - K w(t)$ vanishes at $t = x$ . Then $g$ has at least $n + 2$ zeros in $[a, b]$ (the nodes and $x$ ). By Rolle's theorem $g'$ has at least $n + 1$ zeros between them, $g''$ at least $n$ , …, and $g^{(n+1)}$ at least one, $\xi$ . Since $p^{(n+1)} = 0$ and $w^{(n+1)} = (n+1)!$ , we get $0 = f^{(n+1)}(\xi) - K(n+1)!$ , and $g(x) = 0$ is the claimed formula.

∎

More nodes need not help. For Runge's function $1/(1 + 25x^2)$ on $[-1, 1]$ , interpolants at equally spaced nodes oscillate wildly near the ends and diverge as $n \to \infty$ (Runge, 1901), because $\prod(x - x_i)$ is huge there and the derivatives of $f$ grow fast. Choosing Chebyshev nodes $x_i = \cos\frac{(2i + 1)\pi}{2n + 2}$ makes $\max|\prod(x - x_i)|$ as small as possible, namely $2^{-n}$ , and for functions as smooth as this one the interpolants converge.

Numerical differentiation. Taylor's theorem gives

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

(for the central difference, use the intermediate value theorem to combine the two third-derivative terms). So the truncation error of the forward difference is $O(h)$ and of the central difference $O(h^2)$ . But if each value of $f$ carries an error up to $\epsilon \approx u|f(x)|$ , the forward difference also carries a rounding error up to $2\epsilon/h$ , so its total error is at most about

E(h) = \frac{h}{2}M_2 + \frac{2\epsilon}{h}, \qquad \text{minimised at } h^* = 2\sqrt{\epsilon/M_2}, \text{ where } E(h^*) = 2\sqrt{\epsilon M_2}.

With $\epsilon \approx 10^{-16}$ the best step is $h \approx 10^{-8}$ , giving only about $8$ correct digits; smaller $h$ is worse. For the central difference the optimum is $h \approx \epsilon^{1/3} \approx 10^{-5}$ , with about $11$ digits. The third experiment shows both V-shaped error curves.

ProofH

Numerical integration

Divide $[a, b]$ into $n$ strips of width $h = (b - a)/n$ with points $x_k = a + kh$ . The composite trapezoidal rule and, for even $n$ , the composite Simpson rule are

T_h = h\Big(\tfrac12 f(x_0) + f(x_1) + \cdots + f(x_{n-1}) + \tfrac12 f(x_n)\Big), \qquad S_h = \frac{h}{3}\Big(f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + 4f(x_{n-1}) + f(x_n)\Big).

The trapezoidal rule integrates the piecewise-linear interpolant of $f$ ; Simpson's rule integrates a quadratic interpolant on each pair of strips.

Theorem— Theorem 6 (trapezoidal rule)

If $f$ is twice differentiable on $[a, b]$ with $|f''| \le M_2$ , then $|T_h - \int_a^b f| \le \frac{(b - a)h^2}{12}M_2$ .

Proof

Consider one strip, shifted to $[0, h]$ , and define $E(t) = \frac{t}{2}\big(f(0) + f(t)\big) - \int_0^t f(s)\,ds$ for $0 \le t \le h$ : the error of a single trapezoid of width $t$ . Then $E(0) = 0$ and, by the fundamental theorem of calculus,

E'(t) = \tfrac12\big(f(0) + f(t)\big) + \tfrac{t}{2}f'(t) - f(t) = \tfrac12\big(f(0) - f(t)\big) + \tfrac{t}{2}f'(t), \qquad E'(0) = 0,

E''(t) = -\tfrac12 f'(t) + \tfrac12 f'(t) + \tfrac{t}{2}f''(t) = \tfrac{t}{2}f''(t).

So $|E''(t)| \le \frac{M_2}{2}t$ . Integrating from $0$ (using $E'(0) = 0$ , then $E(0) = 0$ ): $|E'(t)| \le \frac{M_2}{4}t^2$ and $|E(t)| \le \frac{M_2}{12}t^3$ . Each of the $n$ strips therefore contributes at most $\frac{M_2}{12}h^3$ , and $n h^3 = (b - a)h^2$ .

∎

Theorem— Theorem 7 (Simpson's rule)

If $f$ is four times differentiable on $[a, b]$ with $|f^{(4)}| \le M_4$ and $n$ is even, then $|S_h - \int_a^b f| \le \frac{(b - a)h^4}{180}M_4$ .

Proof

Consider one pair of strips, shifted to $[-h, h]$ , and for $0 \le t \le h$ let $E(t) = \frac{t}{3}\big(f(-t) + 4f(0) + f(t)\big) - \int_{-t}^{t} f(s)\,ds$ . Differentiating three times (each time some terms cancel):

E'(t) = -\tfrac23\big(f(t) + f(-t)\big) + \tfrac43 f(0) + \tfrac{t}{3}\big(f'(t) - f'(-t)\big),

E''(t) = -\tfrac13\big(f'(t) - f'(-t)\big) + \tfrac{t}{3}\big(f''(t) + f''(-t)\big), \qquad E'''(t) = \tfrac{t}{3}\big(f'''(t) - f'''(-t)\big),

with $E(0) = E'(0) = E''(0) = 0$ . By the mean value theorem $f'''(t) - f'''(-t) = 2t f^{(4)}(\zeta)$ for some $\zeta$ , so $|E'''(t)| \le \frac23 M_4 t^2$ . Integrating three times from $0$ : $|E''(t)| \le \frac29 M_4 t^3$ , $|E'(t)| \le \frac{1}{18}M_4 t^4$ , $|E(t)| \le \frac{1}{90}M_4 t^5$ . There are $n/2$ pairs, each contributing at most $\frac{M_4}{90}h^5$ , and $\frac{n}{2}\cdot h^5 = \frac{(b - a)h^4}{2}$ .

∎

So for smooth integrands, halving $h$ divides the trapezoidal error by about $4$ and Simpson's by about $16$ . The bounds also explain degrees of exactness: the trapezoidal rule is exact for polynomials of degree $\le 1$ ( $f'' = 0$ ), Simpson's rule for degree $\le 3$ ( $f^{(4)} = 0$ ) — one degree more than its quadratic interpolants would suggest, by symmetry — but not for $x^4$ . Two warnings, both visible in the third experiment:

Hypotheses matter. For $\sqrt x$ on $[0, 1]$ , $f''$ is unbounded near $0$ ; both rules then converge only like $h^{1.5}$ .
Bounds are not estimates. For a smooth periodic function integrated over a full period, the trapezoidal rule converges faster than any power of $h$ (the error terms cancel; this follows from the Euler–Maclaurin formula), so it beats Simpson's rule there.

Combining $T_h$ and $T_{h/2}$ as $\frac{4T_{h/2} - T_h}{3}$ cancels the $h^2$ term of the error and gives exactly Simpson's rule; repeating the idea (Richardson extrapolation) gives Romberg integration.

Formal mathematicsF

Numerical analysis in Lean

Lean has a hardware type Float, but its operations are opaque to proofs: Lean can run 0.1 + 0.2 but cannot reason about it. So we formalise a model: binary64 numbers as real numbers of the form $m \cdot 2^e$ , and rounding as "a relative error of at most $u$ ". First, machine epsilon — no binary64 number lies strictly between $1$ and $1 + 2^{-52}$ :

IsBinary64Lean 4Formally verified

/-- A positive **binary64 number** (an IEEE 754 "double"), ignoring the limits on the exponent:
`m · 2ᵉ` with an integer *significand* `m` of exactly 53 bits, `2⁵² ≤ m < 2⁵³`, and an integer
exponent `e`. -/
def IsBinary64 (x : ℝ) : Prop :=
  ∃ m e : ℤ, 2 ^ 52 ≤ m ∧ m < 2 ^ 53 ∧ x = m * (2 : ℝ) ^ e