5.1 Newton-Cotes formulas

The idea is to construct the polynomial interpolant \(P(x)\) and compute \(\int_{a}^{b}P(x)dx\) as an approximation to \(\int_{a}^{b}f(x)dx\). Given nodes \(x_{0},x_{1},...,x_{n}\), the Lagrange form of the interpolant is

\[\begin{equation*} P_{n}(x)=\sum_{i=0}^{n}f(x_{i})l_{i}(x) \end{equation*}\]

and from the interpolation error formula Theorem 22, we have

\[\begin{equation*} f(x)=P_{n}(x)+(x-x_{0})\cdots(x-x_{n})\frac{f^{(n+1)}(\xi(x))}{(n+1)!}, \end{equation*}\]

where \(\xi(x)\in[a,b]\). (We have written \(\xi(x)\) instead of \(\xi\) to emphasize that \(\xi\) depends on the value of \(x\).)

Taking the integral of both sides yields

(27)\[\int_{a}^{b}f(x)dx = \underbrace{\int_{a}^{b}P_{n}(x)dx}_{\text{quadrature rule}}+\underbrace{\frac{1}{(n+1)!}\int_{a}^{b}\prod_{i=0}^{n}(x-x_{i})f^{(n+1)}(\xi(x))dx.}_{\text{error term}} % & =\sum_{i=0}^{n}\left(\underbrace{\int_{a}^{b}l_{i}(x)dx}_{w_i}\right) f(x_{i})+\text{"error"}\nonumber\]

The first integral on the right-hand side gives the quadrature rule:

\[\begin{equation*} \int_{a}^{b}P_{n}(x)dx=\int_a^b \left( \sum_{i=0}^{n}f(x_{i})l_{i}(x) \right) dx =\sum_{i=0}^{n}\left(\underbrace{\int_{a}^{b}l_{i}(x)dx}_{w_i}\right) f(x_{i})=\sum_{i=0}^n w_i f(x_i), \end{equation*}\]

and the second integral gives the error term.

We obtain different quadrature rules by taking different nodes, or number of nodes. The following result is useful in the theoretical analysis of Newton-Cotes formulas.

Theorem 28 (Weighted mean value theorem for integrals)

Suppose \(f\in C^0[a,b],\) the Riemann integral of \(g(x)\) exists on \([a,b],\) and \(g(x)\) does not change sign on \([a,b].\) Then there exists \(\xi\in(a,b)\) with \(\int_{a}^{b}f(x)g(x)dx=f(\xi)\int_{a}^{b}g(x)dx.\)

Two well-known numerical quadrature rules, trapezoidal rule and Simpson’s rule, are examples of Newton-Cotes formulas:

• Trapezoidal rule

Let \(f\in C^2[a,b]\). Take two nodes, \(x_{0}=a, x_{1}=b,\) and use the linear Lagrange polynomial

\[\begin{equation*} P_{1}(x)=\frac{x-x_{1}}{x_{0}-x_{1}}f(x_{0})+\frac{x-x_{0}}{x_{1}-x_{0}}f(x_{1}) \end{equation*}\]

to estimate \(f(x)\). Substitute \(n=1\) in Equation (27) to get

\[\begin{equation*} \int_{a}^{b}f(x)dx=\int_{a}^{b}P_{1}(x)dx+\frac{1}{2}\int_{a}^{b}\prod_{i=0}^{1}(x-x_{i})f''(\xi(x))dx, \end{equation*}\]

and then substitute for \(P_1\) to obtain

\[\begin{align*} \int_{a}^{b}f(x)dx & =\int_{a}^{b}\frac{x-x_{1}}{x_{0}-x_{1}}f(x_{0})dx+\int_{a}^{b}\frac{x-x_{0}}{x_{1}-x_{0}}f(x_{1})dx +\frac{1}{2}\int_{a}^{b}(x-x_{0})(x-x_{1})f''(\xi(x))dx. \end{align*}\]

