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.