Mathematica >

NIntegrate Integration Rules

Introduction

An integration rule computes an estimate of an integral over a region using a weighted sum. In the context of NIntegrate usage, an integration rule object provides both an integral estimate and an error estimate as a measure of the integral estimate's accuracy.
An integration rule samples the integrand with a set of points. These points are called sampling points. In the literature these are usually called abscissas. Corresponding to each sampling point xi there is a weight number wi. An integration rule estimates the integral with the weighted sum wif (xi). An integration rule is a functional, that is, it maps functions over the interval [a, b] (or a more general region) into real numbers.
If a rule is applied over the region V this will be denoted as RV (f), where f is the integrand.
The sampling points of the rules considered below are chosen to compute estimates for integrals either over the interval [0, 1], or the unit cube [0, 1]d, or the "centered" unit cube , where d is the dimension of the integral. So if V is one of these regions, R (f) will be used instead of RV (f). When these rules are applied to other regions, their abscissas and estimates need to be scaled accordingly.
The integration rule R is said to be exact for the function f if .
The application of an integration rule R to a function f will be referred as an integration of f, for example, "when f is integrated by R, we get R (f)."
A one-dimensional integration rule is said to be of degree n if it integrates exactly all polynomials of degree n or less, and will fail to do so for at least one polynomial of degree n+1.
A multidimensional integration rule is said to be of degree n if it integrates exactly all monomials of degree n or less, and will fail to do so for at least one monomial of degree n+1, that is, the rule is exact for all monomials of the form , where d is the dimension, i≥ 0, and in.
A null rule of degree m will integrate to zero all monomials of degree m and will fail to do so for at least one monomial of degree m+1. Each null rule may be thought of as the difference between a basic integration rule and an appropriate integration rule of a lower degree.
If the set of sampling points of a rule R1 of degree n, contains the set of sampling points of a rule R2 of a lower degree m, that is, n>m, then R2 is said to be embedded in R1. This will be denoted as R2R1.
An integration rule of degree n that is a member of a family of rules with a common derivation and properties, but different degrees will be denoted as R (f, n), where R might be chosen to identify the family. (For example, trapezoidal rule of degree 4 might be referred to as T (f, 4).)
If each rule in a family is embedded in another rule in the same family, then the rules of that family are called progressive. (For any given m there exists n, n>m, for which R (f, m)R (f, n)).
An integration rule is of open type if the integrand is not evaluated at the end points of the interval. It is of closed type if it uses integrand evaluations at the interval end points.
An NIntegrate integration rule object has one integration rule for the integral estimate and one or several null rules for the error estimate. The sampling points of the integration rule and the null rules coincide. It should be clear from the context whether "integration rule" or "rule" would mean an NIntegrate integration rule object, or an integration rule in the usual mathematical sense.

Integration Rule Specification

All integration rules described below, except "MonteCarloRule", are to be used by the adaptive strategies of NIntegrate. In NIntegrate, all Monte Carlo strategies, crude and adaptive, use "MonteCarloRule". Changing the integration rule component of an integration strategy will make a different integration algorithm.
The way to specify what integration rule the adaptive strategies in NIntegrate (see "Global Adaptive Strategy" and "Local Adaptive Strategy") should use is through a Method suboption.
Here is an example of using an integration rule with a strategy ("GlobalAdaptive").
In[1]:=
Click for copyable input
Out[1]//InputForm=
Here is an example of using the same integration rule as in the example above through a different strategy ("LocalAdaptive").
In[2]:=
Click for copyable input
Out[2]//InputForm=
If NIntegrate is given a method option that has only an integration rule specification other than "MonteCarloRule", then that rule is used with the "GlobalAdaptive" strategy. The two inputs below are equivalent.
For this integration only integration rule is specified.
In[3]:=
Click for copyable input
Out[3]//InputForm=
For this integration an integration strategy and an integration rule are specified.
In[4]:=
Click for copyable input
Out[4]//InputForm=
Similarly for "MonteCarloRule", the adaptive Monte Carlo strategy is going to be used when the following two equivalent commands are executed.
For this Monte Carlo integration only the "MonteCarloRule" is specified.
In[5]:=
Click for copyable input
Out[5]//InputForm=
For this Monte Carlo integration a Monte Carlo integration strategy and "MonteCarloRule" are specified.
In[6]:=
Click for copyable input
Out[6]//InputForm=