The first two integrals on the right-hand side can be evaluated easily: the first integral is \(\frac{h}{2}f(x_{0})\) and the second one is \(\frac{h}{2}f(x_{1})\), where \(h=x_1-x_0=b-a\). Let’s evaluate

\[\begin{equation*} \int_{a}^{b}(x-x_{0})(x-x_{1})f''(\xi(x))dx. \end{equation*}\]

We will use Theorem 28 for this computation. Note that the function \((x-x_0)(x-x_1)=(x-a)(x-b)\) does not change sign on the interval \([a,b]\) and it is integrable: so this function serves the role of \(g(x)\) in Theorem 28. The other term, \(f''(\xi(x))\), serves the role of \(f(x)\). Applying the theorem, we get

\[\begin{equation*} \int_{a}^{b}(x-a)(x-b)f''(\xi(x))dx=f''(\xi) \int_a^b (x-a)(x-b)dx \end{equation*}\]

where we kept the same notation \(\xi\), somewhat inappropriately, as we moved \(f''(\xi(x))\) from inside to the outside of the integral. Finally, observe that

\[\begin{equation*} \int_a^b (x-a)(x-b)dx=\frac{(a-b)^3}{6}=\frac{-h^3}{6}, \end{equation*}\]

where \(h=b-a\). Putting all the pieces together, we have obtained

\[\begin{equation*} \int_{a}^{b}f(x)dx=\frac{h}{2}\left[f(x_{0})+f(x_{1})\right]-\frac{h^{3}}{12}f''(\xi). \end{equation*}\]

• Simpson’s rule

Let \(f\in C^4[a,b]\). Take three equally-spaced nodes, \(x_{0}=a\), \(x_{1}=a+h\), \(x_{2}=b,\) where \(h=\frac{b-a}{2},\) and use the second degree Lagrange polynomial

\[\begin{equation*} P_{2}(x)=\frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}f(x_{0})+\frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}f(x_{1})+\frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}f(x_{2}) \end{equation*}\]

to estimate \(f(x)\). Substitute \(n=2\) in Equation (27) to get

\[\begin{equation*} \int_{a}^{b}f(x)dx=\int_{a}^{b}P_{2}(x)dx+\frac{1}{3!}\int_{a}^{b}\prod_{i=0}^{2}(x-x_{i})f^{(3)}(\xi(x))dx, \end{equation*}\]

and then substitute for \(P_2\) to obtain

\[\begin{align*} \int_{a}^{b}f(x)dx & =\int_{a}^{b}\frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}f(x_{0})dx+\int_{a}^{b}\frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}f(x_{1})dx+\\ & +\int_{a}^{b}\frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}f(x_{2})dx+\frac{1}{6}\int_{a}^{b}(x-x_{0})(x-x_{1})(x-x_{2})f^{(3)}(\xi(x))dx. \end{align*}\]

The sum of the first three integrals on the right-hand side simplify as: \(\frac{h}{3}\left[f(x_{0})+4f(x_{1})+f(x_{2})\right]\). The last integral cannot be evaluated using Theorem 28 directly, like in the trapezoidal rule, since the function \((x-x_0)(x-x_1)(x-x_2)\) changes sign on \([a,b]\). However, a clever application of integration by parts transforms the integral to an integral where Theorem 28 is applicable (see Atkinson [1989] for details), and the integral simplifies as \(-\frac{h^{5}}{90}f^{(4)}(\xi)\) for some \(\xi \in (a,b)\). In summary, we obtain

\[\begin{equation*} \int_{a}^{b}f(x)dx=\frac{h}{3}\left[f(x_{0})+4f(x_{1})+f(x_{2})\right]-\frac{h^{5}}{90}f^{(4)}(\xi) \end{equation*}\]

where \(\xi\in(a,b)\).

Exercise 5.1-1

Prove that the sum of the weights in Newton-Cotes rules is \(b-a\), for any \(n\).

Definition 14

