15.1 Numerical analysis

A computer can only add, multiply and round, so it replaces hard functions by polynomials through sample points. The error of that replacement has an exact formula, equally spaced points can make it explode (Runge), Chebyshev points tame it, and rounding adds a second, separate kind of error.

A computer can add, subtract, multiply, divide and compare numbers. That is all. So how does it compute $\sin 2.3$ , the area under $e^{-x^2}$ , or $\sqrt 2$ ? The answer in every case is to replace the hard function by a polynomial, which needs only multiplication and addition, and then to ask two questions. How far is the polynomial from the real function? And how much does it matter that every addition and multiplication is rounded? This chapter answers the first with an exact error formula (and shows the surprise that more nodes can make things worse), and the second with a look at how a computer stores numbers.

Before you start

0.6: a polynomial of degree $n$ has at most $n$ roots. 1.4: Rolle’s theorem, the mean value theorem and Taylor’s theorem. 1.5: the integral. 2.7: a matrix is invertible exactly when its determinant is non-zero. 3.5: machine epsilon and the condition number. 4.2: Newton’s method as a contraction.

Quick checkA non-zero polynomial has degree at most 4. At most how many roots can it have?

Quick checkg is differentiable, and g(0) = g(1) = g(2) = 0. What can you say about g'?

Quick checkA square matrix has a non-zero determinant. What does that tell you about Mx = y?

The picture

Pick $n + 1$ points on the $x$ -axis (the nodes), read off the height of a function $f$ above each, and draw the one polynomial of degree at most $n$ through those $n + 1$ points. That polynomial is the interpolating polynomial. The picture below does this for $f(x) = 1/(1 + 25x^2)$ on $[-1, 1]$ , with the nodes equally spaced to begin with. Slide the degree $n$ up from $2$ to $30$ . For small $n$ the orange polynomial hugs the purple curve. Around $n = 8$ it starts to swing at the ends of the interval, and by $n = 20$ the swing is enormous: more nodes, a worse fit. Now switch the nodes to Chebyshev nodes, which crowd toward the ends. The fit now improves as $n$ grows. Finally choose “Drag the heights yourself” and pull one dot: the polynomial through nine points can leave the dots far behind between them.

Two questions come out of the picture. Why is the error between the nodes what it is (the error formula), and why do some node choices behave so much better than others (the Chebyshev nodes)?

The idea, in words

A polynomial of degree $n$ has $n + 1$ coefficients, so it can be made to pass through $n + 1$ points, and (as we prove) in only one way. Polynomials are easy: a computer evaluates them with additions and multiplications, and we can differentiate and integrate them exactly. So replace $f$ by the polynomial $p_n$ through its values at $n + 1$ nodes. What is the error $f(x) - p_n(x)$ ? At the nodes it is $0$ . Between them it depends on two things: how wiggly $f$ is, measured by its $(n+1)$ -st derivative, and how the nodes are placed, measured by the product $w(x) = (x - x_0)\cdots(x - x_n)$ . The error formula says

f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,(x - x_0)(x - x_1)\cdots(x - x_n).

It is the Taylor remainder of 1.4 with the single centre $a$ replaced by $n + 1$ different centres, and it has the same proof. Everything else in the chapter is either a consequence (quadrature rules are integrals of interpolants) or a warning (the second factor can be huge, and the computer rounds).

Why· Two knobs: the function and the nodes

You do not control $f$ , so you control the nodes. The factor $w(x)$ is the only part that depends on them. Equally spaced nodes make $w$ small in the middle and large near the ends. Chebyshev nodes spread the size of $w$ evenly, and we prove that no choice of nodes makes the largest $|w|$ smaller. That is the whole story of Runge’s phenomenon in two lines.

Where this came from· Interpolation: Newton, Waring, Lagrange, Runge

Astronomers and table makers needed values between the entries of a table. Isaac Newton worked out how to pass a polynomial through any number of points using differences; the method appears as Lemma V of Book III of the Principia (1687) and was published in full as Methodus Differentialis in 1711. Edward Waring published the formula we now call Lagrange’s in 1779. Joseph-Louis Lagrange gave it in 1795, in the fifth of his Leçons élémentaires at the École Normale, and it carries his name. It is natural to expect that more points give a better fit. Carl Runge showed in 1901 that this can fail: he interpolated $1/(1 + 25x^2)$ at equally spaced points on $[-1, 1]$ and found that the error grows without bound with the degree. The cure, nodes that cluster at the ends, is built on polynomials that Pafnuty Chebyshev introduced in 1854 while studying the linkages of steam engines. Using their roots as interpolation nodes came later.

The chapter answers six questions. (1) Does a polynomial through the points exist, and is it unique? (2) How big is its error? (3) How do we place nodes to make the error small, and what goes wrong with equal spacing? (4) How do we compute the polynomial efficiently, and update it when a node is added? (5) How do we integrate numerically, and with what error? (6) What does rounding do, and how fast does Newton’s method converge?

Definitions and theorems

Existence and uniqueness

Definition A

Interpolating polynomial

Let $x_0, x_1, \dots, x_n$ be $n + 1$ distinct real numbers (the nodes) and $y_0, \dots, y_n$ any real numbers. An interpolating polynomial is a polynomial $p$ of degree at most $n$ with $p(x_j) = y_j$ for $j = 0, \dots, n$ . For a function $f$ we take $y_j = f(x_j)$ and write $p_n$ for the interpolating polynomial of $f$ .

In one lineGiven distinct nodes $x_0, \dots, x_n$ and values $y_0, \dots, y_n$ , a polynomial $p$ of degree at most $n$ with $p(x_j) = y_j$ for every $j$ .

Why it existsIt is the simplest way to turn a table of values into a formula that can be evaluated, differentiated and integrated between the table entries.

Definition A

Lagrange basis polynomial

For distinct nodes $x_0, \dots, x_n$ the Lagrange basis polynomial for node $j$ is

\ell_j(x) = \prod_{\substack{k = 0 \\ k \ne j}}^{n} \frac{x - x_k}{x_j - x_k}.

It is a product of $n$ linear factors, so it has degree $n$ . At $x = x_j$ every factor is $1$ , so $\ell_j(x_j) = 1$ . At any other node $x_i$ the factor with $k = i$ is $0$ , so $\ell_j(x_i) = 0$ .

In one line $\ell_j(x) = \prod_{k \ne j} \dfrac{x - x_k}{x_j - x_k}$ : the degree- $n$ polynomial that equals $1$ at the node $x_j$ and $0$ at all the other nodes.

Why it existsIt is a building block: scaling it by y_j fixes the value at x_j without disturbing any other node, so adding the pieces solves the whole problem at once.

Theorem A Apostol 15.6; Lagrange 1795

Existence and uniqueness of the interpolating polynomial

Let $x_0, \dots, x_n$ be distinct and $y_0, \dots, y_n$ arbitrary. Then $p(x) = \sum_{j=0}^{n} y_j\,\ell_j(x)$ is an interpolating polynomial, and any polynomial of degree at most $n$ that agrees with it at the $n + 1$ nodes is equal to it.

In one lineFor distinct nodes there is exactly one polynomial of degree at most $n$ through the given $n+1$ points, and it is $p(x) = \sum_j y_j\,\ell_j(x)$ .

Why it is trueExistence is an explicit formula, so we can compute with it. Uniqueness says that every method (Lagrange, Newton, solving a linear system) finds the same polynomial, so the error formula applies to all.

ProofApostol 15.6

The idea: Existence: the formula works because each basis polynomial is a switch that is on at one node and off at the others. Uniqueness: the difference of two interpolants has too many roots.

p = \sum_j y_j\ell_j

has degree at most

n

A sum of polynomials of degree

n

times constants has degree at most

n

p(x_i) = \sum_j y_j\ell_j(x_i) = y_i

for each

i

, so

p

interpolates.

x = x_i

the basis polynomials are

\ell_j(x_i) = 1

j = i

and

0

otherwise.

Uniqueness: let

r = p - q

Suppose

q

also has degree at most

n

and

q(x_i) = y_i

for all

i

. The difference

r = p - q

has degree at most

n

r

has

n + 1

roots.

r(x_i) = y_i - y_i = 0

at each of the

n + 1

distinct nodes.

r = 0

, that is

q = p

A non-zero polynomial of degree at most

n

has at most

n

roots (0.6), so

r

cannot be non-zero.

Worked example the system solved by elimination in 2.6

Find the polynomial of degree at most 2 through the points (1, 2), (2, 3) and (4, 9). Use it to estimate the value at x = 3.

Solution.

Write the three basis polynomials. $\ell_0(x) = \dfrac{(x-2)(x-4)}{(1-2)(1-4)} = \dfrac{(x-2)(x-4)}{3}$ , $\ \ell_1(x) = \dfrac{(x-1)(x-4)}{(2-1)(2-4)} = -\dfrac{(x-1)(x-4)}{2}$ , $\ \ell_2(x) = \dfrac{(x-1)(x-2)}{(4-1)(4-2)} = \dfrac{(x-1)(x-2)}{6}$ . For $\ell_j$ put the other nodes in the numerator, and divide by the same product evaluated at $x_j$ , so that $\ell_j(x_j) = 1$ .
Combine. $p(x) = 2\ell_0(x) + 3\ell_1(x) + 9\ell_2(x)$ . The values at the nodes are the weights.
Expand. $p(x) = \tfrac23(x^2 - 6x + 8) - \tfrac32(x^2 - 5x + 4) + \tfrac32(x^2 - 3x + 2) = \tfrac23 x^2 - x + \tfrac73$ . Collect the $x^2$ terms: $\tfrac23 - \tfrac32 + \tfrac32 = \tfrac23$ . The $x$ terms: $-4 + \tfrac{15}2 - \tfrac92 = -1$ . The constants: $\tfrac{16}3 - 6 + 3 = \tfrac73$ .
Check the nodes. $p(1) = \tfrac23 - 1 + \tfrac73 = 2$ , $p(2) = \tfrac83 - 2 + \tfrac73 = 3$ , $p(4) = \tfrac{32}3 - 4 + \tfrac73 = 9$ . Any interpolating polynomial must pass through the data.
Evaluate at 3. $p(3) = 6 - 3 + \tfrac73 = \tfrac{16}3$ . Direct from the formula. Check with the basis: $\ell_0(3) = -\tfrac13$ , $\ell_1(3) = 1$ , $\ell_2(3) = \tfrac13$ , so $p(3) = -\tfrac23 + 3 + 3 = \tfrac{16}3$ .

Answer. $p(x) = \tfrac13(2x^2 - 3x + 7)$ , the same polynomial that Gaussian elimination found in 2.6, and $p(3) = \tfrac{16}3 \approx 5.33$ . By uniqueness the two methods have to agree.

The same problem is a system of linear equations. Writing $p(x) = c_0 + c_1x + \cdots + c_nx^n$ , the $n + 1$ conditions $p(x_j) = y_j$ say $Vc = y$ with the following matrix. This is the link with 2.7: there the determinant of this matrix was computed for particular numbers, and the pattern is a theorem.

Definition A

Vandermonde matrix

For nodes $x_0, \dots, x_n$ the Vandermonde matrix is $V$ with $V_{jk} = x_j^{\,k}$ , for $j, k = 0, \dots, n$ . Row $j$ is $(1, x_j, \dots, x_j^n)$ . Then $(Vc)_j = c_0 + c_1x_j + \cdots + c_nx_j^n = p(x_j)$ , so the coefficients $c$ of the interpolating polynomial solve $Vc = y$ .

In one lineThe $(n+1) \times (n+1)$ matrix whose row $j$ is $(1, x_j, x_j^2, \dots, x_j^n)$ . Its determinant is the product of all differences $x_j - x_i$ with $i < j$ .

Why it existsIt is the coefficient matrix of the interpolation problem in the monomial basis. It is invertible exactly when the nodes are distinct, which gives uniqueness a second proof.

Theorem A Vandermonde; Apostol 15.6 and exercises

The Vandermonde determinant

For nodes $x_0, \dots, x_n$ , $\ \det V = \prod_{0 \le i < j \le n}(x_j - x_i)$ .

In one line $\det V = \prod_{0 \le i < j \le n} (x_j - x_i)$ . So $V$ is invertible exactly when the nodes are distinct.

Why it is trueIt says again, in the language of 2.7, that distinct nodes are exactly what is needed for the interpolation problem to have a unique answer.

Proofstandard

The idea: Induction on n. Replace the last node by a variable t. The determinant becomes a polynomial in t of degree at most n that vanishes at the earlier nodes, so it is a known constant times a product of linear factors.

Base case

n = 1

For

n = 1

\det\begin{bmatrix}1 & x_0\\1 & x_1\end{bmatrix} = x_1 - x_0

, the product with the single pair

(0, 1)

Treat the last node as a variable

t

Assume the formula for

n - 1

and look at

n

. Let

D(t)

be the determinant of

V

with the last row replaced by

(1, t, \dots, t^n)

, so

D(x_n) = \det V

D

