5.2 Composite Newton-Cotes formulas

If the interval \([a,b]\) in the quadrature is large, then the Newton-Cotes formulas will give poor approximations. The quadrature error depends on \(h=(b-a)/n\) (closed formulas), and if \(b-a\) is large, then so is \(h\), hence error. If we raise \(n\) to compensate for large interval, then we face a problem discussed earlier: error due to the oscillatory behavior of high-degree interpolating polynomials that use equally-spaced nodes. A solution is to break up the domain into smaller intervals and use a Newton-Cotes rule with a smaller \(n\) on each subinterval: this is known as a composite rule.

Example 44

Let’s compute \(\int_{0}^{2}e^{x}\sin xdx\). The antiderivative can be computed using integration by parts, and the true value of the integral to 6 digits is 5.39689. If we apply the Simpson’s rule we get:

\[\begin{equation*} \int_{0}^{2}e^{x}\sin xdx\approx\frac{1}{3}(e^{0}\sin 0+4e\sin 1+e^{2}\sin 2)=5.28942. \end{equation*}\]

If we partition the integration domain \((0,2)\) into \((0,1)\) and \((1,2)\), and apply Simpson’s rule to each domain separately, we get

\[\begin{align*} \int_{0}^{2}e^{x}\sin xdx & =\int_{0}^{1}e^{x}\sin xdx+\int_{1}^{2}e^{x}\sin xdx\\ &\approx \frac{1}{6}(e^{0}\sin 0+4e^{0.5}\sin(0.5)+e\sin 1)+\frac{1}{6}(e\sin 1+4e^{1.5}\sin(1.5)+e^{2}\sin 2)\\ &= 5.38953, \end{align*}\]

improving the accuracy significantly. Note that we have used five nodes, \(0,0.5,1,1.5,2\), which split the domain \((0,2)\) into four subintervals.

The composite rules for midpoint, trapezoidal, and Simpson’s rule, with their error terms, are:

• Composite Midpoint rule

Let \(f\in C^{2}[a,b],\) \(n\) be even, \(h=\frac{b-a}{n+2},\) and \(x_{j}=a+(j+1)h\) for \(j=-1,0,...,n+1.\) The composite Midpoint rule for \(n+2\) subintervals is

(28)\[\begin{align} \int_{a}^{b}f(x)dx=2h\sum_{j=0}^{n/2}f(x_{2j})+\frac{b-a}{6}h^{2}f''(\xi) \end{align}\]

for some \(\xi\in(a,b)\).

• Composite Trapezoidal rule Let \(f\in C^{2}[a,b],\) \(h=\frac{b-a}{n}\), and \(x_{j}=a+jh\) for \(j=0,1,...,n.\) The composite Trapezoidal rule for \(n\) subintervals is

(29)\[\begin{align} \int_{a}^{b}f(x)dx=\frac{h}{2}\left[f(a)+2\sum_{j=1}^{n-1}f(x_{j})+f(b)\right]-\frac{b-a}{12}h^{2}f''(\xi) \end{align}\]

for some \(\xi\in(a,b)\).

• Composite Simpson’s rule Let \(f\in C^{4}[a,b],\) \(n\) be even, \(h=\frac{b-a}{n},\) and \(x_{j}=a+jh\) for \(j=0,1,...,n.\) The composite Simpson’s rule for \(n\) subintervals is

(30)\[\int_{a}^{b}f(x)dx=\frac{h}{3}\left[f(a)+2\sum_{j=1}^{\frac{n}{2}-1}f(x_{2j})+4\sum_{j=1}^{n/2}f(x_{2j-1})+f(b)\right]-\frac{b-a}{180}h^{4}f^{(4)}(\xi)\]

for some \(\xi\in(a,b)\).

Exercise 4.2-1

Show that the quadrature rule in Example 44 corresponds to taking \(n=4\) in the composite Simpson’s formula (30).

Exercise 4.2-2

Show that the absolute error for the composite trapezoidal rule decays at the rate of \(1/n^2\), and the absolute error for the composite Simpson’s rule decays at the rate of \(1/n^4\), where \(n\) is the number of subintervals.

Example 45

Determine \(n\) that ensures the composite Simpson’s rule approximates \(\int_{1}^{2}x\log xdx\) with an absolute error of at most \(10^{-6}\).

Solution.

The error term for the composite Simpson’s rule is \(\frac{b-a}{180}h^{4}f^{(4)}(\xi)\) where \(\xi\) is some number between \(a=1\) and \(b=2\), and \(h=(b-a)/n.\) Differentiate to get \(f^{(4)}(x)=\frac{2}{x^{3}}.\) Then

\[\begin{equation*} \frac{b-a}{180}h^{4}f^{(4)}(\xi)=\frac{1}{180}h^{4}\frac{2}{\xi^{3}}\leq\frac{h^{4}}{90} \end{equation*}\]