"TrapezoidalRule"

The trapezoidal rule for integral estimation is one of the simplest and oldest rules (possibly used by the Babylonians and certainly by the ancient Greek mathematicians):
The compounded trapezoidal rule is a Riemann sum of the form
where .
If the Method option is given the value "TrapezoidalRule", the compounded trapezoidal rule is used to estimate each subinterval formed by the integration strategy.
A "TrapezoidalRule" integration:
In[7]:=
Click for copyable input
Out[7]=
option name
default value
"Points"5number of coarse trapezoidal points
"RombergQuadrature"Trueshould Romberg quadrature be used or not
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"TrapezoidalRule" options.

The trapezoidal rule and its compounded (multipanel) extension are not very accurate. (The compounded trapezoidal rule is exact for linear functions and converges at least as fast as n-2, if the integrand has continuous second derivative [DavRab84].) The accuracy of the multipanel trapezoidal rule can be increased using the "Romberg quadrature".
Since the abscissas of T (f, n) are a subset of T (f, 2n-1), the difference |T (f, 2 n-1)-T (f, n)|, can be taken to be an error estimate of the integral estimate T (f, 2n-1), and can be computed without extra integrand evaluations.
The option "Points"->k can be used to specify how many coarse points are used. The total number of points used by "TrapezoidalRule" is 2k-1.
This verifies that the sampling points are as in (1).
In[8]:=
Click for copyable input
Out[9]=
"TrapezoidalRule" can be used for multidimensional integrals too.
Here is a multidimensional integration with "TrapezoidalRule". The exact result is .
In[10]:=
Click for copyable input
Out[10]=
Remark: NIntegrate has both a trapezoidal rule and a trapezoidal strategy; see "Trapezoidal" Strategy in the tutorial Integration Strategies. All internally implemented integration rules of NIntegrate have the suffix -Rule. So "TrapezoidalRule" is used to specify the trapezoidal integration rule, and "Trapezoidal" is used to specify the trapezoidal strategy.

Romberg Quadrature

The idea of the Romberg quadrature is to use a linear combination of T (f, n) and T (f, 2 n-1) that eliminates the same order terms of truncation approximation errors of T (f, n) and T (f, 2 n-1).
From the Euler-Maclaurin formula [DavRab84] we have
where
Hence we can write
The h2 terms of the equations above can be eliminated if the first equation is subtracted from the second equation four times. The result is
This example shows that a trapezoidal rule using the Romberg quadrature gives better performance than the standard trapezoidal rule. Also, the result of the former is closer to the exact result, .
In[11]:=
Click for copyable input
Out[11]=
Here is an integration with a trapezoidal rule that does not use Romberg quadrature.
In[10]:=
Click for copyable input
Out[10]=

"TrapezoidalRule" Sampling Points and Weights

The following calculates the trapezoidal sampling points, weights, and error weights for a given precision.
In[3]:=
Click for copyable input
Out[4]=
Here is how the Romberg quadrature weights and error weights can be derived.
In[5]:=
Click for copyable input
Out[9]=

"NewtonCotesRule"

Newton-Cotes integration formulas are formulas of interpolatory type with sampling points that are equally spaced.
The Newton-Cotes quadrature for NIntegrate can be specified with the Method option value "NewtonCotesRule".
In[20]:=
Click for copyable input
Out[20]=
option name
default value
"Points"3number of coarse Newton-Cotes points
"Type"Closedtype of the Newton-Cotes rule
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"NewtonCotesRule" options.

Let the interval of integration, [a, b], be divided into n-1 subintervals of equal length by the points
Then the integration formula of interpolatory type is given by
where
with
When n is large, the Newton-Cotes n-point coefficients are large and are of mixed sign.
In[21]:=
Click for copyable input
Out[21]=
Since this may lead to large losses of significance by cancellation, a high-order Newton-Cotes rule must be used with caution.