The degree of accuracy, or precision, of a quadrature formula is the largest positive integer \(n\) such that the formula is exact for \(f(x)=x^{k},\) when \(k=0,1,...,n,\) or equivalently, for any polynomial of degree less than or equal to \(n\).

Observe that the trapezoidal and Simpson’s rules have degrees of accuracy of one and three. These two rules are examples of closed Newton-Cotes formulas; closed refers to the fact that the end points \(a,b\) of the interval are used as nodes in the quadrature rule. Here is the general definition.

Definition 15 (Closed Newton-Cotes)

The \((n+1)\)-point closed Newton-Cotes formula uses nodes \(x_{i}=x_{0}+ih\), for \(i=0,1,...,n,\) where \(x_{0}=a,x_{n}=b,h=\frac{b-a}{n},\) and

\[\begin{equation*} w_{i}=\int_{x_{0}}^{x_{n}}l_{i}(x)dx=\int_{x_{0}}^{x_{n}}\prod_{j=0,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}dx. \end{equation*}\]

The following theorem provides an error formula for the closed Newton-Cotes formula. A proof can be found in Isaacson and Keller [1966].

Theorem 29

For the \((n+1)\)-point closed Newton-Cotes formula, we have:

• if \(n\) is even and \(f\in C^{n+2}[a,b]\)

\[\begin{equation*} \int_{a}^{b}f(x)dx=\sum_{i=0}^{n}w_{i}f(x_{i})+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_{0}^{n}t^{2}(t-1)\cdots(t-n)dt, \end{equation*}\]

• if \(n\) is odd and \(f\in C^{n+1}[a,b]\)

\[\begin{equation*} \int_{a}^{b}f(x)dx=\sum_{i=0}^{n}w_{i}f(x_{i})+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_{0}^{n}t(t-1)\cdots(t-n)dt \end{equation*}\]

where \(\xi\) is some number in \((a,b).\)

Some well known examples of closed Newton-Cotes formulas are the trapezoidal rule (\(n=1)\), Simpson’s rule (\(n=2)\), and Simpson’s three-eighth rule (\(n=3)\). Observe that in the \((n+1)\)-point closed Newton-Cotes formula, if \(n\) is even, then the degree of accuracy is \((n+1\)), although the interpolating polynomial is of degree \(n\). The open Newton-Cotes formulas exclude the end points of the interval.

Definition 16 (Open Newton-Cotes)

The \((n+1)\)-point open Newton-Cotes formula uses nodes \(x_{i}=x_{0}+ih\), for \(i=0,1,...,n,\) where \(x_{0}=a+h,x_{n}=b-h,h=\frac{b-a}{n+2}\) and

\[\begin{equation*} w_{i}=\int_{a}^{b}l_{i}(x)dx=\int_{a}^{b}\prod_{j=0,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}dx. \end{equation*}\]

We put \(a=x_{-1}\) and \(b=x_{n+1}\).

The error formula for the open Newton-Cotes formula is given next; for a proof see Isaacson and Isaacson and Keller [1966].

Theorem 30

For the \((n+1)\)-point open Newton-Cotes formula, we have:

• if \(n\) is even and \(f\in C^{n+2}[a,b]\)

\[\begin{equation*} \int_{a}^{b}f(x)dx=\sum_{i=0}^{n}w_{i}f(x_{i})+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_{-1}^{n+1}t^{2}(t-1)\cdots(t-n)dt, \end{equation*}\]

• if \(n\) is odd and \(f\in C^{n+1}[a,b]\)

\[\begin{equation*} \int_{a}^{b}f(x)dx=\sum_{i=0}^{n}w_{i}f(x_{i})+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_{-1}^{n+1}t(t-1)\cdots(t-n)dt, \end{equation*}\]

where \(\xi\) is some number in \((a,b)\).

The well-known midpoint rule is an example of open Newton-Cotes formula:

• Midpoint rule

Take one node, \(x_0=a+h\), which corresponds to \(n=0\) in the above theorem to obtain

\[\begin{equation*} \int_{a}^{b}f(x)dx=2hf(x_{0})+\frac{h^{3}f''(\xi)}{3} \end{equation*}\]

