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)

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, there is also more roundoff error when we add them, the sum is not exactly one.

8.6Newton-Cotes integration