Numerical integration¶

In calculus you learn that the elegant way to evaluate a definite integral is to apply the Fundamental Theorem of Calculus and find an antiderivative. The connection is so profound and pervasive that it’s easy to overlook that a definite integral is a numerical quantity existing independently of antidifferentiation. In fact, most conceivable integrands have no antiderivative in terms of familiar functions.

Numerical integration1 is done by combining values of the integrand sampled at nodes, much like finite differences. In this section we will assume equally spaced nodes using the definitions

(147)¶\[ t_i = a +i h, \quad h=\frac{b-a}{n}, \qquad i=0,\ldots,n.\]

The integration formulas are expressed as

(148)¶\[\begin{split} \begin{split} I = \int_a^b f(x)\, dx \approx Q &= h \sum_{i=0}^n w_if(t_i) \\ &= h \bigl[ w_0f(t_0)+w_1f(t_1)+\cdots w_nf(t_n) \bigr]. \end{split}\end{split}\]

The constants \(w_i\) appearing in the formula are called weights. As with finite difference formulas, the weights of numerical integration formulas are chosen independently of the function being integrated, and they determine the formula completely. We can apply quadrature formulas to sequences of data values even if no function is explicitly known to generate them, but for presentation and implementations we assume that we can evaluate \(f(x)\) anywhere.

A straightforward way to derive integration formulas is to mimic the approach taken for finite differences: find an interpolant and operate exactly on it. If the interpolant is a piecewise polynomial, the result is a Newton–Cotes formula.

Trapezoid formula¶

One of the most important Newton–Cotes formulas results from integration of the piecewise linear interpolant (see Piecewise linear interpolation). Using the cardinal basis form of the interpolant in (109), we have

\[I \approx \int_a^b \sum_{i=0}^n f(t_i) H_i(x)\, dx = \sum_{i=0}^n f(t_i) \left[ \int_a^b H_i(x)\right]\, dx.\]

Thus we can identify the weights as \(w_i = h^{-1} \int_a^b H_i(x)\, dx\). Using areas of triangles, it’s trivial to derive that

(149)¶\[\begin{split} w_i = \begin{cases} 1, & i=1,\ldots,n-1,\\ \frac{1}{2}, & i=0,n. \end{cases}\end{split}\]

Putting everything together, the resulting quadrature formula is

(150)¶\[\begin{split} I = \int_a^b f(x)\, dx \approx T_f(n) &= h\left[ \frac{1}{2}f(t_0) + f(t_1) + f(t_2) + \cdots + f(t_{n-1}) + \frac{1}{2}f(t_n) \right]. \end{split}\]

This is called the trapezoid formula or trapezoid rule.2 The trapezoid formula results from integration of the piecewise linear interpolant, or equivalently, as illustrated in Fig. 5, from using the area of approximating trapezoids to estimate the area under a curve. The trapezoid formula is the Swiss Army knife of integration formulas. A short implementation is given as trapezoid.

Fig. 5 Trapezoid formula for integration.¶

Function 55 (trapezoid)

Trapezoid formula for numerical integration.

"""
trapezoid(f,a,b,n)

Apply the trapezoid integration formula for integrand `f` over
interval [`a`,`b`], broken up into `n` equal pieces. Returns
estimate, vector of nodes, and vector of integrand values at the
nodes.
"""
function trapezoid(f,a,b,n)
    h = (b-a)/n
    t = LinRange(a,b,n+1)
    y = f.(t)
    T = h * ( sum(y[2:n]) + 0.5*(y[1] + y[n+1]) )

    return T,t,y
end

In the PL convergence theorem we stated that the pointwise error in a piecewise linear interpolant with equal node spacing \(h\) is bounded by \(O(h^2)\) as \(h\rightarrow 0\). Using \(p\) to stand for the piecewise linear interpolant, we obtain

\[\begin{split}\begin{split} I - T_f(n) = I - \int_a^b p(x)\, dx &= \int_a^b \bigl[f(x)-p(x)\bigr] \, dx \\ &\le (b-a) \max_{x\in[a,b]} |f(x)-p(x)| = O(h^2). \end{split}\end{split}\]

Hence the trapezoid formula has second-order error. This fact is embedded rigorously in “one of the most remarkable formulas in mathematics,” the Euler–Maclaurin formula, which may be stated as

(151)¶\[\begin{split} \begin{split} I = \int_a^b f(x)\, dx &= T_f(n) - \frac{h^2}{12} \left[ f'(b)-f'(a) \right] + \frac{h^4}{740} \left[ f'''(b)-f'''(a) \right] + O(h^6) \\ &= T_f(n) - \sum_{k=1}^\infty \frac{B_{2k}h^{2k}}{(2k)!} \left[ f^{(2k-1)}(b)-f^{(2k-1)}(a) \right], \end{split}\end{split}\]

where the \(B_{2k}\) are constants known as Bernoulli numbers. Unless we happen to be fortunate enough to have a function with \(f'(b)=f'(a)\), we should expect truncation error at second order and no better.

Extrapolation¶