"NewtonCotesRule" Sampling Points and Weights

The following calculates the Newton-Cotes sampling points, weights, and error weights for a given precision.
In[22]:=
Click for copyable input
Out[23]=

"GaussBerntsenEspelidRule"

Gaussian quadrature uses optimal sampling points (through polynomial interpolation) to form a weighted sum of the integrand values over these points. On a subset of these sampling points a lower order quadrature rule can be made. The difference between the two rules can be used to estimate the error. Berntsen and Espelid derived error estimation rules by removing the central point of Gaussian rules with odd number of sampling points.
The Gaussian quadrature for NIntegrate can be specified with the Method option value "GaussBerntsenEspelidRule".
In[24]:=
Click for copyable input
Out[24]=
option name
default value
"Points"Automaticnumber of Gauss points
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"GaussBerntsenEspelidRule" options.

A Gaussian rule G (f, n) of n points for integrand f is exact for polynomials of degree 2n-1 (i.e., if f (x) is a polynomial of degree ≤2n-1).
Gaussian rules are of open type since the integrand is not evaluated at the end points of the interval. (Lobatto rules, Clenshaw-Curtis rules, and the trapezoidal rule are of closed type since they use integrand evaluations at the interval end points.)
This defines the divided differences functional [Ehrich2000]
For the Gaussian rule G (f, 2n+1), with sampling points x1, x2, ..., x2n+1, Berntsen and Espelid have derived the following error estimate functional (see [Ehrich2000])
(The original formula in [Ehrich2000] is for sampling points in [-1, 1]. The formula above is for sampling points in [0, 1].)
This example shows the number of sampling points used by NIntegrate with various values of "GaussBerntsenEspelidRule"'s option "Points".
In[25]:=
Click for copyable input
Out[25]=

"GaussBerntsenEspelidRule" Sampling Points and Weights

The following calculates the Gaussian abscissas, weights, and Bernsen-Espelid error weights for a given number of coarse points and precision.
In[26]:=
Click for copyable input
Out[27]=
The Berntsen-Espelid error weights are implemented below.
This implements the divided differences.
In[28]:=
Click for copyable input
This computes the abscissas and the weights of G (f, 2n+1).
In[30]:=
Click for copyable input
This computes the Berntsen-Espelid error weights.
In[31]:=
Click for copyable input
Out[31]=

"GaussKronrodRule"

Gaussian quadrature uses optimal sampling points (through polynomial interpolation) to form a weighted sum of the integrand values over these points. The Kronrod extension of a Gaussian rule adds new sampling points in between the Gaussian points and forms a higher-order rule that reuses the Gaussian rule integrand evaluations.
The Gauss-Kronrod quadrature for NIntegrate can be specified with the Method option value "GaussKronrodRule".
In[32]:=
Click for copyable input
Out[32]=
option name
default value
"Points"Automaticnumber of Gauss points that will be extended with Kronrod points
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

"GaussKronrodRule" options.

A Gaussian rule G (f, n) of n points for integrand f is exact for polynomials of degree 2n-1, that is, if f (x) is a polynomial of degree ≤2n-1.
Gauss-Kronrod rules are of open type since the integrand is not evaluated at the end points of the interval.
The Kronrod extension GK (f, n) of a Gaussian rule with n points G (f, n) adds n+1 points to G (f, n) and the extended rule is exact for polynomials of degree 3n+1 if n is even, or 3n+2 if n is odd. The weights associated with a Gaussian rule change in its Kronrod extension.
Since the abscissas of G (f, n) are a subset of GK (f, n), the difference GK (f, n)-G (f, n) can be taken to be an error estimate of the integral estimate GK (f, n), and can be computed without extra integrand evaluations.
This example shows the number of sampling points used by NIntegrate with various values of "GaussKronrodRule"'s option "Points".
In[33]:=
Click for copyable input
Out[33]=
For an implementation description of Kronrod extensions of Gaussian rules, see [PiesBrand74].

"GaussKronrodRule" Sampling Points and Weights