has degree at most

n

, with

t^n

-coefficient

\det V_{n-1}

Expand the determinant along the last row (2.7). Each entry

t^k

is multiplied by a constant cofactor, so

D

is a polynomial in

t

of degree at most

n

. The coefficient of

t^n

is the cofactor of the corner entry, which is the determinant of the Vandermonde matrix of

x_0, \dots, x_{n-1}

(sign

(-1)^{n+n} = +1

D(x_i) = 0

for

i = 0, \dots, n-1

t = x_i

for

i < n

then the last row equals row

i

, and a determinant with two equal rows is

0

(2.7).

D(t) = \det V_{n-1}\cdot (t - x_0)\cdots(t - x_{n-1})

A polynomial of degree at most

n

with the

n

distinct roots

x_0, \dots, x_{n-1}

is divisible by

(t - x_0)\cdots(t - x_{n-1})

(factor theorem), and the quotient must be a constant, the coefficient of

t^n

\det V = \prod_{0 \le i < j \le n-1}(x_j - x_i)\ \cdot\ \prod_{i=0}^{n-1}(x_n - x_i) = \prod_{0 \le i < j \le n}(x_j - x_i)

Put

t = x_n

and use the inductive hypothesis

\det V_{n-1} = \prod_{0 \le i < j \le n-1}(x_j - x_i)

So $\det V \ne 0$ exactly when all the differences $x_j - x_i$ are non-zero, that is, when the nodes are distinct. By 2.7, $Vc = y$ then has exactly one solution $c$ , a second proof that the interpolating polynomial exists and is unique.

Your answerWhat is the determinant of the Vandermonde matrix for the nodes 0, 2 and 5?

Quick checkWhy is it dangerous to compute an interpolating polynomial by solving Vc = y for the coefficients when the degree is high?

The error formula (Apostol 15.8)

Theorem A Apostol 15.8

Interpolation error formula

Let $f \in C^{n+1}[a, b]$ , let $x_0, \dots, x_n \in [a, b]$ be distinct, and let $p_n$ be the interpolating polynomial of degree at most $n$ with $p_n(x_j) = f(x_j)$ . Write $w(t) = (t - x_0)(t - x_1)\cdots(t - x_n)$ . Then for every $x \in [a, b]$ there is a point $\xi$ in the open interval $(a, b)$ (it depends on $x$ ) with

f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,w(x).

In one lineIf $f$ has $n+1$ continuous derivatives on $[a, b]$ and $p_n$ interpolates it at distinct nodes in $[a, b]$ , then for each $x$ in $[a, b]$ there is a point $\xi$ in $(a, b)$ with $f(x) - p_n(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{j=0}^{n}(x - x_j)$ .

Why it is trueIt says the error is the next derivative, divided by (n+1)!, times a polynomial that vanishes at every node. It is Taylor's remainder with n+1 centres instead of one, and it is proved the same way.

ProofApostol 15.8; the same trick as the Taylor remainder in 1.4

The idea: Build a function g that is zero at the n+1 nodes and also at x. Rolle's theorem, applied again and again, pushes the zeros down to a zero of the (n+1)-st derivative of g, and that zero is the xi.

Assume

x

is not a node.

x

is one of the nodes, then

f(x) = p_n(x)

and

w(x) = 0

, so both sides are

0

and any

\xi

works.

Let

K = \dfrac{f(x) - p_n(x)}{w(x)}

Then

w(x) \ne 0

, so we can divide.

K

is the unknown error constant: the number that makes the equation

f(x) - p_n(x) = K\,w(x)

true at this

x

Define

g(t) = f(t) - p_n(t) - K\,w(t)

This is a function of

t

with

x

fixed. It has

n+1

continuous derivatives because

f

does and

p_n

w

are polynomials.

g

0

at the

n+1

nodes and at

x

: it has

n + 2

distinct zeros in

[a, b]

At a node

x_j

f(x_j) = p_n(x_j)

and

w(x_j) = 0

. At

t = x

g(x) = f(x) - p_n(x) - K\,w(x) = 0

by the choice of

K

g'

has at least

n + 1

distinct zeros in

(a, b)

Rolle's theorem (1.4) on each of the

n+1

gaps between consecutive zeros of

g

gives a zero of

g'

strictly inside the gap.

g^{(n+1)}

has at least

1

zero. Call it

\xi

Repeat: each application loses one zero. After

k

applications there are at least

n + 2 - k

distinct zeros of

g^{(k)}

g^{(n+1)}(t) = f^{(n+1)}(t) - K\,(n+1)!

p_n

has degree at most

n

, so its

(n+1)

-st derivative is

0

. And

w

t^{n+1}

plus lower powers, so

w^{(n+1)}(t) = (n+1)!

K = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}

, and by the definition of

K

this is the formula.

Put

t = \xi

, where

g^{(n+1)}(\xi) = 0

, and solve for

K

Why· The picture behind the proof

Look at the interpolation picture at the top of the chapter. At each of the $n + 1$ nodes the orange curve and the purple one cross, so their difference is zero there. A smooth function that is zero at $n + 2$ places must wiggle at least $n + 1$ times (once between each pair of zeros), and its derivative must wiggle at least $n$ times, and so on: the $(n+1)$ -st derivative still has a zero. The proof is this count. It also explains why the error depends on $f^{(n+1)}$ : if $f$ were a polynomial of degree $n$ then $f^{(n+1)} = 0$ , the error is zero, and interpolation reproduces polynomials exactly.

Proof workshopRebuild the proof of the interpolation error formula

Tap the lines in the order a proof would write them. Each line can only use what is already above it.

Two consequences. First, a bound: if $|f^{(n+1)}| \le M_{n+1}$ on $[a, b]$ then

|f(x) - p_n(x)| \le \frac{M_{n+1}}{(n+1)!}\,\max_{[a,b]}|w|.

Second, for equally spaced nodes $w$ is explicit enough to bound.

Corollary A standard; Burden and Faires 3.1

Error bound for equally spaced nodes

Let nodes be $x_j = a + jh$ for $j = 0, \dots, n$ , with $h = (b - a)/n$ , and $|f^{(n+1)}| \le M_{n+1}$ on $[a, b]$ . Then for all $x \in [a, b]$ , $|f(x) - p_n(x)| \le \dfrac{M_{n+1}}{4(n+1)}\,h^{n+1}$ .

In one lineWith nodes $x_j = a + jh$ , $h = (b-a)/n$ , we have $|w(x)| \le \tfrac14 h^{n+1} n!$ on $[a, b]$ , so $|f(x) - p_n(x)| \le \dfrac{M_{n+1}\,h^{n+1}}{4(n+1)}$ .

Why it is trueIt turns the error formula into a number you can compute before interpolating, and it shows the error shrinks like h to the n+1 if f is smooth.

Proofstandard

The idea: Take x between two neighbouring nodes. The two nearest factors are at most h squared over 4. Each other factor is at most a whole number of widths h.

|(x - x_j)(x - x_{j+1})| \le h^2/4

Let

x \in [x_j, x_{j+1}]

for some

j \in \{0, \dots, n-1\}

. The product of the two nearest factors

(x - x_j)(x_{j+1} - x)

has fixed sum

h

, and a product of two numbers with a fixed sum is largest when they are equal (each

h/2

The factors to the left of

x_j

multiply to at most

h^j\,(j+1)!

For

i < j

|x - x_i| \le x_{j+1} - x_i = (j + 1 - i)h

. Multiplying over

i = 0, \dots, j-1

gives

h^j (j+1)(j)\cdots 2

The factors to the right of

x_{j+1}

multiply to at most

h^{n-j-1}\,(n-j)!

For

i > j + 1

|x - x_i| \le x_i - x_j = (i - j)h

. Multiplying over

i = j+2, \dots, n

gives

h^{n-j-1}\,2 \cdot 3 \cdots (n-j)

|w(x)| \le \tfrac14 h^{n+1}\,(j+1)!\,(n-j)!

Multiply the three bounds. The powers of

h

add up to

2 + j + (n - j - 1) = n + 1

(j+1)!\,(n-j)! \le n!

, hence

|w(x)| \le \tfrac14 h^{n+1}n!

(j+1)!\,(n-j)!\big/n! = (j+1)\big/\binom nj

. For

j = 0

this is

1

, and for

1 \le j \le n-1

we have

\binom nj \ge n \ge j + 1

, so the ratio is at most

1

|f(x) - p_n(x)| \le \dfrac{M_{n+1}}{(n+1)!}\cdot\dfrac{h^{n+1}n!}{4} = \dfrac{M_{n+1}h^{n+1}}{4(n+1)}

Put this in the error bound and use

n!/(n+1)! = 1/(n+1)

Worked example

How many equally spaced nodes on [0, 1] guarantee an error below 10⁻⁶ when interpolating eˣ? How big is the actual error?

Solution.

Bound the derivative. Every derivative of $e^x$ is $e^x \le e$ on $[0, 1]$ , so $M_{n+1} = e$ . The bound needs only a maximum of $|f^{(n+1)}|$ over the interval.
Write the bound. With $n$ intervals, $h = 1/n$ , and the error is at most $\dfrac{e}{4(n+1)\,n^{n+1}}$ . The corollary with $M_{n+1} = e$ and $h^{n+1} = n^{-(n+1)}$ .
Try small $n$ . $n = 5$ : $\dfrac{e}{4 \cdot 6 \cdot 5^6} \approx 7.2 \times 10^{-6}$ , too big. $n = 6$ : $\dfrac{e}{4 \cdot 7 \cdot 6^7} \approx 3.5 \times 10^{-7}$ , below $10^{-6}$ . The bound falls extremely fast, because $n^{n+1}$ grows faster than exponentially.
Compare with the truth. Computing the interpolants and the largest error on a fine grid gives $2.7 \times 10^{-6}$ for $n = 5$ and $1.2 \times 10^{-7}$ for $n = 6$ . The bound is safe (the true errors are below it) but not sharp. Here $n = 6$ is also the first that works in truth.

Answer. Seven nodes ( $n = 6$ ). The guaranteed error is $3.5 \times 10^{-7}$ ; the actual error is $1.2 \times 10^{-7}$ . Smooth functions are approximated spectacularly well by interpolation, provided the nodes are well placed.

Quick checkIn the error formula, the point ξ depends on x. What does that mean for how you use it?

Runge’s phenomenon and Chebyshev nodes (Apostol 15.14 to 15.17)

The error formula has two factors. For $f(x) = 1/(1 + 25x^2)$ the first factor is large: expanding $\frac{1}{1 + 25x^2} = 1 - 25x^2 + 625x^4 - \cdots$ (geometric series, 1.6) shows the Taylor coefficient of $x^{2m}$ is $\pm 25^m = \pm 5^{2m}$ , so $|f^{(2m)}(0)|/(2m)! = 5^{2m}$ and $M_{n+1}/(n+1)!$ is at least $5^{n+1}$ when $n+1$ is even. The derivatives grow because $f$ has poles at $x = \pm i/5$ in the complex plane (a pole is a point where the denominator is $0$ , here $1 + 25x^2 = 0$ ), only $0.2$ away from the interval. The second factor, $\max|w|$ for equally spaced nodes, is at most about $(2/e)^n \approx 0.74^n$ by the corollary and Stirling’s estimate $n! \approx \sqrt{2\pi n}\,(n/e)^n$ . The product of the two factors behaves like $5^n \cdot 0.74^n \approx 3.7^n$ , which grows, so the bound gives no guarantee. Runge found in 1901 that the actual error, not only the bound, grows with $n$ near the ends of the interval. The table of actual errors below (largest error on $[-1, 1]$ ) shows the damage, and what Chebyshev nodes do.

Degree $n$	Equally spaced	Chebyshev nodes
6	0.617	0.264
10	1.916	0.109
14	7.195	0.047
18	29.19	0.0225
30	2388	0.0021

The remedy changes only $w$ . We look for nodes that make the largest value of $|w|$ on $[-1, 1]$ as small as possible. The answer comes from a family of polynomials that Chebyshev found.

Definition A

Chebyshev polynomial

The Chebyshev polynomials are defined by $T_0(x) = 1$ , $T_1(x) = x$ and the recurrence $T_{n+1}(x) = 2x\,T_n(x) - T_{n-1}(x)$ for $n \ge 1$ . The first few: $T_2 = 2x^2 - 1$ , $T_3 = 4x^3 - 3x$ , $T_4 = 8x^4 - 8x^2 + 1$ .

In one line $T_0 = 1$ , $T_1 = x$ , and $T_{n+1}(x) = 2x\,T_n(x) - T_{n-1}(x)$ . The key property is $T_n(\cos\theta) = \cos(n\theta)$ .

Why it existsIt is the polynomial version of 'cos(n theta)': oscillating between -1 and 1 exactly n+1 times on [-1, 1], which is the most evenly spread a polynomial of its size can be.

Theorem A Apostol 15.15

Properties of the Chebyshev polynomials

For $n \ge 1$ :

$T_n(\cos\theta) = \cos(n\theta)$ for all $\theta$ .
$T_n$ has degree exactly $n$ , with leading coefficient $2^{n-1}$ .
The roots of $T_n$ are $x_k = \cos\dfrac{(2k+1)\pi}{2n}$ , $k = 0, \dots, n-1$ : $n$ distinct numbers in $(-1, 1)$ .
$|T_n(x)| \le 1$ for $|x| \le 1$ , and $T_n\big(\cos\frac{k\pi}{n}\big) = (-1)^k$ for $k = 0, \dots, n$ .

In one line $T_n(\cos\theta) = \cos n\theta$ ; $T_n$ has degree $n$ and leading coefficient $2^{n-1}$ (for $n \ge 1$ ); it has $n$ distinct roots $\cos\frac{(2k+1)\pi}{2n}$ in $(-1, 1)$ ; $|T_n| \le 1$ on $[-1, 1]$ and $T_n(\cos\frac{k\pi}n) = (-1)^k$ .

Why it is trueThese give the nodes (the roots) and the extra fact needed for the minimal property: T_n swings fully from -1 to 1 and back, n times.

ProofApostol 15.15

The idea: Part 1 by induction from the identity cos((n+1)t) + cos((n-1)t) = 2 cos t cos(nt). The rest are read off from cosines.

Part 1, base cases.

T_0(\cos\theta) = 1 = \cos 0

and

T_1(\cos\theta) = \cos\theta

Inductive step.

The addition formulas give

\cos((n+1)\theta) + \cos((n-1)\theta) = 2\cos\theta\cos n\theta

. By induction

T_n(\cos\theta) = \cos n\theta

and

T_{n-1}(\cos\theta) = \cos(n-1)\theta

, so

T_{n+1}(\cos\theta) = 2\cos\theta\cos n\theta - \cos(n-1)\theta = \cos(n+1)\theta

Part 2 by induction on the recurrence.

T_1

has degree

1

and coefficient

1 = 2^0

. If

T_n

has leading term

2^{n-1}x^n

then

2xT_n

has leading term

2^nx^{n+1}

, and

T_{n-1}

has lower degree.

Part 3.

T_n(x_k) = \cos(n\theta_k) = \cos\frac{(2k+1)\pi}{2} = 0

with

\theta_k = \frac{(2k+1)\pi}{2n}

and

x_k = \cos\theta_k

. For

k = 0, \dots, n - 1

the angles

\theta_k

are

n

distinct numbers in

(0, \pi)

, and cosine is one-to-one there, so the

x_k

are

n

distinct numbers. A polynomial of degree

n

has at most

n

roots, so these are all of them.

Part 4.

For

|x| \le 1

write

x = \cos\theta

; then

|T_n(x)| = |\cos n\theta| \le 1

. At

\theta = k\pi/n

\cos(n\theta) = \cos k\pi = (-1)^k

Definition A

Chebyshev nodes

For interpolation of degree $n$ on $[-1, 1]$ the Chebyshev nodes are the $n + 1$ roots of $T_{n+1}$ : $x_k = \cos\dfrac{(2k+1)\pi}{2n+2}$ , $k = 0, \dots, n$ . (On a general interval $[a, b]$ , move them by $x \mapsto \frac{a+b}{2} + \frac{b-a}{2}x$ .) Picture $n+1$ equally spaced points around the upper half of a circle, and drop each down to the diameter.

In one lineThe $n+1$ roots of $T_{n+1}$ : $x_k = \cos\dfrac{(2k+1)\pi}{2n+2}$ for $k = 0, \dots, n$ . They are the projections onto the $x$ -axis of equally spaced points on a semicircle, so they crowd toward the ends.

Why it existsThey are the interpolation nodes that make the factor w(x) as small as possible, uniformly over the interval, which is the whole point of the next theorem.

Theorem A Apostol 15.16; Chebyshev 1854

Chebyshev's minimal property

Let $q$ be a monic polynomial of degree $m \ge 1$ (leading coefficient $1$ ). Then $\max_{-1 \le x \le 1}|q(x)| \ge 2^{1-m}$ . The monic polynomial $T_m/2^{m-1}$ has maximum exactly $2^{1-m}$ .

In one lineAmong all monic polynomials of degree $m \ge 1$ , none has a smaller largest absolute value on $[-1, 1]$ than $2^{1-m}$ , and $T_m/2^{m-1}$ achieves it.

Why it is trueWith nodes at the roots of T_(n+1), the polynomial w(x) is exactly T_(n+1)/2^n, so no other choice of n+1 nodes gives a smaller maximum of |w|. Chebyshev nodes are optimal for the error bound.

ProofApostol 15.16

The idea: Suppose some monic q does better. Subtract it from the Chebyshev polynomial; the difference has too low a degree to change sign as often as it must.

The Chebyshev polynomial achieves

2^{1-m}

By part 2 of the properties,

\tilde T = T_m/2^{m-1}

is monic of degree

m

, and by part 4,

|\tilde T| \le 2^{1-m}

with

\tilde T(\cos\frac{k\pi}m) = (-1)^k 2^{1-m}

. So its maximum is exactly

2^{1-m}

Assume

q

is better.

Suppose, for contradiction, that a monic

q

of degree

m

has

|q(x)| < 2^{1-m}

for all

x \in [-1, 1]

r = \tilde T - q

has degree at most

m - 1

\tilde T

and

q

are both monic of degree

m

, so the

x^m

terms cancel.

r(x_0) > 0

r(x_1) < 0

r(x_2) > 0

, and so on, for the

m + 1

points

x_0 > x_1 > \cdots > x_m

x_k = \cos\frac{k\pi}m

\tilde T(x_k) = (-1)^k2^{1-m}

and

|q(x_k)| < 2^{1-m}

, so

r(x_k)

has the sign of

(-1)^k

r

has a root in each of

m

disjoint intervals: at least

m

roots.

Intermediate value theorem on each of the

m

gaps between consecutive

x_k

r

is continuous and changes sign across each gap.

Contradiction.

A non-zero polynomial of degree at most

m - 1

has at most

m - 1

roots, and

r

is non-zero because

r(x_0) \ne 0

Proof workshopRebuild Chebyshev's minimal property

Tap the lines in the order a proof would write them. Each line can only use what is already above it.

With the Chebyshev nodes $x_0, \dots, x_n$ , the polynomial $w(x) = (x - x_0)\cdots(x - x_n)$ is monic of degree $n + 1$ with the same roots as $T_{n+1}$ , so $w = T_{n+1}/2^n$ and $\max|w| = 2^{-n}$ . By the minimal property no other choice of nodes does better.

Corollary A Apostol 15.17

Error bound for Chebyshev nodes

If $f \in C^{n+1}[-1, 1]$ with $|f^{(n+1)}| \le M_{n+1}$ and $p_n$ interpolates at the Chebyshev nodes, then $|f(x) - p_n(x)| \le \dfrac{M_{n+1}}{2^n\,(n+1)!}$ for $x \in [-1, 1]$ . On a general interval $[a, b]$ , moving the nodes as described above multiplies $w$ by $\big(\tfrac{b-a}{2}\big)^{n+1}$ , so the bound becomes $\dfrac{M_{n+1}}{(n+1)!}\Big(\dfrac{b-a}{2}\Big)^{n+1}2^{-n}$ .

In one lineWith the Chebyshev nodes on $[-1, 1]$ , $|f(x) - p_n(x)| \le \dfrac{M_{n+1}}{2^n\,(n+1)!}$ , and no choice of $n+1$ nodes makes the factor $\max|w|$ smaller than $2^{-n}$ .

Why it is trueIt is the clean answer to 'where should the nodes go'. The factor 2^(-n) shrinks geometrically while the equally spaced factor only shrinks like 0.74^n.

Worked example

Interpolate eˣ on [−1, 1] at the three Chebyshev nodes (n = 2) and bound the error. Compare with three equally spaced nodes.

Solution.

Find the nodes. Roots of $T_3 = 4x^3 - 3x$ : $x = 0$ and $x = \pm\frac{\sqrt3}2 \approx \pm 0.866$ . Or from the formula $\cos\frac{(2k+1)\pi}{6}$ for $k = 0, 1, 2$ : $\cos 30^\circ$ , $\cos 90^\circ$ , $\cos 150^\circ$ .
Bound $M_3$ . $f''' = e^x \le e$ on $[-1, 1]$ . The same maximum of the third derivative for both choices of nodes.
Chebyshev bound. $\dfrac{M_3}{2^2\cdot 3!} = \dfrac{e}{24} \approx 0.113$ . The corollary with $n = 2$ .
Equally spaced bound. The nodes $-1, 0, 1$ have $w(x) = x(x^2 - 1)$ , with $\max|w| = \dfrac{2}{3\sqrt3} \approx 0.385$ (at $x = \pm\frac1{\sqrt3}$ ). So the bound is $\dfrac{e}{6}\cdot 0.385 \approx 0.174$ . Compare $0.385$ with the Chebyshev value $2^{-2} = 0.25$ .
Actual errors. Computing on a fine grid: $0.0565$ with Chebyshev nodes and $0.0785$ with equally spaced nodes. Both below their bounds; Chebyshev is better in the bound and in the truth.

Answer. The Chebyshev bound is $e/24 \approx 0.113$ against $0.174$ for equally spaced nodes; the true largest errors are $0.057$ and $0.079$ . The advantage is modest at $n = 2$ and becomes decisive at large $n$ .

Quick checkWhere do Chebyshev nodes sit on [−1, 1]?

Python labRunge's phenomenon with a plot

import numpy as np
import matplotlib.pyplot as plt

f = lambda x: 1 / (1 + 25 * x**2)

def interpolate(nodes, x):
    """Lagrange form: p(x) = sum over j of f(x_j) * l_j(x)."""
    p = np.zeros_like(x)
    for j, xj in enumerate(nodes):
        lj = np.ones_like(x)
        for k, xk in enumerate(nodes):
            if k != j:
                lj *= (x - xk) / (xj - xk)
        p += f(xj) * lj
    return p

x = np.linspace(-1, 1, 2001)
cheb = lambda n: np.cos((2 * np.arange(n + 1) + 1) * np.pi / (2 * n + 2))
print("degree  max error (equally spaced)  max error (Chebyshev)")
for n in (6, 10, 14, 18):
    e1 = abs(interpolate(np.linspace(-1, 1, n + 1), x) - f(x)).max()
    e2 = abs(interpolate(cheb(n), x) - f(x)).max()
    print(f"{n:<7d} {e1:<27.4f} {e2:.4f}")

plt.plot(x, f(x), "k", label="f")
plt.plot(x, interpolate(np.linspace(-1, 1, 15), x), label="14 equally spaced")
plt.plot(x, interpolate(cheb(14), x), label="14 Chebyshev")
plt.ylim(-1, 2); plt.legend()

Try this. The equally spaced errors are 0.617, 1.916, 7.195 and 29.19, growing by a factor of 3 to 4 every 4 degrees, while the Chebyshev errors fall from 0.264 to 0.0225. Change 25 to 1 in f: now both errors collapse towards 0 (0.0068 and 0.0021 at n = 6, below 0.0001 for both at n = 18). The poles of 1/(1 + c x^2) are at plus or minus i over the square root of c, and a pole close to the interval is what hurts. Then try f = np.abs: equally spaced nodes diverge (0.18 at n = 6 to 32.0 at n = 18) while Chebyshev nodes converge slowly (0.087 to 0.0315).

Even Chebyshev nodes are not a universal cure. A theorem of Georg Faber (1914) says that for every fixed scheme of nodes, one set of nodes for each degree $n$ , there is a continuous function whose interpolating polynomials do not converge uniformly; this includes the Chebyshev nodes (label B). What the nodes do buy is convergence for every function smooth enough, and, in practice, stability: the software library Chebfun is built on interpolation at Chebyshev nodes, and NumPy offers it as numpy.polynomial.chebyshev.chebinterpolate.

Newton’s divided differences (Apostol 15.10)

The Lagrange form has a drawback: if you add one more data point, every $\ell_j$ changes and you start again. Newton’s form builds the polynomial one node at a time, so that adding a node adds one term.

Definition A

Divided difference

For distinct nodes $x_0, \dots, x_k$ , the divided difference $f[x_0, \dots, x_k]$ is the leading coefficient (the coefficient of $x^k$ ) of the interpolating polynomial of degree at most $k$ through $(x_0, f(x_0)), \dots, (x_k, f(x_k))$ . For a single node, $f[x_0] = f(x_0)$ .

In one line $f[x_0, \dots, x_k]$ is the coefficient of $x^k$ in the interpolating polynomial of $f$ at $x_0, \dots, x_k$ . It satisfies $f[x_0, \dots, x_k] = \dfrac{f[x_1, \dots, x_k] - f[x_0, \dots, x_{k-1}]}{x_k - x_0}$ .

Why it existsIt is the number that corrects the interpolant when a new node is added, and it behaves like a derivative divided by a factorial.

Theorem A Newton 1687; Apostol 15.10

Newton's divided-difference formula

Let $x_0, \dots, x_n$ be distinct. Then

$p_n(x) = \displaystyle\sum_{k=0}^{n} f[x_0, \dots, x_k]\ \prod_{i=0}^{k-1}(x - x_i)$ (the empty product, for $k = 0$ , is $1$ ).
$f[x_0, \dots, x_k] = \dfrac{f[x_1, \dots, x_k] - f[x_0, \dots, x_{k-1}]}{x_k - x_0}$ for $k \ge 1$ .
If $f \in C^n$ on an interval containing the nodes, then $f[x_0, \dots, x_n] = \dfrac{f^{(n)}(\xi)}{n!}$ for some $\xi$ in that interval.

In one lineWith $p_k$ the interpolant through the first $k+1$ nodes, $p_n(x) = \sum_{k=0}^{n} f[x_0, \dots, x_k]\,(x - x_0)\cdots(x - x_{k-1})$ , and the divided differences satisfy the quotient recursion in the definition.

Why it is trueIt gives an algorithm: fill a triangle of differences, read the coefficients off the top edge, and add a node by adding one row.

Proofstandard

The idea: Part 1: the difference between two consecutive interpolants vanishes at the earlier nodes, so it is a multiple of their product. Part 2: build the higher interpolant out of two lower ones and use uniqueness. Part 3: Rolle.

Part 1. The difference

p_k - p_{k-1}

has the

k

roots

x_0, \dots, x_{k-1}

p_k - p_{k-1}

has degree at most

k

(both do), and at each node

x_i

with

i < k

both equal

f(x_i)

, so it vanishes there.

p_k(x) = p_{k-1}(x) + f[x_0, \dots, x_k]\,(x - x_0)\cdots(x - x_{k-1})

A polynomial of degree at most

k

with

k

known distinct roots is a constant times

(x - x_0)\cdots(x - x_{k-1})

. The constant is the coefficient of

x^k

;

p_{k-1}

has no

x^k

term, so the constant is the leading coefficient of

p_k

, which is

f[x_0, \dots, x_k]

by definition.

The sum for

p_n

follows.

Add up these equations for

k = 1, \dots, n

and start from

p_0 = f(x_0) = f[x_0]

Combine two lower interpolants.

Part 2. Let

q_a

interpolate at

x_1, \dots, x_k

and

q_b

x_0, \dots, x_{k-1}

, both of degree at most

k - 1

. Define

q(x) = \dfrac{(x - x_0)\,q_a(x) - (x - x_k)\,q_b(x)}{x_k - x_0}

, of degree at most

k

q

interpolates

f

at all of

x_0, \dots, x_k

x_0

q = \frac{0 - (x_0 - x_k)f(x_0)}{x_k - x_0} = f(x_0)

. At

x_k

q = \frac{(x_k - x_0)f(x_k)}{x_k - x_0} = f(x_k)

. At a middle node

x_i

: both

q_a

and

q_b

equal

f(x_i)

, so

q = f(x_i)\frac{(x_i - x_0) - (x_i - x_k)}{x_k - x_0} = f(x_i)

f[x_0, \dots, x_k] = \dfrac{f[x_1, \dots, x_k] - f[x_0, \dots, x_{k-1}]}{x_k - x_0}

By uniqueness (the existence and uniqueness theorem),

q = p_k

. Compare the coefficients of

x^k

q_a

and

q_b

have

x^{k-1}

coefficients

f[x_1, \dots, x_k]

and

f[x_0, \dots, x_{k-1}]

, which become the

x^k

coefficient of

q

after multiplying by

x

f^{(n)}(\xi) = p_n^{(n)}(\xi)

Part 3.

g = f - p_n

has the

n + 1

distinct zeros

x_0, \dots, x_n

. By Rolle's theorem applied

n

times,

g^{(n)}

has a zero

\xi

in the interval.

f[x_0, \dots, x_n] = f^{(n)}(\xi)/n!

p_n

has leading term

f[x_0, \dots, x_n]\,x^n

, and the

n

-th derivative of that is

n!\,f[x_0, \dots, x_n]

, a constant.

Part 3 explains the name: a divided difference of order $n$ is a derivative of order $n$ , divided by $n!$ , evaluated at an unknown point. The first one, $f[x_0, x_1] = \dfrac{f(x_1) - f(x_0)}{x_1 - x_0}$ , is the slope of the chord, and part 3 is the mean value theorem.

Worked example

Build the Newton form of the interpolating polynomial through (1, 2), (2, 3), (4, 9) with a divided-difference table. Then add the point (3, 5) and update.

Solution.

Order-0 column. $f[1] = 2$ , $f[2] = 3$ , $f[4] = 9$ . The data values.
Order-1 column. $f[1, 2] = \dfrac{3 - 2}{2 - 1} = 1$ and $f[2, 4] = \dfrac{9 - 3}{4 - 2} = 3$ . The quotient recursion: slope of each chord.
Order-2 entry. $f[1, 2, 4] = \dfrac{3 - 1}{4 - 1} = \tfrac23$ . Same recursion, applied to the two slopes.
Read off the top edge. $p_2(x) = 2 + 1\cdot(x - 1) + \tfrac23(x - 1)(x - 2)$ . The coefficients are the top diagonal $f[x_0], f[x_0, x_1], f[x_0, x_1, x_2]$ .
Expand to check. $2 + x - 1 + \tfrac23(x^2 - 3x + 2) = \tfrac23 x^2 - x + \tfrac73$ . The same polynomial as the Lagrange example above, as uniqueness demands.
Add the node $(3, 5)$ as $x_3$ . New entries: $f[4, 3] = \dfrac{5 - 9}{3 - 4} = 4$ , $f[2, 4, 3] = \dfrac{4 - 3}{3 - 2} = 1$ , $f[1, 2, 4, 3] = \dfrac{1 - \frac23}{3 - 1} = \tfrac16$ . Only a new bottom row of the table is computed; the old entries stay.
Update the polynomial. $p_3(x) = p_2(x) + \tfrac16(x - 1)(x - 2)(x - 4)$ . Adding a node adds one term. Check at $x = 3$ : $p_2(3) = \tfrac{16}3$ and $\tfrac16(2)(1)(-1) = -\tfrac13$ , giving $5$ .

Answer. $p_2(x) = 2 + (x - 1) + \tfrac23(x - 1)(x - 2)$ , and $p_3 = p_2 + \tfrac16(x-1)(x-2)(x-4)$ . The Lagrange form would have needed all four basis polynomials recomputed.

Your answerFor the data (1, 2), (2, 3), (4, 9), what is the divided difference f[1, 2, 4]?

Numerical integration (Apostol 15.19 and 15.20)

Polynomials are easy to integrate, so a natural way to approximate $\int_a^b f$ is to integrate an interpolant.

Definition A

Quadrature rule

A quadrature rule approximates $\int_a^b f(x)\,dx$ by $\sum_{j=0}^{n} w_j f(x_j)$ , with nodes $x_j$ and fixed weights $w_j$ . If $p_n = \sum_j f(x_j)\ell_j$ is the interpolant, then $\int_a^b p_n = \sum_j f(x_j)\int_a^b \ell_j$ , so the weights are $w_j = \int_a^b \ell_j(x)\,dx$ : the rule integrates the interpolant exactly.

In one lineAn approximation $\int_a^b f \approx \sum_j w_j f(x_j)$ by a weighted sum of values at chosen nodes. Integrating the interpolating polynomial gives one.

Why it existsEvaluating f at a handful of points is all a computer can do with a function it has no formula for, so every numerical integral is of this form.

Definition A

Trapezoidal rule

With the two nodes $a$ and $b$ (so $n = 1$ ), the interpolant is the chord, and its integral is the trapezoidal rule $T = \tfrac{b-a}2\big(f(a) + f(b)\big)$ .

In one line $T = \dfrac{b-a}{2}\big(f(a) + f(b)\big)$ : the area of the trapezoid under the chord from $(a, f(a))$ to $(b, f(b))$ , which is the integral of the linear interpolant.

Why it existsIt is the simplest quadrature rule, and it is exact for straight lines.

Definition A

Simpson's rule

With the three nodes $a$ , $m = \frac{a+b}{2}$ and $b$ (so $n = 2$ ), the integral of the quadratic interpolant is Simpson’s rule $S = \tfrac{b-a}{6}\big(f(a) + 4f(m) + f(b)\big)$ .

In one line $S = \dfrac{b-a}{6}\big(f(a) + 4f(m) + f(b)\big)$ with $m = \frac{a+b}2$ : the integral of the quadratic interpolant at $a$ , $m$ , $b$ .

Why it existsFitting a parabola instead of a line gains two orders of accuracy, and by a happy accident the rule is even exact for cubics.

Where does the $\frac16, \frac46, \frac16$ come from? Shift the interval to $[-h, h]$ (so $b - a = 2h$ and $m = 0$ ). The Lagrange basis polynomials for the nodes $-h, 0, h$ are $\ell_{-1}(t) = \dfrac{t(t - h)}{2h^2}$ , $\ell_0(t) = \dfrac{h^2 - t^2}{h^2}$ and $\ell_1(t) = \dfrac{t(t + h)}{2h^2}$ . Integrating over $[-h, h]$ (odd powers of $t$ integrate to $0$ ): $\int \ell_{-1} = \dfrac{1}{2h^2}\cdot\dfrac{2h^3}{3} = \dfrac h3$ , $\int \ell_0 = \dfrac{2h^3 - \frac{2h^3}3}{h^2} = \dfrac{4h}3$ , $\int \ell_1 = \dfrac h3$ . So $S = \dfrac h3\big(f(-h) + 4f(0) + f(h)\big)$ , which is the formula above with $b - a = 2h$ .

Lemma A standard

Cauchy's mean value theorem

Let $g, k$ be continuous on $[0, h]$ and differentiable on $(0, h)$ , with $k'(t) \ne 0$ for $0 < t < h$ . Then $k(h) \ne k(0)$ and there is $\tau \in (0, h)$ with $\big(g(h) - g(0)\big)\,k'(\tau) = \big(k(h) - k(0)\big)\,g'(\tau)$ .

In one lineIf $g, k$ are continuous on $[0, h]$ , differentiable on $(0, h)$ and $k' \ne 0$ there, then $\dfrac{g(h) - g(0)}{k(h) - k(0)} = \dfrac{g'(\tau)}{k'(\tau)}$ for some $\tau \in (0, h)$ .

Why it is trueIt compares how two functions change over the same interval. We use it to read the error of a rule off the derivative of its error.

Proofstandard

The idea: Apply Rolle's theorem to the combination of g and k that has equal end values.

k(h) \ne k(0)

k(h) = k(0)

, Rolle's theorem would give a zero of

k'

(0, h)

, against the hypothesis.

Let

\varphi(t) = \big(g(h) - g(0)\big)k(t) - \big(k(h) - k(0)\big)g(t)

A linear combination of

g

and

k

. Both end values come out equal:

\varphi(0) = g(h)k(0) - g(0)k(h) = \varphi(h)

There is

\tau \in (0, h)

with

\varphi'(\tau) = 0

, which is the claimed equation.

Rolle's theorem for

\varphi

, which is continuous on

[0, h]

and differentiable inside, with

\varphi(0) = \varphi(h)

Theorem A Apostol 15.19 and 15.20; Burden and Faires 4.3

Error of the trapezoidal and Simpson rules

Let $f \in C^2[a, b]$ for the first formula and $f \in C^4[a, b]$ for the second. Then there are $\eta, \eta' \in (a, b)$ with

\int_a^b f - T = -\frac{(b-a)^3}{12}\,f''(\eta), \qquad \int_a^b f - S = -\frac{(b-a)^5}{2880}\,f^{(4)}(\eta').

In one lineOn an interval of width $b - a$ , $\int_a^b f - T = -\dfrac{(b-a)^3}{12}f''(\eta)$ and $\int_a^b f - S = -\dfrac{(b-a)^5}{2880}f^{(4)}(\eta)$ for some $\eta \in (a, b)$ (a different $\eta$ in each).

Why it is trueEach error is the next derivative times a power of the width. The trapezoidal rule is exact for lines (f'' = 0), Simpson's for cubics (f'''' = 0).

Proofstandard

The idea: Centre the interval at 0, and let the half-width t grow from 0. The error G(t) of the rule on [-t, t] starts at 0 and has derivatives that are easy: G' for the trapezoid, G', G'' and G''' for Simpson. Cauchy's mean value theorem then expresses G(h) through the last derivative.

Set

a = -h

b = h

Translate so that the interval is

[-h, h]

with

h = (b-a)/2

. Integrals and the rules are unchanged by translation.

Define the error as a function of the half-width.

Trapezoid: define

G(t) = \int_{-t}^{t} f - t\big(f(-t) + f(t)\big)

for

0 \le t \le h

, the error of the rule on

[-t, t]

. Then

G(0) = 0

and

G(h) = \int f - T

G'(t) = -t\big(f'(t) - f'(-t)\big)

Differentiate:

\frac{d}{dt}\int_{-t}^{t} f = f(t) + f(-t)

(fundamental theorem of calculus, 1.5), and the product rule on

t(f(-t) + f(t))

. The terms

f(t) + f(-t)

cancel.

G'(t) = -2t^2 f''(\eta_t)

Mean value theorem for

f'

[-t, t]

f'(t) - f'(-t) = 2t\,f''(\eta_t)

for some

\eta_t \in (-t, t)

\dfrac{G(h)}{h^3} = \dfrac{G'(\tau)}{3\tau^2} = -\dfrac{2}{3}f''(\eta)

with

\eta = \eta_\tau \in (-h, h)

Cauchy's mean value theorem with

g = G

and

k(t) = t^3

(so

k' = 3t^2 \ne 0

(0, h)

), using

G(0) = 0 = k(0)

This is the trapezoidal formula.

G(h) = -\frac23 h^3 f''(\eta)

, and

\frac23h^3 = \frac{(2h)^3}{12} = \frac{(b-a)^3}{12}

Same plan for Simpson’s rule.

Simpson: define

G(t) = \int_{-t}^{t} f - \frac t3\big(f(-t) + 4f(0) + f(t)\big)

. Again

G(0) = 0

and

G(h) = \int f - S

G(0) = G'(0) = G''(0) = 0

Differentiate three times as before; each time terms cancel.

G'(t) = \frac23\big(f(t) + f(-t)\big) - \frac43 f(0) - \frac t3\big(f'(t) - f'(-t)\big)

, so

G'(0) = 0

. Then

G''(t) = \frac13\big(f'(t) - f'(-t)\big) - \frac t3\big(f''(t) + f''(-t)\big)

, so

G''(0) = 0

G'''(t) = -\frac{2t^2}{3}f^{(4)}(\eta_t)

G'''(t) = -\frac t3\big(f'''(t) - f'''(-t)\big)

(the other terms cancel). The mean value theorem for

f'''

gives

f'''(t) - f'''(-t) = 2t\,f^{(4)}(\eta_t)

\dfrac{G(h)}{h^5} = \dfrac{G'''(t_3)}{60\,t_3^2} = -\dfrac{2}{3\cdot 60}f^{(4)}(\eta') = -\dfrac{1}{90}f^{(4)}(\eta')

Cauchy's mean value theorem three times, with

k(t) = t^5

, then

5t^4

, then

20t^3

, each with the previous function's value

0

0

\frac{G(h)}{h^5} = \frac{G'(t_1)}{5t_1^4} = \frac{G''(t_2)}{20t_2^3} = \frac{G'''(t_3)}{60t_3^2}

with

0 < t_3 < t_2 < t_1 < h

This is Simpson’s formula.

\frac{h^5}{90} = \frac{(2h)^5}{32\cdot 90} = \frac{(b-a)^5}{2880}

Why· The picture behind the proof

For the trapezoid, draw the chord under a convex curve (say $x^2$ ): the curve bulges below the chord, and the trapezoid overestimates the area. The sign in the formula says this: $f'' > 0$ gives $\int f - T < 0$ . The size of the bulge is controlled by how fast the slope changes, which is $f''$ , and by $(b-a)^2$ for the height times $(b - a)$ for the width. That is $(b-a)^3 f''$ up to a constant, and $\frac1{12}$ is the exact constant. For Simpson, a parabola through three points copies a cubic perfectly (the extra symmetry cancels the cubic term), so the first derivative that can create an error is the fourth.

A single interval is too wide for accuracy. Cut $[a, b]$ into $N$ equal pieces of width $h = (b-a)/N$ and apply the rule to each. Adding up the errors gives the composite rules. In signed form, for continuous $f''$ , the sum of the $N$ values $f''(\eta_i)$ divided by $N$ is itself a value $f''(\eta)$ by the intermediate value theorem, so $\int_a^b f - T_N = -\dfrac{(b-a)h^2}{12}f''(\eta)$ . Simpson’s rule gives in the same way $\int_a^b f - S_N = -\dfrac{(b-a)h^4}{180}f^{(4)}(\eta')$ . The theorem below states the absolute bounds, which are what we use.

Theorem A Apostol 15.19 and 15.20

Composite trapezoidal and Simpson rules

Let $f_i = f(a + ih)$ and $M_k = \max_{[a,b]}|f^{(k)}|$ . For $f \in C^2$ : $\big|\int_a^b f - T_N\big| \le \dfrac{(b-a)\,h^2}{12}\,M_2$ . For $f \in C^4$ and $N$ even: $\big|\int_a^b f - S_N\big| \le \dfrac{(b-a)\,h^4}{180}\,M_4$ .

In one lineWith $N$ pieces of width $h$ : $T_N = h\big(\tfrac12 f_0 + f_1 + \cdots + f_{N-1} + \tfrac12 f_N\big)$ has error at most $\dfrac{(b-a)h^2}{12}M_2$ , and (for $N$ even) $S_N = \dfrac h3\big(f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 4f_{N-1} + f_N\big)$ has error at most $\dfrac{(b-a)h^4}{180}M_4$ .

Why it is trueHalving h cuts the trapezoid error by 4 and the Simpson error by 16. That is the 'order of the method', and it is what you see in the error table of the PyLab below.

Proofstandard

The idea: Apply the one-interval error to each piece and add. The piece count is N or N/2, and that turns the cubic and fifth power of the width into the square and fourth.

Trapezoid: total error at most

\dfrac{Nh^3}{12}M_2

On each of the

N

pieces (width

h

) the trapezoid error is

\frac{h^3}{12}|f''(\eta_i)| \le \frac{h^3}{12}M_2

. Adding

N

of them:

N\frac{h^3}{12}M_2

= \dfrac{(b-a)h^2}{12}M_2

Nh = b - a

Simpson: total error at most

\dfrac N2\cdot\dfrac{h^5}{90}M_4

Simpson's rule uses pairs of pieces: each pair has width

2h

(so half-width

h

), and the error on it is

\frac{h^5}{90}|f^{(4)}(\eta_i)|

by the theorem with half-width

h

. There are

N/2

pairs.

= \dfrac{(b-a)h^4}{180}M_4

Nh = b - a

, so

\frac N2\cdot\frac{h^5}{90} = \frac{(b-a)h^4}{180}

Worked example

Estimate ∫₀¹ eˣ dx = e − 1 ≈ 1.718282 with the trapezoidal and Simpson rules using 4 pieces. Compare the error with the bound.

Solution.

Set up. $h = \tfrac14$ , nodes $0, \tfrac14, \tfrac12, \tfrac34, 1$ . Four pieces of width $\frac14$ .
Trapezoid. $T_4 = \tfrac14\big(\tfrac12e^0 + e^{1/4} + e^{1/2} + e^{3/4} + \tfrac12e^1\big) \approx 1.727222$ . The weights are $\frac12, 1, 1, 1, \frac12$ times $h$ .
Simpson. $S_4 = \tfrac1{12}\big(e^0 + 4e^{1/4} + 2e^{1/2} + 4e^{3/4} + e^1\big) \approx 1.718319$ . The weights are $1, 4, 2, 4, 1$ times $h/3$ .
Actual errors. $T_4 - (e - 1) \approx 0.008940$ and $S_4 - (e - 1) \approx 0.000037$ . Subtract from $1.718282$ . The trapezoid overestimates because $e^x$ is convex.
Bounds. $M_2 = M_4 = e$ on $[0, 1]$ . Trapezoid: $\dfrac{1 \cdot (1/16)}{12}e \approx 0.01416$ . Simpson: $\dfrac{1\cdot(1/256)}{180}e \approx 0.0000590$ . The composite bounds with $b - a = 1$ , $h = \frac14$ .

Answer. $T_4 \approx 1.7272$ (error $0.0089$ , bound $0.0142$ ) and $S_4 \approx 1.71832$ (error $0.000037$ , bound $0.000059$ ). With the same four pieces, Simpson is about $240$ times more accurate.

Worked example

How many pieces does each rule need to compute ∫₀¹ e^(−x²) dx to within 10⁻⁶?

Solution.

Derivatives. $f'' = (4x^2 - 2)e^{-x^2}$ and $f'''' = (16x^4 - 48x^2 + 12)e^{-x^2}$ . On $[0, 1]$ the largest values of $|f''|$ and $|f''''|$ are $M_2 = 2$ and $M_4 = 12$ , both at $x = 0$ . Differentiate twice, then twice more. Check the values at $0$ directly: $-2$ and $12$ .
Trapezoid. Need $\dfrac{1\cdot h^2}{12}\cdot 2 \le 10^{-6}$ , so $h^2 \le 6 \times 10^{-6}$ and $N \ge \sqrt{1/(6\times10^{-6})} \approx 408.2$ , so $N = 409$ . $h = 1/N$ , and $N^2 \ge 166667$ .
Simpson. Need $\dfrac{h^4}{180}\cdot 12 \le 10^{-6}$ , so $h^4 \le 1.5 \times 10^{-5}$ and $N^4 \ge 66667$ , giving $N \ge 16.07$ . $N$ must be even, so $N = 18$ . $12/180 = 1/15$ .
Truth. The actual errors are $3.7 \times 10^{-7}$ for the trapezoid with $N = 409$ and $7.8 \times 10^{-8}$ for Simpson with $N = 18$ (already $1.2 \times 10^{-7}$ at $N = 16$ ). The bounds use the worst case of the derivative everywhere and so overestimate.

Answer. The bounds guarantee the trapezoidal rule at $N = 409$ and Simpson’s rule at $N = 18$ . Simpson needs about $23$ times fewer pieces. The integral is $0.746824\ldots$ , the area behind the normal curve.

Your answerSimpson's rule with one piece (b − a = 2, nodes 0, 1, 2) for the integral of x⁴ from 0 to 2. What value does it give? (The exact value is 6.4.)

Your answerHow many pieces N does the composite trapezoidal rule need to guarantee an error below 10⁻⁴ for the integral of eˣ over [0, 1]? (M₂ = e.)

Quick checkYou double the number of pieces N in the composite Simpson rule for a smooth function. By roughly what factor does the error shrink?

Python labTrapezoid versus Simpson: the error as the number of pieces doubles

import numpy as np

f, a, b, exact = np.exp, 0.0, 1.0, np.e - 1

def trapezoid(f, a, b, n):
    x = np.linspace(a, b, n + 1); h = (b - a) / n
    return h * (f(x[0]) / 2 + f(x[1:-1]).sum() + f(x[-1]) / 2)

def simpson(f, a, b, n):                          # n must be even
    x = np.linspace(a, b, n + 1); h = (b - a) / n
    return h / 3 * (f(x[0]) + 4 * f(x[1:-1:2]).sum() + 2 * f(x[2:-1:2]).sum() + f(x[-1]))

print(" n    trap error    ratio   Simpson error   ratio")
old_t = old_s = None
for n in (2, 4, 8, 16, 32, 64):
    et, es = abs(trapezoid(f, a, b, n) - exact), abs(simpson(f, a, b, n) - exact)
    rt = f"{old_t / et:6.2f}" if old_t else "      "
    rs = f"{old_s / es:6.2f}" if old_s else "      "
    print(f"{n:<4d} {et:<13.3e} {rt}  {es:<14.3e} {rs}")
    old_t, old_s = et, es

Try this. Look at the ratio columns: each time n doubles the trapezoid error is divided by 4 (ratios 3.99, 4.00, ...) and the Simpson error by 16 (15.65, 15.91, 15.98, ...), as the error bounds predict. Now change f to lambda x: x**3, with exact = 0.25: Simpson's error is 0 up to rounding for every n, because the fourth derivative is 0, while the trapezoid error is exactly 1/(4 n^2).

Beyond the course: for understanding· Why the ratio is exactly 4 (and when it is not)

The Euler summation (Euler–Maclaurin) formula, which Apostol proves in 15.22, refines the trapezoid error to $T_N - I = \dfrac{h^2}{12}\big(f'(b) - f'(a)\big) - \dfrac{h^4}{720}\big(f'''(b) - f'''(a)\big) + \cdots$ . The leading term is proportional to $h^2$ , so doubling $N$ divides the error by $4$ , as long as $f'(b) \ne f'(a)$ . When the slopes at both ends agree, as for a smooth periodic function integrated over a full period, all the terms vanish to every order and the trapezoidal rule converges faster than any power of $h$ . It is the best rule for that case, and the reason the discrete Fourier transform is accurate. The same formula also gives the closed forms of sums such as $1^2 + 2^2 + \cdots + N^2$ . Apostol also covers the forward difference operator and factorial polynomials (15.11, 15.12): the Newton formula for equally spaced nodes, where every divided difference is a difference divided by a power of $h$ and a factorial.

Rounding: how a computer stores numbers

Everything above was exact arithmetic. A computer stores a number with a fixed number of binary digits, so it works with a finite set of numbers and rounds every result to the nearest one.

Definition B

Floating-point number

A floating-point number with $p$ significant bits is $x = \pm m \cdot 2^{e}$ , where $m$ is an integer with $0 \le m < 2^p$ and $e$ lies between a smallest and a largest value. Numbers of size about $2^k$ are spaced $2^{k - p + 1}$ apart: the gaps grow with the size, but the gap relative to the number stays between $2^{-p}$ and $2^{1-p}$ . Rounding any real number $x$ in range to the nearest one gives $\mathrm{fl}(x) = x(1 + \delta)$ with $|\delta| \le u = 2^{-p}$ . This $u$ is the machine epsilon of 3.5: the largest relative rounding error. The standard model for arithmetic is that each operation is the exact result rounded: $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ , $|\delta| \le u$ , for $\circ \in \{+, -, \times, \div\}$ .

In one lineA number of the form $\pm m \cdot 2^{e}$ with an integer $m$ having at most $p$ binary digits (the significand) and an exponent $e$ in a fixed range. Every real number inside the range is rounded to the nearest one.

Why it existsIt gives constant relative precision across a huge range of sizes, which is what scientific computing needs, at the price that every operation can be off by a small relative amount.

Format	Bits (sign, exponent, stored fraction)	Significant bits $p$	$u = 2^{-p}$	Largest number
float64 (double)	1, 11, 52	53	$1.1 \times 10^{-16}$	$1.8 \times 10^{308}$
float32 (single)	1, 8, 23	24	$6.0 \times 10^{-8}$	$3.4 \times 10^{38}$
float16 (half)	1, 5, 10	11	$4.9 \times 10^{-4}$	$65504$
bfloat16	1, 8, 7	8	$3.9 \times 10^{-3}$	$3.4 \times 10^{38}$

(The leading binary digit of a non-zero significand is always $1$ and is not stored, so $p$ is one more than the stored fraction bits.) The examples below use double precision. Two facts follow at once. A number such as $0.1$ has no finite binary expansion, so even storing it rounds: in float64, $0.1 + 0.2$ equals $0.30000000000000004$ , not $0.3$ . And once the numbers are large enough the gap between neighbours exceeds $1$ : in float32, $2^{24} = 16777216$ is the last point where every integer is representable, so $16777216 + 1$ rounds back to $16777216$ .

Definition B

Catastrophic cancellation

Suppose $\tilde x = x(1 + \delta_1)$ and $\tilde y = y(1 + \delta_2)$ with $|\delta_i| \le u$ . The computed difference has relative error

\frac{|(\tilde x - \tilde y) - (x - y)|}{|x - y|} = \frac{|x\delta_1 - y\delta_2|}{|x - y|} \le u\,\frac{|x| + |y|}{|x - y|}.

When $x$ and $y$ are nearly equal, the right side is far larger than $u$ . This is catastrophic cancellation. The factor $(|x| + |y|)/|x - y|$ is the condition number of the operation $(x, y) \mapsto x - y$ in the sense of 3.5: subtraction of nearly equal numbers is ill-conditioned, however well the subtraction is performed.

In one lineThe loss of most or all correct digits when two nearly equal numbers, each carrying a small rounding error, are subtracted.

Why it existsSubtraction itself is computed exactly, but it exposes errors that earlier steps already made, and the result is small so the same absolute error is a large relative error.

Worked example

Compute (1 − cos x)/x² for x = 10⁻⁸ in double precision. Why is the answer wrong, and how do you fix it?

Solution.

The true value. For small $x$ , $\cos x = 1 - \tfrac{x^2}2 + \cdots$ , so $(1 - \cos x)/x^2 \to \tfrac12$ . Taylor polynomial of cosine (1.4).
What the computer sees. $\cos(10^{-8}) = 1 - 5 \times 10^{-17} + \cdots$ . The doubles immediately below $1$ are spaced $1.1 \times 10^{-16}$ apart, and $5 \times 10^{-17}$ is less than half of that gap ( $5.5 \times 10^{-17}$ ), so the nearest double is $1$ itself. The rounding of the cosine has erased the entire difference from $1$ .
The subtraction. $1 - \cos x$ evaluates to exactly $0$ , and the whole expression to $0$ . The subtraction is exact, but it subtracts a number from the rounded copy of itself.
The cancellation factor. Here $(|1| + |\cos x|)/|1 - \cos x| \approx 2/(5\times10^{-17}) = 4 \times 10^{16}$ , and $u \cdot 4 \times 10^{16} \approx 4$ : the bound on the relative error is above $1$ . The formula for cancellation predicts that no digit is guaranteed.
Fix by algebra. Use $1 - \cos x = 2\sin^2(x/2)$ , so the expression is $\dfrac{2\sin^2(x/2)}{x^2}$ . The half-angle identity removes the subtraction. $\sin(x/2)$ is computed to full relative accuracy.
Evaluate. $\dfrac{2\sin^2(5\times10^{-9})}{10^{-16}} = 0.5$ to about $16$ digits. Checked in the PyLab below.

Answer. The naive formula returns $0.0$ ; the rewritten one returns $0.5$ . The two formulas are equal as mathematics and different as algorithms. This is why library functions such as expm1 and log1p exist: they compute $e^x - 1$ and $\ln(1 + x)$ for small $x$ without cancellation.

Worked example

You estimate the derivative of sin at x = 1 by (sin(1 + h) − sin 1)/h. As h shrinks the error first falls and then rises. Find the best h.

Solution.

Truncation error. By Taylor (1.4), $\dfrac{f(1+h) - f(1)}{h} = f'(1) + \tfrac h2 f''(\xi)$ . So the truncation error is at most $\tfrac h2 M_2$ . It shrinks in proportion to $h$ .
Rounding error. Each computed value of $\sin$ carries an error up to about $u|f|$ , so the difference carries up to about $2u|f|$ , and dividing by $h$ gives about $2uM_0/h$ . It grows as $h$ shrinks: this is cancellation, since $\sin(1+h) \approx \sin 1$ .
Add and minimize. $E(h) \approx \tfrac h2 M_2 + \dfrac{2uM_0}{h}$ is smallest where its derivative is zero: $\tfrac{M_2}{2} = \dfrac{2uM_0}{h^2}$ , so $h = 2\sqrt{uM_0/M_2}$ . With $M_0 \approx M_2 \approx 1$ and $u \approx 1.1 \times 10^{-16}$ : $h \approx 2 \times 10^{-8}$ . The rule of thumb: the best step is about $\sqrt{u} \approx 10^{-8}$ .
Compare with the computer. The measured errors are $4.2 \times 10^{-3}$ at $h = 10^{-2}$ , $4.2 \times 10^{-5}$ at $10^{-4}$ , $3.0 \times 10^{-9}$ at $10^{-8}$ , $5.8 \times 10^{-8}$ at $10^{-10}$ , $4.3 \times 10^{-5}$ at $10^{-12}$ and $3.7 \times 10^{-3}$ at $10^{-14}$ . The first two fall as $h$ ; the minimum is near $h = 10^{-8}$ ; after that rounding dominates.

Answer. The best step is near $10^{-8}$ , giving about $8$ correct digits: half of the $16$ available are lost to cancellation. For the central difference $\frac{f(x+h) - f(x-h)}{2h}$ the truncation error is proportional to $h^2$ , the best step is near $u^{1/3} \approx 5\times10^{-6}$ , and the best error is near $10^{-11}$ . This trade-off is why automatic differentiation (3.8) is used for gradients in machine learning instead of finite differences.

Quick checkIn double precision, 1 + 10⁻¹⁷ evaluates to exactly 1. What does that say about u?

Your answerx = 1.0000001 and y = 1 are stored with relative error at most u = 1.1×10⁻¹⁶ each. The bound on the relative error of x − y is u(|x| + |y|)/|x − y|. What is it, measured in units of 10⁻⁹? (For example, 3 means 3×10⁻⁹.)

Newton’s method, again (link to 4.2)

In 4.2 Newton’s method $x_{k+1} = x_k - f(x_k)/f'(x_k)$ was shown to be a contraction near a simple root, with the derivative of the iteration map equal to $0$ at the root. The Taylor remainder of 1.4 gives the sharper statement.

Theorem A standard; Burden and Faires 2.3

Newton's method converges quadratically

Let $f$ be twice continuously differentiable near a root $r$ , with $f'(x_k) \ne 0$ . Then the Newton step $x_{k+1} = x_k - f(x_k)/f'(x_k)$ satisfies, for some $\xi_k$ between $x_k$ and $r$ ,

x_{k+1} - r = \frac{f''(\xi_k)}{2f'(x_k)}\,(x_k - r)^2.

In one lineIf $f(r) = 0$ , $f'(x_k) \ne 0$ and $f$ is twice differentiable, then $x_{k+1} - r = \dfrac{f''(\xi_k)}{2f'(x_k)}\,(x_k - r)^2$ for some $\xi_k$ between $x_k$ and $r$ .

Why it is trueThe new error is the square of the old error times a constant. Squaring an error below one doubles the number of correct digits at every step, which is why Newton's method is so fast.

Proofstandard

The idea: Expand f around x_k with Taylor's theorem and evaluate at the root.

0 = f(r) = f(x_k) + f'(x_k)(r - x_k) + \tfrac12 f''(\xi_k)(r - x_k)^2

Taylor's theorem with the Lagrange remainder (1.4), centre

x_k

, order

1

, evaluated at

r

0 = \dfrac{f(x_k)}{f'(x_k)} + (r - x_k) + \dfrac{f''(\xi_k)}{2f'(x_k)}(r - x_k)^2

Divide by

f'(x_k) \ne 0

Solve for the Newton correction.

Rearrange:

\dfrac{f(x_k)}{f'(x_k)} = (x_k - r) - C\,(x_k - r)^2

with

C = \dfrac{f''(\xi_k)}{2f'(x_k)}

x_{k+1} - r = C\,(x_k - r)^2

x_{k+1} - r = (x_k - r) - \dfrac{f(x_k)}{f'(x_k)}

; substitute the line above.

The identity holds at each step, but it does not by itself say that the iteration converges. For that we need the root to be a simple one, $f'(r) \ne 0$ , and $f''$ continuous. Then $f'$ is non-zero near $r$ , so $C = f''(\xi_k)/(2f'(x_k))$ is bounded by some number $K$ on an interval $[r - \delta, r + \delta]$ . Write $e_k = |x_k - r|$ . If $x_0$ lies in that interval and $e_0 < 1/K$ , then $e_{k+1} \le K e_k^2 < e_k$ , so every later iterate stays in the interval, and the errors satisfy $e_k \le (Ke_0)^{2^k}/K$ : the exponent doubles at every step. This is quadratic convergence. At a double root ( $f'(r) = 0$ ) the argument fails and Newton’s method only halves the error each step.

Worked example compare the NumCheck in 4.2

Apply Newton's method to f(x) = x² − 2 from x₀ = 3/2. Check the error formula and watch the digits.

Solution.

The step. $x_{k+1} = x_k - \dfrac{x_k^2 - 2}{2x_k} = \dfrac{x_k}{2} + \dfrac{1}{x_k}$ . $f' = 2x$ .
The error formula. Here $f'' = 2$ is constant, so $\xi_k$ does not matter: $x_{k+1} - \sqrt2 = \dfrac{2}{2\cdot 2x_k}(x_k - \sqrt2)^2 = \dfrac{(x_k - \sqrt2)^2}{2x_k}$ . Direct from the theorem. Direct algebra confirms it: $\frac{x}{2} + \frac1x - \sqrt2 = \frac{x^2 - 2\sqrt2 x + 2}{2x}$ .
Iterate. $x_0 = 1.5$ , error $8.6\times10^{-2}$ . $x_1 = \tfrac{17}{12} = 1.41667$ , error $2.5 \times 10^{-3}$ . $x_2 = 1.414215686$ , error $2.1\times10^{-6}$ . $x_3 = 1.41421356237469$ , error $1.6\times10^{-12}$ . Check one step: $e_1 = (0.0858)^2/(2\cdot1.5) = 0.00245$ .
Read the digits. Correct digits: $1, 3, 6, 12$ . They double each time, as $e_{k+1} \approx e_k^2/(2\sqrt2)$ predicts.
The floor. In double precision $x_4 = 1.414213562373095$ is within $2.2 \times 10^{-16}$ of $\sqrt2$ (one gap between neighbouring doubles), and $x_5$ is the same number. The computed value of $f(x_k)$ is only accurate to about $u$ in relative terms, and so the iteration stops improving at the machine epsilon of 3.5.

Answer. The error formula is exact here, and six digits become twelve become the full $16$ . Newton’s method converges to roots as fast as a number format allows, so rounding, not the method, is what stops it.

Your answerNewton's method for x² − 2 has e_(k+1) = e_k²/(2 x_k). If x_k is close to √2 and e_k = 10⁻³, about how large is e_(k+1), measured in units of 10⁻⁷? (For example, 3 means 3×10⁻⁷.)

Quick checkNewton's method has e_(k+1) ≈ C e_k². After a step the error is 10⁻⁴. Roughly what is it after the next step (C near 1)?

Examples and non-examples

More nodes, worse fit. Interpolating $1/(1 + 25x^2)$ at equally spaced nodes: the largest error on $[-1, 1]$ is $0.62$ at $n = 6$ and $29$ at $n = 18$ . Non-example of the rule “more points means better”.

Smooth is not enough, nodes matter. The same function interpolated at Chebyshev nodes has error $0.26$ at $n = 6$ and $0.0225$ at $n = 18$ . The nodes, not the function, were the problem.

A function that equal spacing handles. For $1/(1 + x^2)$ (poles at $\pm i$ , farther away) equally spaced interpolation converges: error $0.0068$ at $n = 6$ , about $0.0001$ at $n = 14$ and $0.000015$ at $n = 18$ .

Exact rules. The trapezoidal rule is exact for lines, Simpson’s for polynomials up to degree three: for $\int_0^2 x^3\,dx = 4$ Simpson with the single piece gives $\frac{2}{6}(0 + 4 + 8) = 4$ exactly. For $x^4$ it does not (the Simpson NumCheck above gives $\frac{20}{3}$ for the true value $\frac{32}5$ ).

Same function, two algorithms. $(1 - \cos x)/x^2$ and $2\sin^2(x/2)/x^2$ are equal for every $x \ne 0$ and give $0$ and $0.5$ at $x = 10^{-8}$ .

Where it is used

Used inside maths· Inside maths: Taylor's theorem is interpolation at one repeated node

Let all $n + 1$ nodes slide together onto a single point $a$ . The product $w(x) = (x - x_0)\cdots(x - x_n)$ becomes $(x - a)^{n+1}$ , and the divided difference $f[x_0, \dots, x_n] = f^{(n)}(\xi)/n!$ becomes $f^{(n)}(a)/n!$ . The Newton form of the interpolant becomes $\sum_k \frac{f^{(k)}(a)}{k!}(x - a)^k$ , the Taylor polynomial, and the interpolation error formula becomes $\frac{f^{(n+1)}(\xi)}{(n+1)!}(x - a)^{n+1}$ , exactly the Lagrange remainder of 1.4. The proof by Rolle’s theorem is the same proof with the zeros all piled on top of each other. So Taylor polynomials are the extreme case of interpolation, with all information gathered at one point, and Chebyshev interpolation is the opposite extreme, with the information spread as evenly as possible.

Used inside maths· Inside maths: Chebyshev polynomials are orthogonal, and the Gram–Schmidt process produces them

Put $x = \cos\theta$ in $\int_{-1}^{1} \dfrac{T_m(x)T_n(x)}{\sqrt{1 - x^2}}\,dx$ . Then $dx = -\sin\theta\,d\theta$ , $\sqrt{1 - x^2} = \sin\theta$ , and $T_m T_n = \cos m\theta\cos n\theta$ , so the integral becomes $\int_0^\pi \cos m\theta\cos n\theta\,d\theta$ . By the product formula $\cos m\theta\cos n\theta = \tfrac12\big(\cos(m-n)\theta + \cos(m+n)\theta\big)$ , this is $0$ for $m \ne n$ , $\pi/2$ for $m = n \ge 1$ , and $\pi$ for $m = n = 0$ . (Numerical check: the integral for $(m, n) = (2, 3)$ is $0$ up to rounding, for $(3, 3)$ it is $1.5708$ , and for $(0, 0)$ it is $3.1416$ .) So with the inner product $\langle f, g\rangle = \int_{-1}^{1} f g/\sqrt{1-x^2}\,dx$ of 2.9, the polynomials $T_0, T_1, T_2, \dots$ are orthogonal. They are what Gram–Schmidt gives when applied to $1, x, x^2, \dots$ (up to scaling). Least squares approximation by polynomials in this inner product (2.10) is nearly the best uniform approximation, which is why Chebyshev expansions are used in practice.

Used in the real world· Why machine learning libraries use float32, float16 and bfloat16

A neural network has billions of numbers. Halving the bits per number halves the memory and, on a GPU or TPU, roughly doubles the speed of the matrix multiplications of 2.5. The cost is the unit roundoff in the table above: $u = 6\times10^{-8}$ for float32, $4.9\times10^{-4}$ for float16 and $3.9\times10^{-3}$ for bfloat16. Three real consequences.

Sums stall. In float32, once a running total reaches $2^{24} = 16777216$ , adding $1$ does nothing (shown in the PyLab). Adding a long list of small numbers into a float16 accumulator stalls far sooner, at $2^{11} = 2048$ . So libraries accumulate dot products in float32 even when the inputs are 16-bit. Micikevicius and colleagues’ paper Mixed Precision Training (ICLR 2018) trains with float16 values and keeps a float32 master copy of the weights and float32 accumulations, so that small updates are not rounded away.
Range. Float16 can only represent numbers between about $6\times10^{-8}$ and $65504$ . Small gradients underflow to $0$ , so the same paper multiplies the loss by a constant (loss scaling) before back-propagation and divides afterwards.
bfloat16. Developed at Google Brain, it keeps float32’s $8$ exponent bits (range up to $3.4\times10^{38}$ ) and gives up fraction bits ( $7$ stored). Kalamkar and colleagues (2019) found that training with bfloat16 tensors matched float32 results without changing hyperparameters: for training, range matters more than the last few digits.

Notice the connection to this chapter: each low-precision format is a choice of the unit roundoff $u$ , and the question ‘how many digits of the answer survive?’ is answered by the condition number of 3.5.

Used in the real world· Why np.polyfit warns about conditioning at high degree

np.polyfit(x, y, d) finds the degree- $d$ polynomial that fits the data best in least squares (2.10) by solving a linear system whose matrix is the Vandermonde matrix of this chapter. Its columns $1, x, x^2, \dots, x^d$ look more and more alike as $d$ grows, so its condition number explodes. For $21$ equally spaced points on $[0, 1]$ (degree $20$ ) the Vandermonde matrix has condition number about $9\times10^{16}$ , which is about $1/u$ : by the rule of thumb of 3.5 (you lose $\log_{10}\kappa$ digits) there are no reliable digits left in the coefficients, and polyfit raises RankWarning: Polyfit may be poorly conditioned. With $50$ data points on $[0, 1]$ , the warning appears at degree $20$ and not at degree $15$ . Rescaling $x$ to $[-1, 1]$ already reduces the condition number to about $8\times10^8$ for the same $21$ points, and using the Chebyshev polynomials as the basis (the numpy.polynomial.chebyshev module) keeps it small at any degree. Two lessons. A high-degree polynomial fit of noisy data oscillates between points for the same reason as Runge’s example. And coefficients in the monomial basis are a bad way to store a polynomial: store the values, or the Chebyshev coefficients.

Python labRounding in action: the surprises of floating point

import numpy as np
import warnings

print(0.1 + 0.2, 0.1 + 0.2 == 0.3)
print(np.float32(16777216) + np.float32(1))            # 2**24 + 1 in float32
x = 1e-8
print((1 - np.cos(x)) / x**2, 2 * np.sin(x / 2)**2 / x**2)
print(np.finfo(np.float16).max, np.finfo(np.float32).max, float(np.finfo(np.float64).max))
print(np.float16(0.1).item(), np.float32(0.1).item())
print("h        error of (sin(1+h) - sin(1)) / h  as an estimate of cos(1)")
for h in (1e-2, 1e-4, 1e-8, 1e-12):
    print(f"{h:<8g} {abs((np.sin(1 + h) - np.sin(1)) / h - np.cos(1)):.2e}")
xs = np.linspace(0, 1, 50); ys = np.sin(2 * np.pi * xs)
for d in (10, 15, 20):
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always"); np.polyfit(xs, ys, d)
    print("polyfit degree", d, "warnings:", [str(i.message) for i in w])

Try this. Read each line against the chapter. 0.1 + 0.2 is not 0.3 (not exactly representable). In float32 adding 1 to 2^24 does nothing. The naive (1 - cos x)/x^2 gives 0.0 and the rewritten one 0.5. The finite-difference error is smallest near h = 1e-8, as computed above. Then change x = 1e-8 to 1e-4 and see that the naive formula is now correct to about 8 digits, because the cancellation is less severe. Note that float16 max prints as 65500.0: NumPy rounds the display, the true value is 65504.

Traps

Watch out· More interpolation nodes do not mean a better fit

With equally spaced nodes the error can grow with the degree (Runge). The error formula has the factor $f^{(n+1)}/(n+1)!$ and the factor $\max|w|$ . If the first grows faster than the second shrinks, the error grows. Choosing the nodes (Chebyshev) fixes the second factor, and only smoothness of $f$ fixes the first.

Watch out· The point ξ in the error formula is not computable

The theorem says some $\xi$ exists. You cannot solve for it, so the formula gives bounds, never exact errors. The bound uses the maximum of $|f^{(n+1)}|$ over the whole interval, and is usually several times too large (the true errors in the examples were about $1.5$ to $3$ times smaller).

Watch out· Simpson's rule needs an even number of pieces, and the weights alternate 4, 2, 4, 2

The composite weights are $1, 4, 2, 4, 2, \dots, 4, 1$ times $h/3$ : each odd node is the middle of a pair and gets $4$ ; each even interior node is shared between two pairs and gets $2$ . Using $1, 4, 1$ over and over double-counts the shared ends.

Watch out· Order of accuracy is a statement about smooth functions

‘Trapezoid error shrinks 4 times when $h$ halves; Simpson 16 times’ needs $f''$ and $f^{(4)}$ to be bounded. For a function with a corner or an integrable singularity the observed ratios are smaller. And the ratios cannot be seen once the error reaches the rounding level $10^{-16}$ : more pieces then add rounding error instead of removing truncation error.

Watch out· Shrinking the step forever makes finite differences worse

Truncation error falls with $h$ ; rounding error grows like $u/h$ . The best step is about $\sqrt u \approx 10^{-8}$ for a forward difference and the best error is not zero. A step of $10^{-14}$ is thousands of times worse than $10^{-8}$ .

Watch out· Never test floating-point numbers for equality, and subtract nearly equal numbers with care

0.1 + 0.2 == 0.3 is False. Compare with a tolerance. A formula that subtracts two nearly equal quantities (the smaller root of a quadratic when $b^2 \gg 4ac$ , $1 - \cos x$ , $\sqrt{x+1} - \sqrt x$ for large $x$ ) should be rewritten algebraically before it is coded.

Watch out· Rounding and truncation are two different errors

Truncation error (replacing $f$ by a polynomial, an integral by a sum) goes to zero as the method is refined. Rounding error does not: it sets a floor. When a result stops improving as you refine, you have reached the floor, and refining further helps nothing.

Problems

check

For the nodes $0, 1, 3$ , write the Lagrange basis polynomial $\ell_1(x)$ and find $\ell_1(2)$ .

Solution

$\ell_1(x) = \dfrac{(x - 0)(x - 3)}{(1 - 0)(1 - 3)} = -\dfrac{x(x-3)}{2}$ . At $x = 2$ : $-\dfrac{2\cdot(-1)}{2} = 1$ . (Check: $\ell_1(1) = -\frac{1\cdot(-2)}{2} = 1$ and $\ell_1(0) = \ell_1(3) = 0$ .)

check

Find the determinant of the Vandermonde matrix for the nodes $1, 2, 3, 5$ .

Solution

The product of all differences $x_j - x_i$ with $i < j$ : $(2-1)(3-1)(5-1)(3-2)(5-2)(5-3) = 1\cdot2\cdot4\cdot1\cdot3\cdot2 = 48$ .

check

Apply the trapezoidal rule with $N = 2$ pieces to $\int_0^2 x^2\,dx$ . Compare with the exact value and with the error formula.

Solution

$h = 1$ , nodes $0, 1, 2$ : $T_2 = 1\cdot(\tfrac12\cdot0 + 1 + \tfrac12\cdot4) = 3$ . The exact value is $\tfrac83$ , so the error $T_2 - \int = \tfrac13$ . The formula: $\int - T = -\dfrac{(b-a)h^2}{12}f''(\eta) = -\dfrac{2\cdot 1}{12}\cdot 2 = -\tfrac13$ (here $f'' = 2$ is constant, so the formula is exact).

check

Take one Newton step for $x^2 - 3 = 0$ from $x_0 = 2$ . Check it against the error formula.

Solution

$x_1 = \dfrac{x_0}{2} + \dfrac{3}{2x_0} = 1 + \dfrac34 = 1.75$ . The error is $1.75 - \sqrt3 \approx 0.0179$ . The formula: $\dfrac{(x_0 - \sqrt3)^2}{2x_0} = \dfrac{(0.2679)^2}{4} = 0.0179$ . They agree.

use

Find the Newton form of the polynomial through $(0, 1)$ , $(1, 3)$ , $(3, 13)$ using divided differences, and evaluate it at $x = 2$ .

Solution

$f[0,1] = \dfrac{3 - 1}{1} = 2$ , $f[1,3] = \dfrac{13 - 3}{2} = 5$ , $f[0,1,3] = \dfrac{5 - 2}{3 - 0} = 1$ . So $p(x) = 1 + 2x + 1\cdot x(x - 1) = x^2 + x + 1$ . At $x = 2$ : $7$ . (Check: $p(3) = 13$ .)

use

Bound the error of interpolating $\sin x$ on $[0, \pi]$ at the four equally spaced nodes $0, \tfrac\pi3, \tfrac{2\pi}3, \pi$ . The true largest error is about $0.0436$ ; is your bound consistent?

Solution

$n = 3$ , $h = \pi/3$ , and $|f^{(4)}| = |\sin| \le 1 = M_4$ . The corollary gives $\dfrac{M_4h^4}{4\cdot4} = \dfrac{(\pi/3)^4}{16} \approx \dfrac{1.2026}{16} \approx 0.0752$ . The true error $0.0436$ is below it, so the bound is consistent (and about $1.7$ times too large).

use

Use Simpson’s rule with $N = 4$ pieces for $\int_0^\pi \sin x\,dx$ . Compare with the exact value and with the error bound.

Solution

$h = \pi/4$ . $S_4 = \dfrac{\pi}{12}\big(0 + 4\sin\tfrac\pi4 + 2\sin\tfrac\pi2 + 4\sin\tfrac{3\pi}4 + 0\big) = \dfrac{\pi}{12}\big(4\sqrt2 + 2\big) \approx \dfrac{\pi}{12}(7.6569) \approx 2.00456$ . The exact value is $2$ , so the error is $0.00456$ . The bound: $\dfrac{(b-a)h^4}{180}M_4 = \dfrac{\pi\,(\pi/4)^4}{180} \approx 0.00664$ . The actual error is below it.

use

The smaller root of $x^2 - 10^8x + 1 = 0$ is very close to $10^{-8}$ . Computing it as $\big(10^8 - \sqrt{10^{16} - 4}\big)/2$ in double precision gives about $7.45\times10^{-9}$ . Explain the error and give a better formula.

Solution

$\sqrt{10^{16} - 4} \approx 10^8 - 2\times10^{-8}$ , so the subtraction $10^8 - \sqrt{\cdot}$ cancels about $16$ digits and leaves the rounding error of the square root (the spacing of doubles near $10^8$ is $1.5\times10^{-8}$ , comparable with the difference $2\times10^{-8}$ that we are trying to compute). The relative cancellation factor is $2\times10^8/(2\times10^{-8}) = 10^{16} \approx 1/u$ . Fix: the product of the roots is $c/a = 1$ , so the small root is $1/r_{\text{big}}$ with $r_{\text{big}} = \big(10^8 + \sqrt{10^{16} - 4}\big)/2$ , an addition, with no cancellation. This gives $1.0\times10^{-8}$ to full precision.

prove

Prove that the Lagrange basis polynomials add up to $1$ : $\sum_{j=0}^{n}\ell_j(x) = 1$ for every $x$ . Deduce that if $f$ is a polynomial of degree at most $n$ then its interpolant $p_n$ equals $f$ .

Solution

The constant function $1$ has degree $0 \le n$ and takes the value $1$ at every node. By the existence formula, $\sum_j 1\cdot\ell_j$ is an interpolating polynomial of degree at most $n$ with these values, and by uniqueness it is equal to the constant $1$ . For the second part, if $f$ has degree at most $n$ then $f$ itself is a polynomial of degree at most $n$ with $f(x_j) = f(x_j)$ at every node, so by uniqueness $p_n = f$ . (This agrees with the error formula: $f^{(n+1)} = 0$ .) $\blacksquare$

prove

Prove that Simpson’s rule on $[-h, h]$ is exact for every cubic polynomial, without using the error theorem.

Solution

A cubic is $c_0 + c_1t + c_2t^2 + c_3t^3$ , and both sides of the rule are linear in $f$ , so it is enough to test $1, t, t^2, t^3$ . Write $S = \dfrac h3\big(f(-h) + 4f(0) + f(h)\big)$ . For $f = 1$ : $S = \tfrac h3\cdot 6 = 2h = \int_{-h}^{h}1\,dt$ . For $f = t$ : $S = \tfrac h3(-h + 0 + h) = 0 = \int t\,dt$ . For $f = t^2$ : $S = \tfrac h3(h^2 + 0 + h^2) = \tfrac{2h^3}3 = \int t^2\,dt$ . For $f = t^3$ : $S = \tfrac h3(-h^3 + 0 + h^3) = 0 = \int t^3\,dt$ (odd function on a symmetric interval). So the rule is exact for $1, t, t^2, t^3$ and hence for every cubic. By translation the same holds on any interval. $\blacksquare$

prove

Prove that if $f'' \ge 0$ on $[a, b]$ (a convex function) then the trapezoidal rule overestimates: $T \ge \int_a^b f$ .

Solution

By the error theorem, $\int_a^b f - T = -\dfrac{(b-a)^3}{12}f''(\eta)$ for some $\eta \in (a, b)$ . Since $f''(\eta) \ge 0$ and $(b - a)^3 > 0$ , the right side is $\le 0$ , so $\int_a^b f \le T$ . $\blacksquare$ (Picture: the chord lies above a convex curve, so the trapezoid under the chord has more area.)

prove

For the Chebyshev nodes $x_0, \dots, x_n$ show that $w(x) = \prod_j(x - x_j) = T_{n+1}(x)/2^n$ , and conclude that $\max_{[-1,1]}|w| = 2^{-n}$ .

Solution

The nodes are exactly the $n+1$ distinct roots of $T_{n+1}$ , a polynomial of degree $n + 1$ . So $T_{n+1}(x) = c\prod_j(x - x_j)$ where $c$ is the leading coefficient of $T_{n+1}$ , which is $2^n$ . Hence $w = T_{n+1}/2^n$ . Since $|T_{n+1}| \le 1$ on $[-1, 1]$ with the value $\pm1$ attained (at $x = 1$ ), $\max|w| = 2^{-n}$ . $\blacksquare$

prove

Newton’s method for $x^2 - a$ ( $a > 0$ ) is $x_{k+1} = \tfrac12(x_k + a/x_k)$ . Prove that if $x_0 > 0$ then $x_k \ge \sqrt a$ for all $k \ge 1$ , and that $x_{k+1} \le x_k$ for $k \ge 1$ . Deduce convergence to $\sqrt a$ .

Solution

By the arithmetic-geometric mean inequality, $\tfrac12(x + a/x) \ge \sqrt{x\cdot a/x} = \sqrt a$ for any $x > 0$ (equivalently, $\tfrac12(x + a/x) - \sqrt a = \dfrac{(x - \sqrt a)^2}{2x} \ge 0$ ). So $x_k \ge \sqrt a$ for $k \ge 1$ . Then $x_k - x_{k+1} = \tfrac12\big(x_k - a/x_k\big) = \dfrac{x_k^2 - a}{2x_k} \ge 0$ . So $(x_k)_{k \ge 1}$ is decreasing and bounded below by $\sqrt a$ ; by the monotone convergence theorem (1.2) it converges to a limit $L \ge \sqrt a$ , and passing to the limit in the recursion gives $L = \tfrac12(L + a/L)$ , i.e. $L^2 = a$ , so $L = \sqrt a$ . $\blacksquare$

stretch

Find weights $w_0, w_1, w_2$ so that $\int_{-1}^{1} f \approx w_0f(-1) + w_1f(0) + w_2f(1)$ is exact for $f = 1, t, t^2$ (the method of undetermined coefficients). Show that the rule is then also exact for $t^3$ , and identify it.

Solution

Exactness for $1$ : $w_0 + w_1 + w_2 = 2$ . For $t$ : $-w_0 + w_2 = 0$ . For $t^2$ : $w_0 + w_2 = \tfrac23$ . The last two give $w_0 = w_2 = \tfrac13$ , then $w_1 = 2 - \tfrac23 = \tfrac43$ . This is Simpson’s rule with $h = 1$ . For $t^3$ : the rule gives $\tfrac13(-1) + 0 + \tfrac13(1) = 0$ , and $\int_{-1}^1 t^3 = 0$ . It is exact for all cubics, which is one degree more than the three conditions forced. (The method of undetermined coefficients is how rules of higher order are derived.)

stretch

Suppose $e_{k+1} \le K e_k^2$ for all $k$ and $K e_0 \le \tfrac12$ . Prove that $e_k \le \dfrac1K\left(\tfrac12\right)^{2^k}$ . How many steps until $e_k < 10^{-16}/K$ ?

Solution

Let $a_k = Ke_k$ . Multiplying the hypothesis by $K$ gives $a_{k+1} \le a_k^2$ . By induction $a_k \le a_0^{2^k}$ : true for $k = 0$ , and if $a_k \le a_0^{2^k}$ then $a_{k+1} \le a_k^2 \le a_0^{2^{k+1}}$ (all quantities non-negative). Since $a_0 \le \tfrac12$ , $a_k \le (\tfrac12)^{2^k}$ , i.e. $e_k \le \tfrac1K(\tfrac12)^{2^k}$ . We need $(\tfrac12)^{2^k} < 10^{-16}$ , that is $2^k > 16\log_2 10 \approx 53.2$ , so $k = 6$ steps ( $2^6 = 64$ ). Six Newton steps from a start with $Ke_0 \le \frac12$ reach double precision.

Summary

The interpolating polynomial of degree at most $n$ through $n + 1$ points with distinct nodes exists and is unique. Existence: the Lagrange form $\sum y_j\ell_j(x)$ . Uniqueness: two interpolants differ by a polynomial with too many roots. The coefficient system has the Vandermonde matrix, with determinant $\prod_{i<j}(x_j - x_i)$ .
The error formula is $f(x) - p_n(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_j(x - x_j)$ , proved by Rolle’s theorem applied $n + 1$ times. It is Taylor’s remainder with $n + 1$ centres. Equally spaced nodes give the bound $M_{n+1}h^{n+1}/(4(n+1))$ .
Runge’s phenomenon: for $1/(1 + 25x^2)$ the equally spaced interpolation error grows with $n$ (from $0.6$ at $n = 6$ to $29$ at $n = 18$ ). Chebyshev nodes $\cos\frac{(2k+1)\pi}{2n+2}$ make $\max|w| = 2^{-n}$ , the smallest possible (Chebyshev’s minimal property), and the error shrinks.
Newton’s form adds one divided difference per node and works as a table. A divided difference of order $n$ equals $f^{(n)}(\xi)/n!$ . Taylor’s polynomial is the case of one repeated node.
Quadrature integrates the interpolant. The trapezoidal rule has error $-\frac{(b-a)^3}{12}f''(\eta)$ and Simpson’s $-\frac{(b-a)^5}{2880}f^{(4)}(\eta)$ ; composite versions have errors proportional to $h^2$ and $h^4$ , so doubling the pieces divides the error by $4$ and $16$ .
A floating-point number has $p$ significant bits and unit roundoff $u = 2^{-p}$ ( $1.1\times10^{-16}$ in double). Every operation is correct up to a factor $1 + \delta$ , $|\delta| \le u$ . Subtracting nearly equal numbers (catastrophic cancellation) multiplies the relative error by $(|x| + |y|)/|x - y|$ . Truncation error shrinks with the method, rounding error sets a floor.
Newton’s method has $x_{k+1} - r = \frac{f''(\xi)}{2f'(x_k)}(x_k - r)^2$ : quadratic convergence, correct digits double each step, until rounding stops it.
Used: the choice of float32, float16 or bfloat16 in machine learning, the conditioning warning of np.polyfit, finite-difference step sizes, and, inside maths, Taylor’s theorem as repeated-node interpolation and Chebyshev polynomials as an orthogonal family.

Watch and read

The same ideas explained by someone else. Use these when an explanation here does not land.

The picture

The idea, in words

Definitions and theorems

Existence and uniqueness

Interpolating polynomial

Lagrange basis polynomial

Existence and uniqueness of the interpolating polynomial

Vandermonde matrix

The Vandermonde determinant

The error formula (Apostol 15.8)

Interpolation error formula

Error bound for equally spaced nodes

Runge’s phenomenon and Chebyshev nodes (Apostol 15.14 to 15.17)

Chebyshev polynomial

Properties of the Chebyshev polynomials

Chebyshev nodes

Chebyshev's minimal property

Error bound for Chebyshev nodes

Newton’s divided differences (Apostol 15.10)

Divided difference

Newton's divided-difference formula

Numerical integration (Apostol 15.19 and 15.20)

Quadrature rule

Trapezoidal rule

Simpson's rule

Cauchy's mean value theorem

Error of the trapezoidal and Simpson rules

Composite trapezoidal and Simpson rules

Rounding: how a computer stores numbers

Floating-point number

Catastrophic cancellation

Newton’s method, again (link to 4.2)

Newton's method converges quadratically

Examples and non-examples

Where it is used

Traps

Problems

Summary

Watch and read

Where this sits

Comes from

Leads to

Same idea elsewhere