If evaluations of \(f\) are computationally expensive, we want to get as much accuracy as possible from them by using a higher-order formula. There are many routes for doing so; for example, we could integrate a not-a-knot cubic spline interpolant. However, splines are difficult to compute by hand, and as a result different methods were developed before computers came on the scene.

Knowing the structure of the error allows the use of extrapolation to improve accuracy. Suppose a quantity \(A_0\) is approximated by an algorithm \(A(h)\) with an error expansion

(152)¶\[ A_0 = A(h) + c_1 h + c_2 h^2 + c_3 h^3 + \cdots.\]

Crucially, it is not necessary to know the values of the error constants \(c_k\), merely that they exist and are independent of \(h\). In the case of the trapezoid formula, we have

\[ I = T_f(n) + c_2 h^2 + c_4 h^{4} + \cdots,\]

as proved by the Euler–Maclaurin formula (151). The error constants depend on \(f\) and can’t be evaluated in general, but we know that this expansion holds.

For convenience we recast the error expansion in terms of \(n=O(h^{-1})\):

(153)¶\[ I = T_f(n) + c_2 n^{-2} + c_4 n^{-4} + \cdots,\]

We make the simple observation that

(154)¶\[ I = T_f(2n) + \tfrac{1}{4} c_2 n^{-2} + \tfrac{1}{16} c_4 n^{-4} + \cdots.\]

It follows that if we combine (153) and (154) correctly, we can cancel out the second-order term in the error. Specifically, define

(155)¶\[ S_f(2n) = \frac{1}{3} \Bigl[ 4 T_f(2n) - T_f(n) \Bigr].\]

(We associate \(2n\) rather than \(n\) with the extrapolated result because of the total number of nodes needed.) Then

(156)¶\[ I = S_f(2n) + O(n^{-4}) = b_4 n^{-4} + b_6 n^{-6} + \cdots.\]

The formula (155) is called Simpson’s formula. A different presentation and derivation are considered in an exercise.

Equation (156) is another particular error expansion in the form (152), so we can extrapolate again! The details change only a little. Considering that

\[ I = S_f(4n) = \tfrac{1}{16} b_4 n^{-4} + \tfrac{1}{64} b_6 n^{-6} + \cdots,\]

the proper combination this time is

(157)¶\[ R_f(4n) = \frac{1}{15} \Bigl[ 16 S_f(4n) - S_f(2n) \Bigr],\]

which is sixth-order accurate. Clearly the process can be repeated to get eighth-order accuracy and beyond. Doing so goes by the name of Romberg integration, which we will not present in full generality.

Node doubling¶

Note in (157) that \(R_f(4n)\) depends on \(S_f(2n)\) and \(S_f(4n)\), which in turn depend on \(T_f(n)\), \(T_f(2n)\), and \(T_f(4n)\). There is a useful benefit realized by doubling of the nodes in each application of the trapezoid formula. For simplicity, suppose that \([a,b]=[0,1]\) and that \(n=2m\) for some positive integer \(m\). The nodes are

\[\Bigl[\, 0, \;\frac{1}{2m}, \;\frac{2}{2m}, \;\frac{3}{2m}, \;\frac{4}{2m}, \;\ldots \;\frac{2m-3}{2m}, \;\frac{2m-2}{2m}, \;\frac{2m-1}{2m}, \; 1 \,\Bigr].\]

Suppose we delete every other node:

\[\Bigl[\, 0, \;\frac{2}{2m}, \;\frac{4}{2m}, \;\ldots \;\frac{2m-2}{2m}, \; 1 \,\Bigr] \quad = \quad \Bigl[\, 0, \;\frac{1}{m}, \;\frac{2}{m}, \;\ldots \;\frac{m-1}{m}, \; 1 \,\Bigr] .\]

What remains are the nodes with \(n=m\). That is, if we have computed \(T_f(m)\) and want to compute \(T_f(2m)\), we begin with half of the evaluations of \(f\) already in our pocket. More specifically,

(158)¶\[\begin{split}\begin{split} T_f(2m) & = \frac{1}{2m} \left[ \tfrac{1}{2} f(0) + \tfrac{1}{2} f(1) + \sum_{k=1}^{m-1} f\Bigl( \tfrac{2k-1}{2m} \Bigr) + f\Bigl( \tfrac{2k}{2m} \Bigr) \right] \\ &= \frac{1}{2m} \left[ \tfrac{1}{2} f(0) + \tfrac{1}{2} f(1) + \sum_{k=1}^{m-1} f\Bigl( \tfrac{k}{m} \Bigr) \right] + \frac{1}{2m} \sum_{k=1}^{m} f\Bigl( \tfrac{2k-1}{2m} \Bigr) \\ &= \frac{1}{2} T_f(m) + \frac{1}{2m} \sum_{k=1}^{m-1} f\left( t_{2k-1} \right), \end{split}\end{split}\]

where the nodes referenced in the last line are relative to \(n=2m\). To summarize: when \(n\) is doubled, new integrand evaluations are needed only at the odd-numbered nodes of the finer grid. Although we derived this result in the particular interval \([0,1]\), it is valid for any interval.