The following calculates the Gauss-Kronrod abscissas, weights, and error weights for a given number of coarse points and precision.
In[34]:=
Click for copyable input
Out[35]=
The calculations below demonstrate the degree of the Gauss-Kronrod integration rule (see above).
This computes the degree of the Gauss-Kronrod integration rule.
In[36]:=
Click for copyable input
Out[36]=
This defines a function.
In[37]:=
Click for copyable input
The command below implements the integration rule weighted sums for the integral estimate, wi f (xi), and the error estimate, ei f (xi), where are the abscissas, are the weights, and are the error weights.
These are the integral and error estimates for computed with the rule.
In[38]:=
Click for copyable input
Out[38]=
The integral estimate coincides with the exact result.
In[39]:=
Click for copyable input
Out[39]=
The error estimate is not zero since the embedded Gauss rule is exact for polynomials of degree ≤2 n-1. If we integrate a polynomial of that degree, the error estimate becomes zero.
This defines a function.
In[40]:=
Click for copyable input
These are the integral and error estimates for computed with the rule.
In[41]:=
Click for copyable input
Out[41]=
Here is the exact result using Integrate.
In[42]:=
Click for copyable input
Out[42]=

"LobattoKronrodRule"

The Lobatto integration rule is a Gauss type rule with preassigned abscissas. It uses the end points of the integration interval and optimal sampling points inside the interval to form a weighted sum of the integrand values over these points. The Kronrod extension of a Lobatto rule adds new sampling points in between the Lobatto rule points and forms a higher-order rule that reuses the Lobatto rule integrand evaluations.
NIntegrate uses the Kronrod extension of the Lobatto rule if the Method option is given the value "LobattoKronrodRule".
In[43]:=
Click for copyable input
Out[43]=
option name
default value
"Points"5number of Gauss-Lobatto points that will be extended with Kronrod points
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"LobattoKronrodRule" options.

A Lobatto rule L (f, n) of n points for integrand f is exact for polynomials of degree 2n-3, (i.e., if f (x) is a polynomial of degree ≤2n-3).
The Kronrod extension LK (f, n) of a Lobatto rule with n points L (f, n) adds n-1 points to L (f, n) and the extended rule is exact for polynomials of degree 3n-2 if n is even, or 3n-3 if n is odd. The weights associated with a Lobatto rule change in its Kronrod extension.
As with "GaussKronrodRule", the number of Gauss points is specified with the option "GaussPoints". If "LobattoKronrodRule" is invoked with "Points"->n, the total number of rule points will be 2 n -1.
This example shows the number of sampling points used by NIntegrate with various values of "LobattoKronrodRule"'s option "Points".
In[44]:=
Click for copyable input
Out[44]=
Since the Lobatto rule is a closed rule, the integrand needs to be evaluated at the end points of the interval. If there is a singularity at these end points, NIntegrate will ignore it.
For an implementation description of Kronrod extensions of Lobatto rules, see [PiesBrand74].

"LobattoKronrodRule" Sampling Points and Weights

The following calculates the Lobatto-Kronrod abscissas, weights, and error weights for a given number of coarse points and precision.
In[45]:=
Click for copyable input
Out[46]=
The calculations below demonstrate the degree of the Lobatto-Kronrod integration rule (see above).
This computes the degree of the Lobatto-Kronrod integration rule.
In[47]:=
Click for copyable input
Out[47]=
This defines a function.
In[48]:=
Click for copyable input
The command below implements the integration rule weighted sums for the integral estimate, wi f (xi), and the error estimate, ei f (xi), where are the abscissas, are the weights, and are the error weights.
These are the integral and error estimates for computed with the rule.
In[49]:=
Click for copyable input
Out[49]=
The preceding integral estimate coincides with the exact result.
In[50]:=
Click for copyable input
Out[50]=
The preceding error estimate is not zero since the embedded Lobatto rule is exact for polynomials of degree ≤2 n-3. If we integrate a polynomial of that degree, the error estimate becomes zero.
This defines a function.
In[51]:=
Click for copyable input
These are the integral and error estimates for computed with the rule.
In[52]:=
Click for copyable input
Out[52]=
The exact result using Integrate.
In[53]:=
Click for copyable input
Out[53]=

