Skip to main content

Section 14.3 Newton-Cotes rules

One can create rules of any given order of accuracy, using enough sample points. Usually such rules are derived for the convenient interval \([-1,1]\text{,}\) and then applied to other intervals by change of variable.

Let \((x_1,\dots, x_n)\) be some points in the interval \([-1,1]\text{.}\) We want to come up with an integration rule

\begin{equation} \int_{-1}^1 f(x)\,dx \approx w_1 f(x_1) + \dots + w_n f(x_n)\label{eq-cotes1}\tag{14.3.1} \end{equation}

where \(w_k\) are coefficients (weights). For example, Simpson's rule (in its simple form) has evaluation points \(-1,0,1\) and weights \(1/3, 4/3, 1/3\text{.}\)

How to find the weights that make the rule (14.3.1) as accurate as possible? We can require both sides to be equal when \(f\) is each of the following: \(x^0, x^1, \dots, x^{n-1}\text{.}\) This gives \(n\) linear equations with \(n\) unknowns \(w_j\text{.}\) Then solve the system, perhaps using Matlab. The resulting rules are called Newton-Cotes rules when the points are distributed at equal intervals. They include trapezoidal rule (2 points) and Simpson's rule (3 points) as special cases.

Let evaluation points on \([-1, 1]\) be \(-1, -1/2, 0, 1/2, 1\text{.}\) What will their weights be?

Solution

We need 5 linear equations involving integrals of \(x^0, x^1, \dots, x^4\text{.}\) If we number them by index \(i\) running from 1 to 5, then the equation has \(\int_{-1}^1 x^{i-1}\,dx\) on the right hand side, which is \((1-(-1)^i)/i\text{.}\) On the left we have the sum \(\sum_j w_j x_j^{i-1}\) which means the coefficients of the unknown \(w_j\) is \(x_j^{i-1}\text{.}\) We can create a matrix using outer exponentiation of a row vector by a column vector.

n = 5;
x = linspace(-1, 1, n);
i = (1:n)'; 
A = x.^(i-1);
b = (1-(-1).^i)./i;
w = A\b;
disp(rats(w'))

Here rats is used to display the solution as rational numbers (which they are) rather than the usual decimal approximation.

Having computed the weights as above, we can integrate any function f on [-1, 1] simply by executing f(x)*w, which is the dot product of a vector with function values with the vector of weights. Try this for the exponential function.

Although one can get rules of arbitrarily high order in this way, the gain in practice does not justify the increasing complexity. (Try computing the 9-point rule, for example). It turns out that using many equally spaced points on an interval is generally a bad idea. We need a smarter way of choosing the points, which we will start working on next time.