Exercises¶

⌨ For each integral below, use trapezoid to estimate the integral for \(n=10\cdot 2^k\) nodes for \(k=1,2,\ldots,10\). Make a log–log plot of the errors and confirm or refute second-order accuracy. (These integrals were taken from [BJL05].)

(a) \(\displaystyle \int_0^1 x\log(1+x)\, dx = \frac{1}{4}\)

(b) \(\displaystyle \int_0^1 x^2 \tan^{-1}x\, dx = \frac{\pi-2+2\log 2}{12}\)

(c) \(\displaystyle \int_0^{\pi/2}e^x \cos x\, dx = \frac{e^{\pi/2}-1}{2}\)

(d) \(\displaystyle \int_0^1 \sqrt{x} \log(x) \, dx = -\frac{4}{9}\) (Note: Although the integrand has the limiting value zero as \(x\to 0\), it cannot be evaluated naively at \(x=0\). You can start the integral at \(x=\macheps\) instead.)

(e) \(\displaystyle \int_0^1 \sqrt{1-x^2}\,\, dx = \frac{\pi}{4}\)
✍ The Euler–Maclaurin error expansion (151) for the trapezoid formula implies that if we could cancel out the term due to \(f'(b)-f'(a)\), we would obtain fourth-order accuracy. We should not assume that \(f'\) is available, but approximating it with finite differences can achieve the same goal. Suppose the forward difference formula (138) is used for \(f'(a)\), and its reflected backward difference is used for \(f'(b)\). Show that the resulting modified trapezoid formula is

(159)¶\[ G_f(h) = T_f(h) - \frac{h}{24} \left[ 3\Bigl( f(x_n)+f(x_0) \Bigr) -4\Bigr( f(x_{n-1}) + f(x_1) \Bigr) + \Bigl( f(x_{n-2})+f(x_2) \Bigr) \right],\]

which is known as a Gregory integration formula.
⌨ Repeat each integral in exercise 1 above using Gregory integration (159) instead of the trapezoid formula. Compare the observed errors to fourth-order convergence.
✍ Simpson’s formula can be derived without appealing to extrapolation.

(a) Show that

\[p(x) = \beta + \frac{\gamma-\alpha}{2h}\, x + \frac{\alpha-2\beta+\gamma}{2h^2}\, x^2\]

interpolates the three points \((-h,\alpha)\), \((0,\beta)\), and \((h,\gamma)\).

(b) Find

\[ \int_{-h}^h p(s)\, ds,\]

where \(p\) is the quadratic polynomial from part~(a), in terms of \(h\), \(\alpha\), \(\beta\), and \(\gamma\).

(c) Assume equally spaced nodes in the form \(t_i=a+ih\), for \(h=(b-a)/n\) and \(i=0,\ldots,n\). Suppose \(f\) is approximated by \(p(x)\) over the subinterval \([t_{i-1},t_{i+1}]\). Apply the result from part (b) to find

\[ \int_{t_{i-1}}^{t_{i+1}} f(x)\, dx \approx \frac{h}{3} \bigl[ f(t_{i-1}) + 4f(t_i) + f(t_{i+1}) \bigr].\]

(Use the change of variable \(s=x-t_i\).)

(d) Now also assume that \(n=2m\) for an integer \(m\). Derive Simpson’s formula,

(160)¶\[\begin{split} \begin{split} \int_a^b f(x)\, dx \approx \frac{h}{3}\bigl[ &f(t_0) + 4f(t_1) + 2f(t_2) + 4f(t_3) + 2f(t_4) + \cdots\\ &+ 2f(t_{n-2}) + 4f(t_{n-1}) + f(t_n) \bigr]. \end{split}\end{split}\]
✍ Show that the Simpson formula (160) is equivalent to \(S_f(n/2)\), given the definition of \(S_f\) in (155).
⌨ For each integral in exercise 1 above, apply the Simpson formula (160) and compare the errors to fourth-order convergence.
⌨ For \(n=10,20,30,\ldots,200\), compute the trapezoidal approximation to

\[\int_{0}^1 \frac{1}{2.01+\sin (6\pi x)-\cos(2\pi x)} \,d x \approx 0.9300357672424684.\]

Make two separate plots of the absolute error as a function of \(n\), one using log–log scales and the other using log only for the \(y\)-axis. The graphs suggest that the error asymptotically behaves as \(C \alpha^n\) for some \(C>0\) and some \(0<\alpha<1\). How does this result relate to (151)?
⌨ For each integral in exercise 1 above, extrapolate the trapezoidal results two levels to get sixth-order accurate results, and compare the expected convergence rate to the observed errors.
✍ Find a formula like (157) that extrapolates two values of \(R_f\) to obtain an 8th-order accurate one.

1: Numerical integration also goes by the older name quadrature.
2: Some texts distinguish between a formula for a single subinterval \([t_{k-1},t_k]\) and a “composite” formula that adds them up over the whole interval to get something like our (150).

Fundamentals of Numerical Computation

Numerical integration¶

Trapezoid formula¶

Extrapolation¶

Node doubling¶

Exercises¶