"ClenshawCurtisRule"

A Clenshaw-Curtis rule uses sampling points derived from the Chebyshev polynomial approximation of the integrand.
The Clenshaw-Curtis quadrature for NIntegrate can specified with the Method option value "ClenshawCurtisRule".
In[54]:=
Click for copyable input
Out[54]=
option name
default value
"Points"5number of coarse Clenshaw-Curtis points
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

"ClenshawCurtisRule" options.

Theoretically a Clenshaw-Curtis rule with n sampling points is exact for polynomials of degree n or less. In practice, though, Clenshaw-Curtis rules achieve the accuracy of the Gaussian rules [Evans93][OHaraSmith68]. The error of the Clenshaw-Curtis formula is analyzed in [OHaraSmith68].
The sampling points of the classical Clenshaw-Curtis rule are zeros of Chebyshev polynomials. The sampling points of a practical Clenshaw-Curtis rule are chosen to be Chebyshev polynomial extremum points. The classical Clenshaw-Curtis rules are not progressive, but the practical Clenshaw-Curtis rules are [DavRab84][KrUeb98].
Let PCC (f, n) denote a practical Clenshaw-Curtis rule of n sampling points for the function f.
The progressive property means that the sampling points of PCC (f, n) are a subset of the sampling points of PCC (f, 2 n-1). Hence the difference |PCC (f, 2 n-1)-PCC (f, n) can be taken to be an error estimate of the integral estimate PCC (f, 2 n-1), and can be computed without extra integrand evaluations.
The NIntegrate option Method->{"ClenshawCurtisRule", "Points"->k} uses a practical Clenshaw-Curtis rule with 2n-1 points PCC (f, 2 n-1).
In[55]:=
Click for copyable input
Out[55]=
This example shows the number of sampling points used by NIntegrate with various values of "ClenshawCurtisRule"'s option "Points".
In[56]:=
Click for copyable input
Out[56]=

"ClenshawCurtisRule" Sampling Points and Weights

Here are the sampling points and the weights of the Clenshaw-Curtis rule for a given coarse number of points and precision.
In[57]:=
Click for copyable input
Out[58]=
Here is another way to compute the sampling points of PCC (f, 2 n-1).
In[59]:=
Click for copyable input
Out[60]=
This defines a function.
In[61]:=
Click for copyable input
These are the integral and error estimates for computed with the rule.
In[62]:=
Click for copyable input
Out[62]=
The exact value by Integrate.
In[63]:=
Click for copyable input
Out[63]=

"MultiPanelRule"

"MultiPanelRule" combines into one rule the applications of a one-dimensional integration rule over two or more adjacent intervals. An application of the original rule to any of the adjacent intervals is called a panel.
Here is an example of an integration with "MultiPanelRule".
In[64]:=
Click for copyable input
Out[64]=
option name
default value
Method"NewtonCotesRule"integration rule specification that provides the abscissas, weights, and error weights for a single panel
"Panels"5number of panels
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

"MultiPanelRule" options.

Let the unit interval [0, 1] be partitioned into k sub-intervals with the points 0=y0<y1<...<yk=1.
If we have the rule
it can be transformed into a rule for the interval [yj-1, yj],
Let xi j=xi (yj-yj-1)+yj-1, and yj-yj-1=1/k, j =1, ..., k. Then the k-panel integration rule based on R (f) can be written explicitly as
If R (f) is closed, that is, R (f) has 0 and 1 as sampling points, then xn j-1=x1 j, and the number of sampling points of k×R (f) can be reduced to k (n-1)+1. (This is done in the implementation of "MultiPanelRule".)
More about the theory of multi-panel rules, also referred to as compounded or composite rules, can be found in [KrUeb98] and [DavRab84].

"MultiPanelRule" Sampling Points and Weights

