6.4: 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. 6.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 (6.5) with grid size $h$ yields

$$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

$$\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)), \\ &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$$

Now we ask whether $S_{2}$ 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$$

Then since

$$S_{1}+E_{1}=S_{2}+E_{2} \nonumber$$

and

$$E_{1}=16 E_{2} \nonumber$$

we have for our estimate for the error term $E_{2}$,

$$E_{2}=\frac{1}{15}\left(S_{2}-S_{1}\right) \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 tolerance is not achieved for $I$, 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.