Chapter 4 — Numerical Differentiation and Integration

Differentiation and integration are both limits, and a computer cannot take a limit. The strategy throughout this chapter is the same one as in Chapter 3: replace \(f\) by a polynomial that interpolates it, then differentiate or integrate the polynomial exactly. Where the two operations differ is in how they behave. Differentiation is unstable — it divides by a small number, so round-off eventually dominates — while integration is stable and averages errors out. That asymmetry is why §4.1 is short and §4.3 onwards is long.

4.1 Numerical Differentiation

Forward and backward differences

Start with the smallest possible case: two nodes, \(x_0\) and \(x_1 = x_0 + h\). The Lagrange polynomial with its error term is \[ f(x) = f(x_0)\frac{x - x_0 - h}{-h} + f(x_0+h)\frac{x - x_0}{h} + \frac{(x-x_0)(x-x_0-h)}{2}f''(\xi(x)). \]

Differentiating with respect to \(x\), \[ f'(x) = \frac{f(x_0+h) - f(x_0)}{h} + \frac{2(x-x_0) - h}{2}f''(\xi(x)) + \frac{(x-x_0)(x-x_0-h)}{2}\,D_x\!\left[f''(\xi(x))\right]. \]

Setting \(x = x_0\) makes the last term vanish, because the factor \((x-x_0)\) is zero there — which is convenient, since we know nothing about the derivative of \(\xi(x)\). What remains is the forward-difference formula: \[ f'(x_0) = \frac{f(x_0+h) - f(x_0)}{h} - \frac{h}{2}f''(\xi). \]

Taking \(h > 0\) gives the forward difference; \(h < 0\) gives the backward difference. Either way the error is \(O(h)\): first order, and not very good.

Three-point formulas

Using three nodes instead of two gives a better formula, by the same route: write down the quadratic interpolant and differentiate it.

Differentiating the Lagrange form at a node \(x_j\) gives \[ \begin{aligned} f'(x_j) = f(x_0)\frac{2x_j - x_1 - x_2}{(x_0-x_1)(x_0-x_2)} &+ f(x_1)\frac{2x_j - x_0 - x_2}{(x_1-x_0)(x_1-x_2)} \\ &+ f(x_2)\frac{2x_j - x_0 - x_1}{(x_2-x_0)(x_2-x_1)} + \frac{1}{6}f^{(3)}(\xi_j)\prod_{k\ne j}(x_j - x_k). \end{aligned} \]

This is general, but the useful cases are equally spaced nodes, \(x_1 - x_0 = x_2 - x_1 = h\). Substituting and simplifying at \(x_j = x_0\), \[ \begin{aligned} f'(x_0) &= \frac{1}{h}\left[-\tfrac32 f(x_0) + 2f(x_1) - \tfrac12 f(x_2)\right] + \frac{h^2}{3}f^{(3)}(\xi_0) \\ &= \frac{1}{h}\left[-\tfrac32 f(x_0) + 2f(x_0+h) - \tfrac12 f(x_0+2h)\right] + \frac{h^2}{3}f^{(3)}(\xi_0), \end{aligned} \] which is the three-point endpoint formula. Replacing \(h\) by \(-h\) gives the version that looks backwards: \[ f'(x_0) = \frac{1}{h}\left[\tfrac12 f(x_0-2h) - 2f(x_0-h) + \tfrac32 f(x_0)\right] + \frac{h^2}{3}f^{(3)}(\xi_0). \]

Evaluating at the middle node instead gives the three-point midpoint formula, in which the \(f(x_0)\) term drops out altogether: \[ f'(x_0) = \frac{1}{2h}\left[f(x_0+h) - f(x_0-h)\right] - \frac{h^2}{6}f^{(3)}(\xi_1). \]

Both are \(O(h^2)\), but the midpoint formula has half the error constant and uses one fewer evaluation. It is the one to reach for when the data allows it.

Second derivatives

For higher derivatives it is easier to work directly from Taylor’s theorem than from an interpolant. Expanding in both directions about \(x_0\): \[ \begin{aligned} f(x_0+h) &= f(x_0) + f'(x_0)h + \tfrac12 f''(x_0)h^2 + \tfrac16 f^{(3)}(x_0)h^3 + \tfrac{1}{24}f^{(4)}(\xi_1)h^4, \\ f(x_0-h) &= f(x_0) - f'(x_0)h + \tfrac12 f''(x_0)h^2 - \tfrac16 f^{(3)}(x_0)h^3 + \tfrac{1}{24}f^{(4)}(\xi_{-1})h^4 . \end{aligned} \]

Adding them cancels every odd-order term: \[ f(x_0+h) + f(x_0-h) = 2f(x_0) + f''(x_0)h^2 + \tfrac{1}{24}\left[f^{(4)}(\xi_1) + f^{(4)}(\xi_{-1})\right]h^4 . \]

The bracket is an average of two values of \(f^{(4)}\), so by the intermediate value theorem it equals \(2f^{(4)}(\xi)\) for some \(\xi\) between them. Solving for \(f''(x_0)\) gives the centred second-difference formula: \[ f''(x_0) = \frac{1}{h^2}\left[f(x_0+h) - 2f(x_0) + f(x_0-h)\right] - \frac{h^2}{12}f^{(4)}(\xi). \]

The bound above suggests taking \(h\) as small as possible. Round-off says otherwise: the numerator is a difference of nearly equal numbers divided by \(h^2\), so the cancellation of §1.2 is amplified as \(h\) shrinks. The total error is the sum of a truncation term that decreases with \(h\) and a round-off term that increases, and the best \(h\) is where the two balance — not zero.

4.2 Richardson’s Extrapolation

Richardson’s extrapolation is a way of getting a high-order formula out of a low-order one, for free, by combining results at two step sizes. Suppose an approximation \(N(h)\) to an unknown \(M\) has an error expansion \[ M = N(h) + K_1h + K_2h^2 + K_3h^3 + \cdots \tag{$*$} \] with constants \(K_i\) we do not know. Halving the step, \[ M = N\!\left(\tfrac{h}{2}\right) + K_1\frac{h}{2} + K_2\frac{h^2}{4} + K_3\frac{h^3}{8} + \cdots \tag{$**$} \]

Now form \(2\cdot(**) - (*)\). The \(K_1\) terms cancel exactly: \[ M = \underbrace{2N\!\left(\tfrac{h}{2}\right) - N(h)}_{\textstyle N_2(h)} + K_2\!\left(\frac{h^2}{2} - h^2\right) + K_3\!\left(\frac{h^3}{4} - h^3\right) + \cdots \]

So defining \[ N_2(h) = N\!\left(\tfrac{h}{2}\right) + \left[N\!\left(\tfrac{h}{2}\right) - N(h)\right] \] gives an \(O(h^2)\) approximation built from two \(O(h)\) ones. The bracketed term is a computable estimate of the error in \(N(h/2)\), which is what makes the trick useful in practice as well as in theory.

The same step applies again. Writing the expansion for \(N_2\), \[ M = N_2(h) - \frac{K_2}{2}h^2 - \frac{3K_3}{4}h^3 + \cdots, \qquad M = N_2\!\left(\tfrac{h}{2}\right) - \frac{K_2}{2}\cdot\frac{h^2}{4} - \frac{3K_3}{4}\cdot\frac{h^3}{8} + \cdots \] and forming \(4\cdot(\text{second}) - (\text{first})\) eliminates \(K_2\): \[ N_3(h) = N_2\!\left(\tfrac{h}{2}\right) + \tfrac13\left[N_2\!\left(\tfrac{h}{2}\right) - N_2(h)\right], \] which is \(O(h^3)\). In general, \[ N_j(h) = N_{j-1}\!\left(\tfrac{h}{2}\right) + \frac{N_{j-1}\!\left(\tfrac{h}{2}\right) - N_{j-1}(h)}{2^{\,j-1} - 1}. \]

The computation is organised as a tableau, each column one order better than the last, and each new entry built from two entries to its left:

Each row costs one new evaluation of \(N\) at half the previous step; everything else is arithmetic.

4.3 Elements of Numerical Integration

The recipe is the same as for differentiation: interpolate, then integrate the interpolant exactly. Choose nodes \(x_0, \ldots, x_n\) and write \[ P_n(x) = \sum_{i=0}^n f(x_i)L_i(x). \]

Integrating \(f\) in its Lagrange-plus-remainder form, \[ \int_a^b f(x)\,dx = \int_a^b \sum_{i=0}^n f(x_i)L_i(x)\,dx + \int_a^b \prod_{i=0}^n (x - x_i)\frac{f^{(n+1)}(\xi(x))}{(n+1)!}\,dx, \] and pulling the constants out of the first integral gives a quadrature rule \[ \int_a^b f(x)\,dx = \sum_{i=0}^n a_i f(x_i) + E(f), \qquad a_i = \int_a^b L_i(x)\,dx . \]

The weights are the integrals of the basis functions — nothing more. And because the rule reproduces \(\int P_n\) exactly, it is exact for every polynomial of degree \(\le n\).

The trapezoidal rule

Take \(n = 1\), with \(x_0 = a\), \(x_1 = b\) and \(h = x_1 - x_0\). The basis functions are the two lines \[ L_0(x) = \frac{x - x_1}{x_0 - x_1}, \qquad L_1(x) = \frac{x - x_0}{x_1 - x_0}, \] and each weight is the area of a triangle of base \(h\) and height \(1\): \[ a_0 = \int_{x_0}^{x_1} L_0(x)\,dx = \frac{h}{2}, \qquad a_1 = \frac{h}{2}. \]

For the error, \(g(x) = (x-x_0)(x-x_1)\) does not change sign on \([x_0,x_1]\), so the weighted mean value theorem for integrals (Theorem 1.6) lets us pull \(f''\) out: \[ E(f) = \frac12\int_{x_0}^{x_1} f''(\xi(x))\,(x-x_0)(x-x_1)\,dx = \frac{f''(\xi)}{2}\int_{x_0}^{x_1}(x-x_0)(x-x_1)\,dx . \] Evaluating that integral gives \(-h^3/6\), so \[ \int_a^b f(x)\,dx = \frac{h}{2}\left[f(x_0) + f(x_1)\right] - \frac{h^3}{12}f''(\xi). \]

Note the consequence: the rule is exact whenever \(f'' \equiv 0\), that is for all linear functions. Its degree of precision is \(1\).

Simpson’s rule

Taking \(n = 2\) with equally spaced nodes gives Simpson’s rule. Rather than integrate three basis functions, it is quicker to impose exactness and solve for the weights.

The error involves \(f^{(4)}\), so Simpson’s rule is exact for cubics as well as quadratics: its degree of precision is \(3\), one better than interpolation alone would suggest. This bonus order is why Simpson’s rule, and not the quadratic rule one might expect, is the workhorse of elementary quadrature.

4.4 Composite Numerical Integration

The obvious way to improve accuracy is to use more nodes in a single rule. It does not work well.

The coefficients become unwieldy, some of them turn negative — which invites cancellation — and the interpolant oscillates, exactly as in the Runge example of §3.5. The fix is the same as it was there: go piecewise.

The composite trapezoidal rule

Divide \([a,b]\) into \(n\) equal subintervals of width \(h = (b-a)/n\), with \(x_j = a + jh\), and apply the trapezoidal rule on each: \[ \begin{aligned} \int_a^b f(x)\,dx &\approx \frac{h}{2}\bigl(f(x_0) + f(x_1)\bigr) + \frac{h}{2}\bigl(f(x_1) + f(x_2)\bigr) + \cdots + \frac{h}{2}\bigl(f(x_{n-1}) + f(x_n)\bigr) \\ &= \frac{h}{2}\left[f(x_0) + 2\sum_{j=1}^{n-1}f(x_j) + f(x_n)\right]. \end{aligned} \]

Its error works out to \[ E(f) = -\frac{b-a}{12}h^2 f''(\mu) \] for some \(\mu \in (a,b)\): second order in \(h\).

The composite Simpson rule

Simpson’s rule consumes two subintervals at a time, so \(n\) must be even. Applying it to each consecutive pair and collecting terms, \[ \begin{aligned} \int_a^b f(x)\,dx &\approx \frac{h}{3}\left[f(x_0) + 4f(x_1) + f(x_2)\right] + \frac{h}{3}\left[f(x_2) + 4f(x_3) + f(x_4)\right] + \cdots \\ &= \frac{h}{3}\left[f(x_0) + 2\sum_{j=1}^{n/2-1}f(x_{2j}) + 4\sum_{j=1}^{n/2}f(x_{2j-1}) + f(x_n)\right]. \end{aligned} \]

The error is a sum of the individual panel errors, \[ E(f) = -\sum_{j=1}^{n/2}\frac{h^5}{90}f^{(4)}(\xi_j), \qquad \xi_j \in (x_{2j-2}, x_{2j}), \] and this can be collapsed into a single term. Since every \(f^{(4)}(\xi_j)\) lies between the minimum and maximum of \(f^{(4)}\) on \([a,b]\), so does their average: \[ \min_{x\in[a,b]} f^{(4)}(x) \;\le\; \frac{2}{n}\sum_{j=1}^{n/2}f^{(4)}(\xi_j) \;\le\; \max_{x\in[a,b]} f^{(4)}(x), \] so by the intermediate value theorem there is a \(\mu \in (a,b)\) with \[ f^{(4)}(\mu) = \frac{2}{n}\sum_{j=1}^{n/2}f^{(4)}(\xi_j). \] Substituting, and using \(n = (b-a)/h\), \[ E(f) = -\frac{h^5}{90}\cdot\frac{n}{2}\cdot f^{(4)}(\mu) = -\frac{h^4}{180}(b-a)f^{(4)}(\mu) = -\frac{b-a}{180}h^4 f^{(4)}(\mu), \] fourth order: two orders better than the composite trapezoidal rule for the same number of evaluations.

Example 4.4

Composite Simpson on \(\displaystyle\int_0^\pi \sin^2 x\,dx\) with \(n = 4\).

Here \(h = \pi/4\) and the nodes are \(0, \pi/4, \pi/2, 3\pi/4, \pi\), with values \(0, \tfrac12, 1, \tfrac12, 0\).

Applying the rule, \[ \begin{aligned} I &\approx \frac{h}{3}\left[f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4)\right] \\ &= \frac{\pi}{4\cdot 3}\left[0 + 4\cdot\tfrac12 + 2\cdot 1 + 4\cdot\tfrac12 + 0\right] = \frac{\pi}{12}\cdot 6 = \frac{\pi}{2}. \end{aligned} \] This is the exact answer — a coincidence, not a general phenomenon.

The bound does not know that. With \(f(x) = \sin^2 x = \tfrac{1-\cos 2x}{2}\) we get \(f^{(4)}(x) = -\tfrac{2^4\cos 2x}{2} = -8\cos 2x\), so \(|f^{(4)}| \le 8\) and \[ |E(f)| \le \frac{b-a}{180}h^4\left|f^{(4)}(\mu)\right| \le \frac{\pi}{180}\left(\frac{\pi}{4}\right)^{4}\cdot 8 = 0.053 . \]

4.5 Romberg Integration

Romberg integration is Richardson’s extrapolation applied to the composite trapezoidal rule. What makes the pairing so effective is that the trapezoidal error expansion contains only even powers of \(h\), so each extrapolation step gains two orders rather than one.

Write \(R_{k,1}\) for the composite trapezoidal estimate with \(h_k = (b-a)/2^{k-1}\). Successive levels nest: level \(1\) uses only \(a\) and \(b\), level \(2\) adds the midpoint, level \(3\) adds the two quarter points, and so on. Each level therefore costs only the new points: \[ R_{k,1} = \frac12\left[R_{k-1,1} + h_{k-1}\sum_{i=1}^{2^{k-2}} f\!\left(a + \left(i - \tfrac12\right)h_{k-1}\right)\right]. \]

Extrapolating across each row, \[ R_{k,j} = R_{k,j-1} + \frac{R_{k,j-1} - R_{k-1,j-1}}{4^{\,j-1} - 1}, \] where the \(4^{j-1}\) — rather than \(2^{j-1}\) — is precisely the even-powers bonus. The entry \(R_{k,j}\) is \(O(h_k^{2j})\).

function R = romberg(f, a, b, n)
% Compute integral of f(x) from a to b using Romberg integration.

h = b - a;
R = zeros(n, n);
R(1, 1) = h / 2 * (f(a) + f(b));
for i = 2:n
    x = a + ((1:2^(i - 2)) - 0.5) * h;
    R(i, 1) = (R(i - 1, 1) + h * sum(f(x))) / 2;
    for j = 2:i
        R(i, j) = R(i, j - 1) + ...
            (R(i, j - 1) - R(i - 1, j - 1)) / (4^(j - 1) - 1);
    end
    h = h / 2;
end

4.6 Adaptive Quadrature Methods

Composite rules spread their nodes uniformly, which is wasteful when \(f\) is smooth in most of the interval and difficult in a small part of it.

To place points intelligently we need to know where the error is large — but the error formula involves \(f^{(4)}(\mu)\), about which we know nothing in practice. The way out is to estimate the error from two computations at different step sizes.

Let \(S(a,b)\) denote Simpson’s rule applied on \([a,b]\). Then \[ \int_a^b f(x)\,dx = S(a,b) - \frac{h^5}{90}f^{(4)}(\mu), \] and halving the interval, applying Simpson to each half, \[ \int_a^b f(x)\,dx = S\!\left(a, \tfrac{a+b}{2}\right) + S\!\left(\tfrac{a+b}{2}, b\right) - \frac{1}{16}\cdot\frac{h^5}{90}f^{(4)}(\tilde\mu), \] the \(1/16\) arising because halving \(h\) divides \(h^5\) by \(32\) while there are now two panels.

Assuming \(f^{(4)}(\mu) \approx f^{(4)}(\tilde\mu)\) — the one leap of faith — and subtracting, \[ \frac{h^5}{90}f^{(4)}(\mu) \approx \frac{16}{15} \left[S(a,b) - S\!\left(a, \tfrac{a+b}{2}\right) - S\!\left(\tfrac{a+b}{2}, b\right)\right]. \] The error in the better (two-panel) estimate is \(1/16\) of this, so \[ \text{error} \approx \frac{1}{15} \left[S(a,b) - S\!\left(a, \tfrac{a+b}{2}\right) - S\!\left(\tfrac{a+b}{2}, b\right)\right]. \]

This is computable from quantities we already have. If it exceeds the tolerance, subdivide and recurse; if not, accept the panel.

4.7 Gaussian Quadrature

Every rule so far has used equally spaced nodes, chosen for convenience. If we are willing to choose the nodes too, we get far more accuracy per evaluation.

Choosing the nodes

Suppose instead we allow both nodes and weights to be free: \[ \int_{-1}^{1} f(x)\,dx \approx c_1 f(x_1) + c_2 f(x_2). \] That is four unknowns, so we can hope for exactness on four monomials — degree of precision \(3\). Imposing them: \[ \begin{aligned} f = 1:\quad & 2 = c_1 + c_2, &(1)\\ f = x:\quad & 0 = c_1x_1 + c_2x_2, &(2)\\ f = x^2:\quad & \tfrac23 = c_1x_1^2 + c_2x_2^2, &(3)\\ f = x^3:\quad & 0 = c_1x_1^3 + c_2x_2^3. &(4) \end{aligned} \]

Symmetry suggests \(c_1 = c_2\) and \(x_1 = -x_2\), which satisfies \((2)\) and \((4)\) identically. Then \((1)\) gives \(2 = 2c_1\), so \(c_1 = c_2 = 1\), and \((3)\) becomes \[ \tfrac23 = x_1^2 + x_2^2 = 2x_1^2 \qquad\Longrightarrow\qquad x_1 = -\frac{1}{\sqrt3},\quad x_2 = \frac{1}{\sqrt3}, \] giving the two-point Gauss rule \[ \int_{-1}^{1} f(x)\,dx \approx f\!\left(-\tfrac{1}{\sqrt3}\right) + f\!\left(\tfrac{1}{\sqrt3}\right), \] with degree of precision \(3\) — the same as Simpson’s rule, from two evaluations instead of three.

Legendre polynomials

Solving for the nodes by hand does not scale, but there is a pattern: the optimal nodes are the roots of the Legendre polynomials. \(P_n\) is the monic polynomial of degree \(n\) with \[ \int_{-1}^{1} P(x)P_n(x)\,dx = 0 \qquad\text{whenever } \deg P < n . \]

The first few are \[ P_0(x) = 1, \qquad P_1(x) = x, \qquad P_2(x) = x^2 - \tfrac13, \] and indeed the roots of \(P_2\) are \(\pm 1/\sqrt3\) — the nodes we just derived.

function [x, c] = gaussquad(n)
% Compute Gaussian quadrature points and coefficients.
% For education only: the monomial basis is unstable for high n.

P = zeros(n + 1, n + 1);
P([1, 2], 1) = 1;
for k = 1:n - 1
    P(k + 2, 1:k + 2) = ...
        ((2 * k + 1) * [P(k + 1, 1:k + 1), 0] - ...
        k * [0, 0, P(k, 1:k)]) / (k + 1);
end
x = sort(roots(P(n + 1, 1:n + 1)));

A = zeros(n, n);
for i = 1:n
    A(i, :) = polyval(P(i, 1:i), x)';
end
c = A \ [2; zeros(n - 1, 1)];

Arbitrary intervals

Nodes and weights are tabulated only for \([-1,1]\), so an affine change of variable does the rest.

Setting \[ x = \tfrac12\left[(b-a)t + a + b\right], \qquad\text{equivalently}\qquad t = \frac{2x - a - b}{b - a}, \] we get \[ \int_a^b f(x)\,dx = \int_{-1}^{1}f\!\left(\frac{(b-a)t + (b+a)}{2}\right)\frac{b-a}{2}\,dt . \]

4.8 Multiple Integrals

A double integral over a rectangle is two nested one-dimensional integrals, so the rules carry over directly: apply a rule in \(y\) for each fixed \(x\), then a rule in \(x\).

Take \(R = [a,b]\times[c,d]\) with \(h = (b-a)/n\) and \(k = (d-c)/m\). Writing the inner integral first, \[ \int_a^b\!\left(\int_c^d f(x,y)\,dy\right)dx \approx \int_a^b \frac{k}{2}\Bigl(f(x,c) + 2f(x,c+k) + f(x,d)\Bigr)dx \] for \(m = 2\), and then applying the trapezoidal rule in \(x\) as well and collecting terms gives \[ \begin{aligned} \approx \frac{hk}{4}\Bigl( &f(a,c) + 2f(a{+}h,c) + f(b,c) \\ &+ 2f(a,c{+}k) + 4f(a{+}h,c{+}k) + 2f(b,c{+}k) \\ &+ f(a,d) + 2f(a{+}h,d) + f(b,d)\Bigr). \end{aligned} \]

The weights are products of the one-dimensional weights, which is why they form a recognisable stencil.

Example 4.11

Gaussian double integration on the same integral.

Using the two-point rule in each direction needs only four evaluations.

The maps are \(x = \tfrac12(t+1)\) and \(y = t+1\), so \[ x_1 = \tfrac12\!\left(1 - \tfrac{1}{\sqrt3}\right),\quad x_2 = \tfrac12\!\left(1 + \tfrac{1}{\sqrt3}\right),\quad y_1 = 1 - \tfrac{1}{\sqrt3},\quad y_2 = 1 + \tfrac{1}{\sqrt3}, \] and \[ I \approx \frac{(1-0)(2-0)}{4} \left[e^{x_1+y_1} + e^{x_2+y_1} + e^{x_1+y_2} + e^{x_2+y_2}\right] = 10.940 . \] The exact value is \((e-1)(e^2-1) \approx 11.026\), so four Gauss points beat nine Simpson points.

Non-rectangular regions

When the region is bounded by curves \(c(x) \le y \le d(x)\), the inner limits depend on \(x\), so each vertical slice gets its own grid.

4.9 Improper Integrals

An integrand that blows up at an endpoint defeats every rule so far, since they all evaluate \(f\) there. The integral may still exist as a limit: \[ \int_0^1 \frac{1}{x^{1/2}}\,dx = \lim_{t\to 0}\int_t^1 x^{-1/2}\,dx = \lim_{t\to 0}\left[\frac{x^{1/2}}{1/2}\right]_t^1 = \lim_{t\to 0} 2\left[1 - t^{1/2}\right] = 2 . \]

The trick is to subtract off the part that causes the trouble. Replace the smooth factor by its Taylor polynomial, integrate that piece exactly — the singularity is then only \(x^{-p}\) against a polynomial, which we can do in closed form — and handle the remainder numerically, where it is well behaved.

Infinite limits of integration

An infinite interval is handled by the substitution \(t = 1/x\), which maps \([1,\infty)\) onto \((0,1]\).