3: Integration

We want to construct numerical algorithms that can perform definite integrals of the form

\[I=\int_{a}^{b} f(x) d x . \nonumber \]

Calculating these definite integrals numerically is called numerical integration, numerical quadrature, or more simply quadrature.

Elementary formulas

We first consider integration from 0 to \(h\), with \(h\) small, to serve as the building blocks for integration over larger domains. We here define \(I_{h}\) as the following integral:

\[I_{h}=\int_{0}^{h} f(x) d x . \nonumber \]

To perform this integral, we consider a Taylor series expansion of \(f(x)\) about the value \(x=h / 2\) :

\[\begin{aligned} f(x)=f(h / 2)+(x-h / 2) f^{\prime}(h / 2) &+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ +& \frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots \end{aligned} \nonumber \]

Midpoint rule

The midpoint rule makes use of only the first term in the Taylor series expansion. Here, we will determine the error in this approximation. Integrating,

\[\begin{aligned} I_{h}=h f(h / 2)+\int_{0}^{h}((x-&h / 2) f^{\prime}(h / 2)+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ &\left.+\frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d x . \end{aligned} \nonumber \]

Changing variables by letting \(y=x-h / 2\) and \(d y=d x\), and simplifying the integral depending on whether the integrand is even or odd, we have

\[\begin{aligned} I_{h}=& h f(h / 2) \\ &+\int_{-h / 2}^{h / 2}\left(y f^{\prime}(h / 2)+\frac{y^{2}}{2} f^{\prime \prime}(h / 2)+\frac{y^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{y^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y \\ =& h f(h / 2)+\int_{0}^{h / 2}\left(y^{2} f^{\prime \prime}(h / 2)+\frac{y^{4}}{12} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y . \end{aligned} \nonumber \]

The integrals that we need here are

\[\int_{0}^{\frac{h}{2}} y^{2} d y=\frac{h^{3}}{24}, \quad \int_{0}^{\frac{h}{2}} y^{4} d y=\frac{h^{5}}{160} . \nonumber \]

Therefore,

\[I_{h}=h f(h / 2)+\frac{h^{3}}{24} f^{\prime \prime}(h / 2)+\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber \]

Trapezoidal rule

From the Taylor series expansion of \(f(x)\) about \(x=h / 2\), we have

\[f(0)=f(h / 2)-\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)-\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber \]

and

\[f(h)=f(h / 2)+\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)+\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber \]

Adding and multiplying by \(h / 2\) we obtain

\[\frac{h}{2}(f(0)+f(h))=h f(h / 2)+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber \]

We now substitute for the first term on the right-hand-side using the midpoint rule formula:

\[\begin{aligned} \frac{h}{2}(f(0)+f(h))=\left(I_{h}-\frac{h^{3}}{24} f^{\prime \prime}(h / 2)-\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)\right) &\, \\ &+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \end{aligned} \nonumber \]

and solving for \(I_{h}\), we find

\[I_{h}=\frac{h}{2}(f(0)+f(h))-\frac{h^{3}}{12} f^{\prime \prime}(h / 2)-\frac{h^{5}}{480} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber \]

Simpson’s rule

To obtain Simpson’s rule, we combine the midpoint and trapezoidal rule to eliminate the error term proportional to \(h^{3}\). Multiplying (3.4) by two and adding to (3.8), we obtain

\[3 I_{h}=h\left(2 f(h / 2)+\frac{1}{2}(f(0)+f(h))\right)+h^{5}\left(\frac{2}{1920}-\frac{1}{480}\right) f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber \]

or

\[I_{h}=\frac{h}{6}(f(0)+4 f(h / 2)+f(h))-\frac{h^{5}}{2880} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber \]

Usually, Simpson’s rule is written by considering the three consecutive points \(0, h\) and \(2 h\). Substituting \(h \rightarrow 2 h\), we obtain the standard result

\[I_{2 h}=\frac{h}{3}(f(0)+4 f(h)+f(2 h))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(h)+\ldots \nonumber \]

Composite rules

We now use our elementary formulas obtained for (3.2) to perform the integral given by (3.1).

Trapezoidal rule

We suppose that the function \(f(x)\) is known at the \(n+1\) points labeled as \(x_{0}, x_{1}, \ldots, x_{n}\), with the endpoints given by \(x_{0}=a\) and \(x_{n}=b\). Define

\[f_{i}=f\left(x_{i}\right), \quad h_{i}=x_{i+1}-x_{i} . \nonumber \]

Then the integral of (3.1) may be decomposed as

\[\begin{aligned} \int_{a}^{b} f(x) d x &=\sum_{i=0}^{n-1} \int_{x_{i}}^{x_{i+1}} f(x) d x \\ &=\sum_{i=0}^{n-1} \int_{0}^{h_{i}} f\left(x_{i}+s\right) d s, \end{aligned} \nonumber \]

where the last equality arises from the change-of-variables \(s=x-x_{i}\). Applying the trapezoidal rule to the integral, we have

\[\int_{a}^{b} f(x) d x=\frac{1}{2} \sum_{i=0}^{n-1} h_{i}\left(f_{i}+f_{i+1}\right) . \nonumber \]

If the points are not evenly spaced, say because the data are experimental values, then the \(h_{i}\) may differ for each value of \(i\) and (3.13) is to be used directly.

However, if the points are evenly spaced, say because \(f(x)\) can be computed, we have \(h_{i}=h\), independent of \(i\). We can then define

\[x_{i}=a+i h, \quad i=0,1, \ldots, n ; \nonumber \]

and since the end point \(b\) satisfies \(b=a+n h\), we have

\[h=\frac{b-a}{n} . \nonumber \]

The composite trapezoidal rule for evenly space points then becomes

\[\begin{align} \nonumber \int_{a}^{b} f(x) d x &=\frac{h}{2} \sum_{i=0}^{n-1}\left(f_{i}+f_{i+1}\right) \\ &=\frac{h}{2}\left(f_{0}+2 f_{1}+\cdots+2 f_{n-1}+f_{n}\right) . \end{align} \nonumber \]

The first and last terms have a multiple of one; all other terms have a multiple of two; and the entire sum is multiplied by \(h / 2\)

Figure 3.1: Adaptive Simpson quadrature: Level \(1 .\)

Simpson’s rule

We here consider the composite Simpson’s rule for evenly space points. We apply Simpson’s rule over intervals of \(2 h\), starting from \(a\) and ending at \(b\) :

\[\begin{aligned} \int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+f_{2}\right)+\frac{h}{3}\left(f_{2}+4 f_{3}+f_{4}\right)+\ldots & \\ &+\frac{h}{3}\left(f_{n-2}+4 f_{n-1}+f_{n}\right) . \end{aligned} \nonumber \]

Note that \(n\) must be even for this scheme to work. Combining terms, we have

\[\int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+2 f_{2}+4 f_{3}+2 f_{4}+\cdots+4 f_{n-1}+f_{n}\right) . \nonumber \]

