The Newton-Cotes formulas, also called as the Newton-Cotes quadrature rules or simply the Newton-Cotes rules, are a group of formulas for numerical integration which are based on evaluating the integrand at equally spaced points. They are named after Isaac Newton and Roger Cotes respectively. Newton-Cotes formulas can be very useful if the values of the integrand at equally spaced points are available. If the set of points are not equally spaced, then other methods such as Gaussian quadrature are probably more accurate. There are many derivations of the Newton-Cotes formulas. The two most famous ones are the trapezoidal quadrature and the Simpson quadrature formulas.
The Trapezoidal Quadrature Formula
The trapezoidal quadrature formula works by replacing the function using a straight line, basically approximating the region under the function as a trapezoid and calculating its area.
The line equation passing through the two points \( a \) and \( b \) is
\( \displaystyle L_1(x)=\frac{f(b)-f(a)}{b-a}(x-a)+f(a) \)
The area under the line is the area of the trapezoid, which is simply multiplying the sum of the bases which are the parallel lines by the height of the trapezoid, which is the perpendicular distance between the bases, and then divide by 2. We can also do a proper integration to get the same result. Now replacing \( f(x) \) with \( L_1(x) \), we have
\( \int_{a}^{b}f(x)dx\approx\int^{b}_{a}L_1(x)dx \)
\( \displaystyle\int_{a}^{b}f(x)dx=\int_{a}^{b}\Big[\frac{f(b)-f(a)}{b-a}(x-a)+f(a)\Big] \)
\( \displaystyle\int_{a}^{b}f(x)dx=\frac{f(a)+f(b)}{2}(b-a) \)
For a function \( f(x) \) and an interval \( [a,b] \), we can divide the interval into \( N \) sub-intervals. Then we apply the trapezoidal rule to each sub-interval and then sum all of then up.
Let
\(\displaystyle \Delta x=\frac{b-a}{N} \)
\( \displaystyle\int_{a}^{b}f(x)dx\approx\frac{\Delta x}{2}\sum_{k=1}^{N}(f(x_{k-1})+f(x_k)) \)
The error form of trapezoidal rule is given by
\( \displaystyle \varepsilon=-\frac{(b-a)^3}{12N}f”(\xi) \)
for some \( \xi\epsilon[a,b] \).
The following example solves the problem in the above figure of the trapezoidal formula.
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
return -(x * x- 4 * x + 6); // -(x^2 -4x + 6)
}
};
//the limit
double a = 0, b = 4;
// an integrator using the trapezoidal rule
Integrator integrator = new Trapezoidal(1e-8, 20); // precision, max number of iterations
// the integration
double I = integrator.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %f", a, b, I));
The output is:
S_[0,4] f(x) dx = -13.333333
The Simpson Quadrature Formula
The Simpson quadrature formula, instead of replacing the function by a line, replaces \( f(x) \) by a quadratic function or a parabola. teh interval is divided into two prtitions giving three points: \( a,\,b \) and a mid-point \( m=\frac{a+b}{2} \). The quadratic function is given by
\( \displaystyle L_2(x)=\frac{(x-m)(x-b)}{(a-m)(a-b)}f(a)+\frac{(x-a)(x-b)}{(m-a)(m-b)}f(m)+\frac{(x-m)(x-a)}{(b-m)(b-a)}f(b) \)
From a geometrical point of view, the Simpson formula calculates the area of the curved trapezoid surrounded by a parabola to approximate the area enclosed by \( f(x) \).
Replacing \( f(x) \) with \( L_2(x) \), we have
\( \int_{a}^{b}f(x)dx\approx\int_{a}^{b}L_2(x)dx \)
\( \displaystyle \int_{a}^{b}f(x)dx=\frac{b-a}{6}\Big[f(a)+4f\Big(\frac{a+b}{2}\Big)+f(b)\Big] \)
\( \displaystyle \int_{a}^{b}f(x)dx=\frac{h}{3}\Big[f(a)+4f\Big(\frac{a+b}{2}\Big)+f(b)\Big] \)
\( h=\frac{b-a}{2} \) is the step size in each partition. Because of the \( \frac{1}{3} \) in the formula, Simpson’s quadrature formula is also called as Simpson’s 1/3 rule.
For a function \( f(x) \) and an interval \( [a,b] \), we can divide the interval into \( N \) sub-intervals. Then we can apply the Simpson rule to each of the subinterval and then sum them all up. We have
\( \displaystyle \int_{a}^{b}f(x)dx\approx\frac{h}{3}\sum^{N/2}_{j=1}[f(x_{2j-2})+4f(x_{2j-1})+f(x_{2j})] \)
Simson’s rule is expected to improve on the trapezoidal rule for the functions which are twice continuously differentiable. However, for the rougher functions the trapezoidal rule is likely to be more preferable than the Simpson’s rule. The error term for the Simpson’s rule is given by
\( \displaystyle \varepsilon=-\frac{h^4}{180}(b-a)f^{(4)}(\xi) \)
for some \( \xi\epsilon[a,b] \).
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
return -(x * x- 4 * x + 6); // -(x^2 -4x + 6)
}
};
//the limit
double a = 0, b = 4;
// an integrator using the Simpson rule
Integrator integrator = new Simpson(1e-8, 20); // precision, max number of iterations
// the integration
double I = integrator.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %f", a, b, I));
The output is:
S_[0,4] f(x) dx = -13.333333
The Newton-Cotes Quadrature Formula
It is possible that we can divide the interval \( [a,b] \) into more than 1 partitions, as in the case of trapezoidal rule or more than 2 partitions, as in the case of the Simpson’s rule. One can divide it into \( n \) intervals. One can also use a closed formula where the endpoints \( a \) an \( b \) are used, or an open formula where the endpoints are not used because \( f(a) \) and \( f(b) \) may not exist. The Newton-Cotes formula is such a generalization. It divides the interval \( [a,b] \) into \( n \) equal partitions. For a closed formula, the equally spaced points are given by
\( x_i-a+ih,\,h=\frac{b-a}{n+2} \)
For an open formula, the equally spaced points are
\( x_i=a+(i+1),\,h=\frac{b-a}{n+2} \)
For the closed formula, an n-degree difference polynomial can be constructed,
\( \displaystyle L_n(x)=\sum_{k=0}^{n}\frac{w(x)}{(x-x_k)w'(x_k)}f(x_k) \)
where \( w(x)=(x-x_0)(x-x_1)\hdots(x-x_n) \)
Now replacing \( f(x) \) with \( L_n(x) \), we have
\( \displaystyle\int_{a}^{b}f(x)dx\approx\int_{a}^{b}L_n(x)dx=\int_{a}^{b}\Big(\sum_{k=0}^{n}\frac{w(x)}{(x-x_k)w'(x_k)}f(x_k)\Big)dx \)
\( \displaystyle \int_{a}^{b}L_n(x)dx=\sum_{i=0}^{n}\Big(\int_{a}^{b}\frac{w(x)}{(x-x_k)w'(x_k)}dx\Big)f(x_k) \)
\( \displaystyle\int_{a}^{b}L_n(x)dx=\sum_{i=0}^{n}A_kf(x_k) \)
where the coefficients \( \displaystyle A_k=\int_{a}^{b}\frac{w(x)}{(x-x_k)w'(x_k)}dx \)
The above formula is called the Newton-Cotes formula. The coeffients are called the Newton-Cotes coefficients. To use the formuls, we need to calculate the coefficients \( A_k \). Substituting \( x=a+th \), we have
\( w(x)=w(a+th)=h^{n+1}t(t-1)\hdots(t-n) \)
\( w'(x)=h^n(-1)^{n-1}(n!)(n-k)! \)
\( \displaystyle A_k=\int_{a}^{b}\frac{w(x)}{(x-x_k)w'(x_k)}dx \)
\( \displaystyle A_k=\int_{0}^{n}\frac{h^{n+1}t(t-1)\hdots(t-n)}{h^n(-1)^{n-1}(k!)(n-k)!h(t-k)}hdt \)
Let,
\( \displaystyle C_k^{(n)}=\frac{(-1)^{n-k}}{n(k!)(n-k)!}\int_{0}^{n}\frac{t(t-1)\hdots(t-n)}{(t-k)}dt \)
So, we have
\( A_k=(b-a)C_k^{(n)} \)
From the above equation we can say that, \( A_k \) is a constant that depends only on the interval \( [a,b] \) and the number of positions \( n \), but not on the function. Hence, the Newton-Cotes coefficients can eb computed in advanced because of this particular reason.
The Trapezoidal and Simpson’s formula are special cases of the Newton-Cotes formula. When \( n=1 \), The Newton-Cotes formula is given as
\( \displaystyle L_1(x)=\frac{b-a}{2}[f(a)+f(b)] \)
The above equation is the Trapezoidal Formula.
When \( n=2 \), the Newton-Cotes formula becomes
\( \displaystyle L_2(x)=\frac{b-a}{6}[f(a)+4f\Big(\frac{a+b}{2}\Big)+f(b)] \)
This is the Simpson’s formula.
When \( n=4 \), the Newrton-Cotes formula is
\( \displaystyle C=\frac{b-a}{90}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)] \)
and \( x_i=a+kh\,(k=0,1,2,3,4),\,h=\frac{b-a}{4} \)
Note that the Newton-Cotes coefficienst \( A_k \) depends on \( x \) but not the integrand function iteself. For \( n=1 \), the trapeoidal formulas is polynomial of degree 1. If \( f(x) \) in the Newton-Cotes formula is a polynomial with degree nop higher than \( n \), then the equal sign holds.
\( \int_{a}^{b}f(x)dx\approx\sum_{k=0}^{n}A_kf(x_k) \)
If \( f(x) \) is a polynomial of degree \( n+1 \), the right hand side is only an approximation. the quadrature formula is said to have an \( n \) degree algebraic accuracy.
For \( n=1 \), the trapezoidal quadrature rule has thsi error and is said to have the primary algebraic precision.
\( \displaystyle \varepsilon=-\frac{(b-a)^3}{12N}f”(\xi),\,a\leq\xi\leq b \)
For \( n=2 \), the Simpson’s quadrature rule has this error and si said to have the cubic algebraic precisio.
\( \displaystyle \varepsilon=-\frac{h^4}{180}(b-a)f^{(4)}(\xi),\,a\leq\xi\leq b \)
For a general \( n \), the error term is
\( \displaystyle \varepsilon=\frac{f^{(n+1)}(\xi)}{(n+1)!}w(x),\,a\leq\xi\leq b \)
Assuming that the rounding error of \( f(x_i) \) where \( x_i \) are the grid points is \( \varepsilon_i \), \( \varepsilon_{max}=\underset{o\leq i\leq n}{\mathrm{max}}|\varepsilon_i| \), the error of the Newton-Cotes formula is given as,
\( \displaystyle\varepsilon=|b-a|.\Big|\sum_{i=0}^{n}C_i^{(n)}\varepsilon_i\Big|\leq|(b-a)|.\varepsilon_{max}.\Big|\sum_{i=0}^{n}C_i^{(n)}\Big| \)
The last term \( \sum_{i=0}^{n}|C_i^{(n)}| \) increases with \( n \). This will cause the eror of the approximation to increase, so the Newton-Cotes formula is nstable for big \( n \). therefore, in practice the Newton-Cotes formula is rarely used for \( n\geq 8 \).
Consider the following example of integrating the following function,
\( \displaystyle I+\int_{0}^{1}ydx=\int_{0}^{1}\frac{4}{1+x^2}dx=\pi \)
Using the trapezoid formula, we have
\( \displaystyle \int_{0}^{1}ydx\approx\frac{1-0}{2}\Big[\frac{4}{1+1}+\frac{4}{1+0}\Big]=3 \)
using the Simpson’s formula, we have
\( \displaystyle \int_{0}^{1}ydx\approx\frac{1-0}{6}\Big[\frac{4}{1+1}+4\frac{4}{1+\frac{1}{4}}+\frac{4}{1+0}\Big]\approx 3.1333 \)
using the Newton-Cotes formula with \( n=4 \), we have
\( x_0=0,\quad x_1=0.2500,\quad x_2=0.5000\quad x_3=0.7500, \quad x_4=1.0000 \)
\( \displaystyle \int_{0}^{1}ydx\approx\frac{b-a}{90}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)]\approx 3.1421 \)
System.out.println("integrate using the Newton-Cotes formulas");
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
return 4 / (1 + x * x); // 4/(1+x^2)
}
};
//the limit
double a = 0, b = 1;
Integrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator2 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new NewtonCotes(3, Type.CLOSED, 1e-8, 20); // using the Newton-cotes rule
Integrator integrator4 = new NewtonCotes(3, Type.OPEN, 1e-8, 20); // using the Newton-cotes rule
// the integrations
double I1 = integrator1.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f", a, b, I3));
double I4 = integrator4.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f", a, b, I4));
The output is:
S_[0,1] f(x) dx = 3.1415926436556850, using the trapezoidal rule
S_[0,1] f(x) dx = 3.1415926535528365, using the Simpson rule
S_[0,1] f(x) dx = 3.1415926497180000, using the Newton_Cotes closed rule
S_[0,1] f(x) dx = 3.1415926674370604, using the Newton_Cotes open rule
The accuracy and performance of an integration depends not only on which quadrature formula to use, but also on the integrand function itself. In general, the trapezoidal formula is not as accurate as the Simpson’s formula. The low-order Newton-Cotes formula are simple to compute, convenient to use and have high degree of accuracy. On the other hand, the higher-order Newton-Cotes formulas are not only complicated to calculate, but also have poor stability. So, they are rarely used. In practice, we in general prefer low-order formulas like the Simpson’s rule.