8.14Adaptive quadrature

#%config InlineBackend.figure_format = 'svg'
from pylab import *

The error in some numerical quadrature depends on the number of nodes used which the user has to select. It is desirable to have an automatic algorithm where the user specifies the error ε that is acceptable and the algorithm should select the number of nodes $n$ so that

The error clearly depends on the smoothness of the integrand $f$ also. If $f$ has a singularity, e.g., $f(x) = \sqrt{x}$ has singular derivatives at $x=0$, then the error will be large. To have an acceptable error we must use a large $n$ (small $h$). But this will increase the number of nodes. Actually, most of the error will come from the region near the singularity and we can use a small $h$ only in those regions, and not everywhere. This requires the ability to use non-uniformly spaced nodes in the numerical integration formula, and all the composite rules allow us to do this.

8.14.1Adaptive Simpson quadrature¶

Let $n$ be even and we use $n+1$ points. The composite Simpson quadrature is

with error formula

Note that we have a decomposition of the total error into error contributions from each of the sub-intervals. If $f^{(4)}$ is large in some sub-interval, we should choose the length of that sub-interval $x_{2j} - x_{2j-2}$ to be small, in order to reduce the error.

The best strategy is to equi-distribute the error among all the sub-intervals.

Define

then

We want

This can be achieved if

If this condition is not satisfied, then we iteratively refine the sub-intervals until desired error level is reached.

8.14.1.1Algorithm I¶

We will assume that we can estimate $M_j$ and hence we can estimate the error. The algorithm will iteratively increase the number of sub-intervals until the desired error level is reached.

1. Start with three nodes, i.e., one Simpson rule, and compute $I_2(f)$ and estimate the error. If $|I(f) - I_2(f)| \le \epsilon$, then we stop and the answer is $I_2(f)$.

2. Otherwise, insert a point mid-way between the three points and use two Simpson rules. Estimate the error from each rule. If both errors are $\le \epsilon/2$, then stop.

3. If not, suppose first error $> \epsilon/2$ and second error $\le \epsilon/2$. Refine the first Simpson rule by inserting two points.

4. Now the error from the first two rules must be $\le \epsilon/4$. If this is true, then we stop.

5. Otherwise, continue with sub-division.

8.14.1.2Algorithm II¶

We may not be able to estimate the derivatives, and the only information available is the function values. We can approximate the derivatives by first constructing a local approximation and differentiating it. But we have seen other ways to estimate the error in integral, e.g., using Richardson extrapolation.

Set $\alpha=a$ and $\beta=b$ and define

The Simpson rule is

To estimate the error, we compute a refined estimate as

and use error estimate based on Richarson extrapolation, see Section 8.13.2.1

as the error estimate.

1. If $E_{\alpha,\beta} \le \epsilon$ then keep interval $[\alpha,\beta]$ and use $I_{\alpha,\beta}^{(2)}$ as estimate of $I_{\alpha,\beta}$.

2. Otherwise, divide $[\alpha,\beta]$ into $[\alpha,\gamma]$ and $[\gamma,\beta]$. Now we want to satisfy

   $$E_{\alpha,\gamma} \le \frac{\epsilon}{2}, \qquad E_{\gamma,\beta} \le \frac{\epsilon}{2}$$

   (12)

3. If both conditions are satisfied, then the integral is $I_{\alpha,\gamma}^{(2)} + I_{\gamma,\beta}^{(2)}$.

4. Otherwise, suppose

   $$E_{\alpha,\gamma} > \frac{\epsilon}{2}, \qquad E_{\gamma,\beta} \le \frac{\epsilon}{2}$$

   (13)

   Then we keep $I_{\gamma,\beta}^{(2)}$ and divide $[\alpha,\gamma]$ into $[\alpha,\delta]$ and $[\delta,\gamma]$ with $\delta=\half(\alpha + \gamma)$. Now try to satisfy

   $$E_{\alpha,\delta} \le \frac{\epsilon}{4}, \qquad E_{\delta,\gamma} \le \frac{\epsilon}{4}$$

   (14)

   If both conditions are satisfied, then integral is $I_{\alpha,\delta}^{(2)} + I_{\delta,\gamma}^{(2)} + I_{\gamma,\beta}^{(2)}$.

5. Continue this process until the errors in all sub-intervals is less than the tolerance.

This algorithm is implemented in the following code.

• simpson1 performs Simpson quadrature in one interval $[a,b]$.
• asimpson performs adaptive quadrature in a recursive manner.

from scipy.integrate import simpson

# Apply Simpson using points {a, (a+b)/2, c}
def simpson1(f,a,b):
    h = 0.5*(b-a)
    c = 0.5*(a+b)
    return (h/3.0)*(f(a) + 4*f(c) + f(b))

# Performs adaptive Simpson to make error < eps
def asimpson(f,a,b,eps):
    c = 0.5*(a+b)
    i1 = simpson1(f,a,b)
    i2 = simpson1(f,a,c) + simpson1(f,c,b)
    err = (1.0/15.0) * abs(i2 - i1)
    if err < eps:
        print('[%12.6g %12.6g] %12.6g %12.6g %12.6g' % (a,b,err,eps,i2))
        return i2
    else:
        eps1 = 0.5*eps
        return asimpson(f,a,c,eps1) + asimpson(f,c,b,eps1)