The sampling points and the weights of the "MultiPanelRule" can be obtained with this command.
In[65]:=
Click for copyable input
Out[66]=
Here are the abscissas and weights of a Gauss-Kronrod rule.
In[67]:=
Click for copyable input
Out[67]=
The multi-panel rule abscissas can be obtained using Rescale.
In[68]:=
Click for copyable input
Out[68]=
This shows how to derive the multi-panel rule weights from the original weights.
In[69]:=
Click for copyable input
Out[69]=

"CartesianRule"

A d-dimensional Cartesian rule has sampling points that are a Cartesian product of the sampling points of d one-dimensional rules. The weight associated with a Cartesian rule sampling point is the product of the one-dimensional rule weights that correspond to its coordinates.
The Cartesian product integration for NIntegrate can be specified with the Method option value "CartesianRule".
In[70]:=
Click for copyable input
Out[70]=
option name
default value
Method"GaussKronrodRule"a rule or a list of rules with which the Cartesian product rule will be formed
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"CartesianRule" options.

For example, suppose we have the formulas:
that are exact for polynomials of degree d1, d2, and d3, respectively. Then it is not difficult to see that the formula with n1×n2×n3 points,
is exact for polynomials in x1, x2, x3 of degree min (d1, d2, d3). Note that the weight associated with the abscissa is .
The general Cartesian product formula for D one-dimensional rules the i of which has ni sampling points and weights is
Clearly (2) can be written as
where n=nk , and for each integer k[1, n], and .
Here is a visualization of a Cartesian product rule integration. Along the x axis "TrapezoidalRule" is used; along the y axis "GaussKronrodRule" is used.
In[71]:=
Click for copyable input
Out[72]=
Cartesian rules are applicable for relatively low dimensions (≤ 4), since for higher dimensions they are subject to "combinatorial explosion." For example, a 5-dimensional Cartesian product of 5 identical one-dimensional rules each having 10 sampling points would have 10^5 sampling points.
NIntegrate uses Cartesian product rule if the integral is multidimensional and the Method option is given a one-dimensional rule or a list of one-dimensional rules.
Here is an example specifying Cartesian product rule integration with GaussKronrodRule.
In[73]:=
Click for copyable input
Out[73]=
Here is an example specifying Cartesian product rule integration with a list of one-dimensional integration rules.
In[74]:=
Click for copyable input
Out[74]=
Another example specifying Cartesian product rule integration with a list of one-dimensional integration rules.
In[75]:=
Click for copyable input
Out[75]=
More about Cartesian rules can be found in [Stroud71].

"CartesianRule" Sampling Points and Weights

The sampling points and the weights of the "CartesianRule" rule can be obtained with the command NIntegrate`CartesianRuleData.
In[76]:=
Click for copyable input
Out[76]=
NIntegrate`CartesianRuleData keeps the abscissas and the weights of each rule separated. Otherwise, as it can be seen from (3) the result might be too big for higher dimensions.
The results of NIntegrate`CartesianRuleData can be put into the form of (4) with this function.
In[77]:=
Click for copyable input
In[78]:=
Click for copyable input
Out[78]=

"MultiDimensionalRule"

A fully symmetric integration rule for the cube , d, d>1 consists of sets of points with the properties: (i) all points in a set can be generated by permutations and/or sign changes of the coordinates of any fixed point from that set, (ii) all points in a set have the same weight associated with them.
The fully symmetric multidimensional integration (fully symmetric cubature) for NIntegrate can be specified with the Method option value "MultiDimensionalRule".
In[79]:=
Click for copyable input
Out[79]=
A set of points of a fully symmetric integration rule that satisfies the preceding properties is called an orbit. A point of an orbit, {x1, x2, ..., xd}, for the coordinates of which the inequality x1x2≥...≥ xd holds, is called a generator. (See [KrUeb98][GenzMalik83].)
option name
default value
"Generators"5number of generators of the fully symmetric rule
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

"MultiDimensionalRule" options.

If an integration rule has K orbits denoted 1, 2, ..., K, and the ith of them, i, has a weight wi associated with it, then the integral estimate is calculated with the formula
A null rule of degree m will integrate to zero all monomials of degree m and will fail to do so for at least one monomial of degree m+1. Each null rule may be thought of as the difference between a basic integration rule and an appropriate integration of lower degree.