where \(h=(b-a)/2.\) This rule interpolates \(f\) by a constant (the values of \(f\) at the midpoint), that is, a polynomial of degree 0, but it has degree of accuracy 1.

Remark 9

Both closed and open Newton-Cotes formulas using odd number of nodes (\(n\) is even), gain an extra degree of accuracy beyond that of the polynomial interpolant on which they are based. This is due to cancellation of positive and negative error.

There are some drawbacks of Newton-Cotes formulas:

• In general, these rules are not of the highest degree of accuracy possible for the number of nodes used.

• The use of large number of equally spaced nodes may incur the erratic behavior associated with high-degree polynomial interpolation. Weights for a high-order rule may be negative, potentially leading to loss of significance errors.

• Let \(I_{n}\) denote the Newton-Cotes estimate of an integral based on \(n\) nodes. \(I_{n}\) may not converge to the true integral as \(n\rightarrow\infty\) for perfectly well-behaved integrands.

Example 42

Estimate \(\int_{0.5}^1 x^x dx\) using the midpoint, trapezoidal, and Simpson’s rules.

Solution.

Let \(f(x)=x^x\). The midpoint estimate for the integral is \(2hf(x_{0})\) where \(h=(b-a)/2=1/4\) and \(x_{0}=0.75.\) Then the midpoint estimate, using 6 digits, is \(f(0.75)/2=0.805927/2=0.402964\). The trapezoidal estimate is \(\frac{h}{2}\left[f(0.5)+f(1)\right]\) where \(h=1/2\), which results in \(1.707107/4=0.426777.\) Finally, for Simpson’s rule, \(h=(b-a)/2=1/4,\) and thus the estimate is

\[\begin{equation*} \frac{h}{3}\left[f(0.5)+4f(0.75)+f(1)\right]=\frac{1}{12}\left[0.707107+4(0.805927)+1\right]=0.410901. \end{equation*}\]

Here is a summary of the results:

\[\begin{equation*} \begin{array}{ccc} \text{Midpoint} & \text{Trapezoidal} & \text{Simpson's} \\ \hline 0.402964 & 0.426777 & 0.410901 \\ \end{array} \end{equation*}\]

Example 43

Find the constants \(c_{0},c_{1},x_{1}\) so that the quadrature formula

\[\begin{equation*} \int_{0}^{1}f(x)dx=c_{0}f(0)+c_{1}f(x_{1}) \end{equation*}\]

has the highest possible degree of accuracy.

Solution.

We will find how many of the polynomials \(1,x,x^{2},...\) the rule can integrate exactly. If \(p(x)=1,\) then

\[\begin{equation*} \int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\Rightarrow1=c_{0}+c_{1}. \end{equation*}\]

If \(p(x)=x,\) we get

\[\begin{equation*} \int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\Rightarrow\frac{1}{2}=c_{1}x_{1} \end{equation*}\]

and \(p(x)=x^{2}\) implies

\[\begin{equation*} \int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\Rightarrow\frac{1}{3}=c_{1}x_{1}^{2}. \end{equation*}\]

We have three unknowns and three equations, so we have to stop here. Solving the three equations we get: \(c_{0}=1/4,c_{1}=3/4,x_{1}=2/3.\) So the quadrature rule is of precision two and it is:

\[\begin{equation*} \int_{0}^{1}f(x)dx=\frac{1}{4}f(0)+\frac{3}{4}f(\frac{2}{3}). \end{equation*}\]

Exercise 5.1-2

Find \(c_{0},c_{1},c_{2}\) so that the quadrature rule \(\int_{-1}^{1}f(x)dx=c_{0}f(-1)+c_{1}f(0)+c_{2}f(1)\) has degree of accuracy 2.

Atk89

K.E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons, 2nd edition, 1989.

IK66(1,2)

Isaacson, E. and H.B. Keller. Analysis of Numerical Methods. John Wiley & Sons, 1966.