where we used the fact that \(\frac{2}{\xi^{3}}\leq\frac{2}{1}=2\) when \(\xi\in(1,2)\). Now make the upper bound less than \(10^{-6}\), that is,

\[\begin{equation*} \frac{h^{4}}{90}\leq10^{-6}\Rightarrow\frac{1}{n^{4}(90)}\leq10^{-6}\Rightarrow n^{4}\geq\frac{10^{6}}{90}\approx11111.11 \end{equation*}\]

which implies \(n\geq 10.27.\) Since \(n\) must be even for Simpson’s rule, this means the smallest value of \(n\) to ensure an error of at most \(10^{-6}\) is 12.

Using the Python code for the composite Simpson’s rule that will be introduced next, we get 0.6362945608 as the estimate, using 10 digits. The correct integral to 10 digits is \(0.6362943611\), which means an absolute error of \(2\times10^{-7}\), better than the expected \(10^{-6}\).

Python codes for Newton-Cotes formulas

We write codes for the trapezoidal and Simpson’s rules, and the composite Simpson’s rule. Coding trapezoidal and Simpson’s rule is straightforward.

Trapezoidal rule

def trap(f, a, b):
    return (f(a)+f(b))*(b-a)/2

Let’s verify the calculations of Example 42:

trap(lambda x: x**x, 0.5, 1)

0.42677669529663687

Simpson’s rule

def simpson(f, a, b):
    return (f(a)+4*f((a+b)/2)+f(b))*(b-a)/6

simpson(lambda x: x**x, 0.5, 1)

0.4109013813880978

Recall that the degree of accuracy of Simpson’s rule is 3. This means the rule integrates polynomials \(1, x, x^2,x^3\) exactly, but not \(x^4\). We can use this as a way to verify our code:

simpson(lambda x: x, 0, 1)

0.5

simpson(lambda x: x**2, 0, 1)

0.3333333333333333

simpson(lambda x: x**3, 0, 1)

0.25

simpson(lambda x: x**4, 0, 1)

0.20833333333333334

Composite Simpson’s rule

Next we code the composite Simpson’s rule, and verify the result of Example 45.

import numpy as np

def compsimpson(f, a, b, n):
    h = (b-a)/n
    nodes = np.zeros(n+1)
    for i in range(n+1):
        nodes[i] = a+i*h
    sum = f(a)+f(b)
    for i in range(2, n-1, 2):
        sum += 2*f(nodes[i])
    for i in range(1, n, 2):
        sum += 4*f(nodes[i])
    return sum*h/3

compsimpson(lambda x: x*np.log(x), 1, 2, 12)

0.636294560831306

Exercise 4.2-3

Determine the value of \(n\) required to approximate

\[\begin{equation*} \int_{0}^{2}\frac{1}{x+1}dx \end{equation*}\]

to within \(10^{-4}\), and compute the approximation, using the composite trapezoidal and composite Simpson’s rule.

Composite rules and roundoff error

As we increase \(n\) in the composite rules to lower error, the number of function evaluations increases, and a natural question to ask would be whether roundoff error could accumulate and cause problems. Somewhat remarkably, the answer is no. Let’s assume the roundoff error associated with computing \(f(x)\) is bounded for all \(x,\) by some positive constant \(\epsilon\). And let’s try to compute the roundoff error in composite Simpson rule. Since each function evaluation in the composite rule incorporates an error of (at most) \(\epsilon\), the total error is bounded by

\[\begin{align*} \frac{h}{3}\left[\epsilon+2\sum_{j=1}^{\frac{n}{2}-1}\epsilon+4\sum_{j=1}^{n/2}\epsilon+\epsilon\right] & \leq\frac{h}{3}\left[\epsilon+2\left(\frac{n}{2}-1\right)\epsilon+4\left(\frac{n}{2}\right)\epsilon+\epsilon\right] =\frac{h}{3}(3n\epsilon)=hn\epsilon. \end{align*}\]

However, since \(h=(b-a)/n\), the bound simplifies as \(hn\epsilon=(b-a)\epsilon.\) Therefore no matter how large \(n\) is, that is, how large the number of function evaluations is, the roundoff error is bounded by the same constant \((b-a)\epsilon\) which only depends on the size of the interval.

Exercise 4.2-4

(This problem shows that numerical quadrature is stable with respect to error in function values.) Assume the function values \(f(x_{i})\) are approximated by \(\tilde{f}(x_{i})\), so that \(|f(x_{i})-\tilde{f}(x_{i})|<\epsilon\) for any \(x_{i}\in(a,b).\) Find an upper bound on the error of numerical quadrature \(\sum w_{i}f(x_{i})\) when it is actually computed as \(\sum w_{i}\tilde{f}(x_{i})\).