Newton-Cotes integration - Numerical Analysis

from pylab import *
from scipy.integrate import newton_cotes

The Trapezoidal and Simpson rules were based on integrating a polynomial approximation. Such methods are called Newton-Cotes integration methods. If we choose $n+1$ distinct points in $[a,b]$ and approximate the function $f : [a,b] \to \re$ by interpolation

p_n(x) = \sum_{j=0}^n f(x_{j,n}) \ell_{j,n}(x), \qquad x \in [a,b]

(1)

the integral can be approximated by

I_n(f) = \int_a^b p_n(x) \ud x = \sum_{j=0}^n f(x_{j,n}) \clr{red}{ \int_a^b \ell_{j,n}(x) \ud x } = \sum_{j=0}^n \clr{red}{w_{j,n}} f(x_{j,n})

(2)

where the weights are given by

w_{j,n} = \int_a^b \ell_{j,n}(x) \ud x

(3)

As an example for $n=3$ and using uniformly spaced points, we interpolate by a cubic polynomial and the integral is approximated as

I_3(f) = \frac{3h}{8}[f_0 + 3 f_1 + 3 f_2 + f_3]

(4)

which is known as Boole’s rule. For general $n$ , we can state the following error estimates.

The accuracy of a quadrature formula can be characterized by the largest set of polynomials which it can integrate exactly.

Example 1

On uniformly spaced nodes, the interpolating polynomial $p_n(x)$ does not always converge. Hence the Newton-Cotes formulae may not converge also. Consider

I = \int_{-4}^4 \frac{\ud x}{1 + x^2} = 2 \tan^{-1}(4) \approx 2.6516

(9)

xmin, xmax = -4.0, 4.0
f = lambda x: 1/(1 + x**2)
exact = 2.0 * arctan(4.0)

for n in arange(2,12,2):
    w,B = newton_cotes(n, equal=1)
    x = linspace(xmin,xmax,n+1)
    dx = (xmax-xmin)/n
    quad = dx * sum(w * f(x))
    error = abs(quad - exact)
    print('{:4d}  {:12.9f}  {:12.5e}'.format(n, quad, error))

   2   5.490196078   2.83856e+00
   4   2.277647059   3.73988e-01
   6   3.328798127   6.77163e-01
   8   1.941094304   7.10541e-01
  10   3.595560400   9.43925e-01

The error shown in last column is increasing.

Hence it is not advisable to use Newton-Cotes rule for large degree $n$ . However, we can use moderate but fixed degree of the interpolating polynomial and build a composite rule. The convergence here is obtained as the length of the sub-intervals tends to zero but the polynomial degree in each sub-interval remains fixed.

We now look at necessary and sufficient conditions for a numerical integration formula to converge.

Note that for an integration formula

\sum_{j=0}^n w_{j,n} = b - a

(15)

Proof 2

(a) Condition (14) implies (1) is obvious since $\mathcal{F}$ is a subset of $\cts[a,b]$ . That it implies (2) also is an application of the uniform boundedness principle.

(b) Now we prove the reverse implication. We have to use the density of $\mathcal{F}$ and the given two conditions. Let $\epsilon > 0$ be arbitrary and decompose the integration error into three parts;

\begin{align*} I(f) - I_n(f) &= \clr{red}{I(f) - I(f_\epsilon)} + \clr{blue}{I(f_\epsilon) - I_n(f_\epsilon) } + \clr{magenta}{ I_n(f_\epsilon) - I_n(f) } \\ |I(f) - I_n(f)| &\le |I(f) - I(f_\epsilon)| + |I(f_\epsilon) - I_n(f_\epsilon)| + |I_n(f_\epsilon) - I_n(f)| \\ &\le (b-a)\norm{f - f_\epsilon}_\infty + |I(f_\epsilon) - I_n(f_\epsilon)| + \norm{f - f_\epsilon}_\infty \underbrace{\sum_{j=0}^n |w_{j,n}|}_{\le B} \end{align*}

(16)

Using density, choose $f_\epsilon$ such that

\norm{f - f_\epsilon}_\infty \le \frac{\epsilon}{2(b-a + B)}

(17)

Using property (1), we can choose $n$ large enough so that

|I(f_\epsilon) - I_n(f_\epsilon)| \le \frac{\epsilon}{2}

(18)

Then

\begin{align} |I(f) - I_n(f)| &\le \frac{(b-a) \epsilon}{2(b - a + B)} + \frac{\epsilon}{2} + \frac{B \epsilon}{2 (b - a + B)} \\ &= \epsilon \end{align}

(19)

Since ε was arbitrary, we conclude that $I_n(f) \to I(f)$ .

Example 3

For the function $f(x) = 1/(1+x^2)$ , $x \in [-4,4]$ , the Newton-Cotes on uniform points does not converge.

Take $\mathcal{F}$ to be set of all polynomials. Then given any $f \in \mathcal{F}$ , clearly $I_n(f) = I(f)$ when $n \ge$ degree of $f$ , so condition (1) of theorem is satisfied.

However, property (2) does not hold for uniform points. The weights

w_{j,n} = \int_a^b \ell_{j,n}(x) \ud x

(20)

increase with $n$ taking both positive and negative values, so that $B = \infty$ .

Here are the weights for $n=20$ intervals in $[0,1]$ .

n = 20
w,B = newton_cotes(n, equal=1)
w = w / n
for w1 in w:
    print("{:18.12f}".format(w1))
print("sum(w)   = ", sum(w))
print("sum(|w|) = ", sum(abs(w)))

    0.011825273249
    0.114137717645
   -0.236478370511
    1.206186893482
   -3.771031726716
   10.336798219898
  -22.708815844147
   41.828057422070
  -64.075279490324
   82.797283471795
  -90.005367135629
   82.797283472028
  -64.075279489974
   41.828057422070
  -22.708815843856
   10.336798219942
   -3.771031726694
    1.206186893480
   -0.236478370511
    0.114137717645
    0.011825273249
sum(w)   =  0.9999999981875201
sum(|w|) =  544.1771559949144

Because the weights are large in size and take both signs, there is also more roundoff error when we add them, and the sum is not exactly one. This error increases as we use more points.

References¶

Atkinson, K. E. (2004). An Introduction to Numerical Analysis (2nd ed.). Wiley.

Numerical Analysis

Peano kernel error formulas

Numerical Analysis

Mid-point rule

8.6Newton-Cotes integration