The first and last terms have a multiple of one; the even indexed terms have a multiple of 2; the odd indexed terms have a multiple of 4 ; and the entire sum is multiplied by \(h / 3\).

Adaptive integration

The useful MATLAB function quad.m performs numerical integration using adaptive Simpson quadrature. The idea is to let the computation itself decide on the grid size required to achieve a certain level of accuracy. Moreover, the grid size need not be the same over the entire region of integration.

We begin the adaptive integration at what is called Level 1. The uniformly spaced points at which the function \(f(x)\) is to be evaluated are shown in Fig. 3.1. The distance between the points \(a\) and \(b\) is taken to be \(2 h\), so that

\[h=\frac{b-a}{2} . \nonumber \]

Integration using Simpson’s rule (3.11) with grid size \(h\) yields for the integral \(I\),

\[I=\frac{h}{3}(f(a)+4 f(c)+f(b))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \nonumber \]

where \(\xi\) is some value satisfying \(a \leq \xi \leq b\). Integration using Simpson’s rule twice with grid size \(h / 2\) yields

\[I=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b))-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{l}\right)-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{r}\right), \nonumber \]

with \(\xi_{l}\) and \(\xi_{r}\) some values satisfying \(a \leq \xi_{l} \leq c\) and \(c \leq \xi_{r} \leq b\).

We now define the two approximations to the integral by

\[\begin{aligned} &S_{1}=\frac{h}{3}(f(a)+4 f(c)+f(b)), \\ &S_{2}=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b)), \end{aligned} \nonumber \]

and the two associated errors by

\[\begin{aligned} &E_{1}=-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \\ &E_{2}=-\frac{h^{5}}{2^{5} \cdot 90}\left(f^{\prime \prime \prime \prime}\left(\xi_{l}\right)+f^{\prime \prime \prime \prime}\left(\xi_{r}\right)\right) . \end{aligned} \nonumber \]

We now ask whether the value of \(S_{2}\) for the integral is accurate enough, or must we further refine the calculation and go to Level 2? To answer this question, we make the simplifying approximation that all of the fourth-order derivatives of \(f(x)\) in the error terms are equal; that is,

\[f^{\prime \prime \prime \prime}(\xi)=f^{\prime \prime \prime \prime}\left(\xi_{l}\right)=f^{\prime \prime \prime \prime}\left(\xi_{r}\right)=C . \nonumber \]

Then

\[\begin{aligned} &E_{1}=-\frac{h^{5}}{90} C, \\ &E_{2}=-\frac{h^{5}}{2^{4} \cdot 90} C=\frac{1}{16} E_{1} . \end{aligned} \nonumber \]

Now since the integral is equal to the approximation plus its associated error,

\[S_{1}+E_{1}=S_{2}+E_{2} \text {, } \nonumber \]

and since

\[E_{1}=16 E_{2} \text {, } \nonumber \]

we can derive an estimate for the error term \(E_{2}\) :

\[E_{2}=\frac{1}{15}\left(S_{2}-S_{1}\right) \text {. } \nonumber \]

Therefore, given some specific value of the tolerance tol, if

\[\left|\frac{1}{15}\left(S_{2}-S_{1}\right)\right|<\text { tol, } \nonumber \]

then we can accept \(S_{2}\) as \(I\). If the error estimate is larger in magnitude than tol, then we proceed to Level 2 . The computation at Level 2 further divides the integration interval from \(a\) to \(b\) into the two integration intervals \(a\) to \(c\) and \(c\) to \(b\), and proceeds with the above procedure independently on both halves. Integration can be stopped on either half provided the tolerance is less than tol/2 (since the sum of both errors must be less than tol). Otherwise, either half can proceed to Level 3 , and so on.

As a side note, the two values of \(I\) given above (for integration with step size \(h\) and \(h / 2)\) can be combined to give a more accurate value for I given by

\[I=\frac{16 S_{2}-S_{1}}{15}+\mathrm{O}\left(h^{7}\right), \nonumber \]

where the error terms of \(\mathrm{O}\left(h^{5}\right)\) approximately cancel. This free lunch, so to speak, is called Richardson’s extrapolation.