The Mathematica function NIntegrate uses as one of its methods a fairly sophisticated Gauss-Kronrod based algorithm. Other types of quadrature formulas exist, each with their own advantages. For example, Gaussian quadrature uses values of the integrand at oddly spaced abscissas. If you want to integrate a function presented in tabular form at equally spaced abscissas, it won't work very well. An alternative is to use Newton-Cotes quadrature.
The basic idea behind Newton-Cotes quadrature is to approximate the value of an integral as a linear combination of values of the integrand evaluated at equally spaced points:
In addition, there is the question of whether or not to include the end points in the sum. If they are included, the quadrature formula is referred to as a "closed" formula. If not it is an "open" formula. If the formula is open there is some ambiguity as to where the first abscissa is to be placed. The open formulas given in this package have the first abscissa one half step from the lower end point.
Since there are free parameters to be chosen (the weights) and since both integration and the sum are linear operations, we can expect to be able to make the formula correct for all polynomials of degree less than about
. In addition to knowing what the weights are, it is often desirable to know how large the error in the approximation will be. This package allows you to answer both of these questions.
Finding formulas for Newton-Cotes quadrature.
This loads the package.
Option for NewtonCotesWeights and NewtonCotesError.
Here are the weights and abscissas for the five-point closed Newton-Cotes quadrature formula on the interval
In:= NewtonCotesWeights[5, -3, 7]
Here is the error in that formula. Unfortunately it involves the sixth derivative of
at an unknown point so we don't really know what the error itself is.
In:= NewtonCotesError[5, f, -3, 7]
We can see that the error decreases rapidly with the length of the interval.
In:= NewtonCotesError[5, f, a, a+h]
This gives the weights and abscissas for the five-point open Newton-Cotes quadrature formula on the interval (
In:= NewtonCotesWeights[5, -3, 7, QuadratureType -> Open]
Here is the error in that formula.
In:= NewtonCotesError[5, f, -3, 7, QuadratureType -> Open]