Newton-Cotes rules

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)\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.

🔗

Example 14.3.1. Find the weights for a 5-point Newton-Cotes formula.

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.

🔗

Prev Top Next