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
x_{i} there is a weight number
w_{i}. An integration rule estimates the integral
with the weighted sum
w_{i}f (x_{i}). 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
R^{V} (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
R^{V} (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 onedimensional 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
_{i}≤n.
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
R_{1} of degree
n, contains the set of sampling points of a rule
R_{2} of a lower degree
m, that is,
n>m, then
R_{2} is said to be embedded in
R_{1}. This will be denoted as
R_{2}R_{1}.
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").
Out[1]//InputForm= 
 

Here is an example of using the same integration rule as in the example above through a different strategy ( "LocalAdaptive").
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 an integration strategy and an integration rule are specified.
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.
Out[5]//InputForm= 
 

For this Monte Carlo integration a Monte Carlo integration strategy and "MonteCarloRule" are specified.
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:
Out[7]=  

  
"Points"  5  number of coarse trapezoidal points 
"RombergQuadrature"  True  should Romberg quadrature be used or not 
"SymbolicProcessing"  Automatic  number 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, 2n1), the difference
T (f, 2 n1)T (f, n), can be taken to be an error estimate of the integral estimate
T (f, 2n1), 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
2k1.
This verifies that the sampling points are as in ( 1).
Out[9]=  

"TrapezoidalRule" can be used for multidimensional integrals too.
Here is a multidimensional integration with "TrapezoidalRule". The exact result is .
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 n1) that eliminates the same order terms of truncation approximation errors of
T (f, n) and
T (f, 2 n1).
From the EulerMaclaurin formula [
DavRab84] we have
The
h^{2} 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, .
Out[11]=  

Here is an integration with a trapezoidal rule that does not use Romberg quadrature.
Out[10]=  

"TrapezoidalRule" Sampling Points and Weights
The following calculates the trapezoidal sampling points, weights, and error weights for a given precision.
Out[4]=  

Here is how the Romberg quadrature weights and error weights can be derived.
Out[9]=  

"NewtonCotesRule"
NewtonCotes integration formulas are formulas of interpolatory type with sampling points that are equally spaced.
  
"Points"  3  number of coarse NewtonCotes points 
"Type"  Closed  type of the NewtonCotes rule 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic preprocessing 
"NewtonCotesRule" options.
Let the interval of integration,
[a, b], be divided into
n1 subintervals of equal length by the points
Then the integration formula of interpolatory type is given by
When
n is large, the NewtonCotes
npoint coefficients are large and are of mixed sign.
Out[21]=  
Since this may lead to large losses of significance by cancellation, a highorder NewtonCotes rule must be used with caution.
"NewtonCotesRule" Sampling Points and Weights
The following calculates the NewtonCotes sampling points, weights, and error weights for a given precision.
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".
Out[24]=  

  
"Points"  Automatic  number of Gauss points 
"SymbolicProcessing"  Automatic  number 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
2n1 (i.e.,
if
f (x) is a polynomial of degree
≤2n1).
Gaussian rules are of open type since the integrand is not evaluated at the end points of the interval. (
Lobatto rules,
ClenshawCurtis 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
x_{1}, x_{2}, ..., x_{2n+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".
Out[25]=  

"GaussBerntsenEspelidRule" Sampling Points and Weights
The following calculates the Gaussian abscissas, weights, and BernsenEspelid error weights for a given number of coarse points and precision.
Out[27]=  

The BerntsenEspelid error weights are implemented below.
This implements the divided differences. 
This computes the abscissas and the weights of G (f, 2n+1). 
This computes the BerntsenEspelid error weights.
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 higherorder rule that reuses the Gaussian rule integrand evaluations.
  
"Points"  Automatic  number of Gauss points that will be extended with Kronrod points 
"SymbolicProcessing"  Automatic  number 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
2n1, that is,
if
f (x) is a polynomial of degree
≤2n1.
GaussKronrod 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".
Out[33]=  

For an implementation description of Kronrod extensions of Gaussian rules, see [
PiesBrand74].
"GaussKronrodRule" Sampling Points and Weights
The following calculates the GaussKronrod abscissas, weights, and error weights for a given number of coarse points and precision.
Out[35]=  

The calculations below demonstrate the degree of the GaussKronrod integration rule (see
above).
This computes the degree of the GaussKronrod integration rule.
Out[36]=  

The command below implements the integration rule weighted sums for the integral estimate,
w_{i} f (x_{i}), and the error estimate,
e_{i} f (x_{i}), where
are the abscissas,
are the weights, and
are the error weights.
These are the integral and error estimates for computed with the rule.
Out[38]=  

The integral estimate coincides with the exact result.
Out[39]=  

The error estimate is not zero since the embedded Gauss rule is exact for polynomials of degree
≤2 n1. If we integrate a polynomial of that degree, the error estimate becomes zero.
These are the integral and error estimates for computed with the rule.
Out[41]=  

"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 higherorder 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".
Out[43]=  

  
"Points"  5  number of GaussLobatto points that will be extended with Kronrod points 
"SymbolicProcessing"  Automatic  number 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
2n3, (i.e.,
if
f (x) is a polynomial of degree
≤2n3).
The Kronrod extension
LK (f, n) of a Lobatto rule with
n points
L (f, n) adds
n1 points to
L (f, n) and the extended rule is exact for polynomials of degree
3n2 if
n is even, or
3n3 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".
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 LobattoKronrod abscissas, weights, and error weights for a given number of coarse points and precision.
Out[46]=  

The calculations below demonstrate the degree of the LobattoKronrod integration rule (see
above).
This computes the degree of the LobattoKronrod integration rule.
Out[47]=  

The command below implements the integration rule weighted sums for the integral estimate,
w_{i} f (x_{i}), and the error estimate,
e_{i} f (x_{i}), where
are the abscissas,
are the weights, and
are the error weights.
These are the integral and error estimates for computed with the rule.
Out[49]=  

The preceding integral estimate coincides with the exact result.
Out[50]=  

The preceding error estimate is not zero since the embedded Lobatto rule is exact for polynomials of degree
≤2 n3. If we integrate a polynomial of that degree, the error estimate becomes zero.
These are the integral and error estimates for computed with the rule.
Out[52]=  

"ClenshawCurtisRule"
A ClenshawCurtis rule uses sampling points derived from the Chebyshev polynomial approximation of the integrand.
  
"Points"  5  number of coarse ClenshawCurtis points 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic processing 
"ClenshawCurtisRule" options.
Theoretically a ClenshawCurtis rule with
n sampling points is exact for polynomials of degree
n or less. In practice, though, ClenshawCurtis rules achieve the accuracy of the
Gaussian rules [
Evans93][
OHaraSmith68]. The error of the ClenshawCurtis formula is analyzed in [
OHaraSmith68].
The sampling points of the classical ClenshawCurtis rule are zeros of Chebyshev polynomials. The sampling points of a practical ClenshawCurtis rule are chosen to be Chebyshev polynomial extremum points. The classical ClenshawCurtis rules are not progressive, but the practical ClenshawCurtis rules are [
DavRab84][
KrUeb98].
Let
PCC (f, n) denote a practical ClenshawCurtis 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 n1). Hence the difference
PCC (f, 2 n1)PCC (f, n) can be taken to be an error estimate of the integral estimate
PCC (f, 2 n1), and can be computed without extra integrand evaluations.
The NIntegrate option Method>{"ClenshawCurtisRule", "Points">k} uses a practical ClenshawCurtis rule with 2n1 points PCC (f, 2 n1).
Out[55]=  

This example shows the number of sampling points used by NIntegrate with various values of "ClenshawCurtisRule"'s option "Points".
Out[56]=  

"ClenshawCurtisRule" Sampling Points and Weights
Here are the sampling points and the weights of the ClenshawCurtis rule for a given coarse number of points and precision.
Out[58]=  

Here is another way to compute the sampling points of PCC (f, 2 n1).
Out[60]=  

These are the integral and error estimates for computed with the rule.
Out[62]=  

"MultiPanelRule"
"MultiPanelRule" combines into one rule the applications of a onedimensional 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".
Out[64]=  

  
Method  "NewtonCotesRule"  integration rule specification that provides the abscissas, weights, and error weights for a single panel 
"Panels"  5  number of panels 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic processing 
"MultiPanelRule" options.
Let the unit interval
[0, 1] be partitioned into
k subintervals with the points
0=y_{0}<y_{1}<...<y_{k}=1.
it can be transformed into a rule for the interval
[y_{j1}, y_{j}],
Let
x_{i j}=x_{i} (y_{j}y_{j1})+y_{j1}, and
y_{j}y_{j1}=1/k,
j =1, ..., k. Then the
kpanel 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
x_{n j1}=x_{1 j}, and the number of sampling points of
k×R (f) can be reduced to
k (n1)+1. (This is done in the implementation of
"MultiPanelRule".)
More about the theory of multipanel 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.
Out[66]=  

Here are the abscissas and weights of a GaussKronrod rule.
Out[67]=  

This shows how to derive the multipanel rule weights from the original weights.
Out[69]=  

"CartesianRule"
A
ddimensional Cartesian rule has sampling points that are a Cartesian product of the sampling points of
d onedimensional rules. The weight associated with a Cartesian rule sampling point is the product of the onedimensional rule weights that correspond to its coordinates.
The Cartesian product integration for NIntegrate can be specified with the Method option value "CartesianRule".
Out[70]=  

  
Method  "GaussKronrodRule"  a rule or a list of rules with which the Cartesian product rule will be formed 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic preprocessing 
"CartesianRule" options.
For example, suppose we have the formulas:
that are exact for polynomials of degree
d_{1},
d_{2}, and
d_{3}, respectively. Then it is not difficult to see that the formula with
n_{1}×n_{2}×n_{3} points,
is exact for polynomials in
x_{1}, x_{2}, x_{3} of degree
min (d_{1}, d_{2}, d_{3}). Note that the weight associated with the abscissa
is
.
The general Cartesian product formula for
D onedimensional rules the
i of which has
n_{i} sampling points
and weights
is
Clearly (
2) can be written as
where
n=n_{k} , 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.
Out[72]=  

Cartesian rules are applicable for relatively low dimensions (
≤ 4), since for higher dimensions they are subject to "combinatorial explosion." For example, a
5dimensional Cartesian product of
5 identical onedimensional 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 onedimensional rule or a list of onedimensional rules.
Here is an example specifying Cartesian product rule integration with GaussKronrodRule.
Out[73]=  

Here is an example specifying Cartesian product rule integration with a list of onedimensional integration rules.
Out[74]=  

Another example specifying Cartesian product rule integration with a list of onedimensional integration rules.
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.
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.
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".
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,
{x_{1}, x_{2}, ..., x_{d}}, for the coordinates of which the inequality
x_{1}≥x_{2}≥...≥ x_{d} holds, is called a generator. (See [
KrUeb98][
GenzMalik83].)
  
"Generators"  5  number of generators of the fully symmetric rule 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic preprocessing 
"MultiDimensionalRule" options.
If an integration rule has
K orbits denoted
_{1}, _{2}, ..., _{K}, and the
i^{th} of them,
_{i}, has a weight
w_{i} 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.
NIntegrate's
"MultiDimensionalRule" object is basically an interface to three different integration rule objects that combine an integration rule and one or several null rules. Their number of generators and orders are summarized in the table below. The rule objects with 6 and 9 generators use three null rules, each of which is a linear combination of two null rules. The null rule linear combinations are used in order to avoid
phase errors. See [
BerntEspGenz91] for more details about how the null rules are used.
Number of generators and orders of NIntegrate's fully symmetric rules:
This is the number of sampling points used by NIntegrate with its fully symmetric multidimensional integration rules for integrals of the form , m=1, ..., 20.
Out[81]=  

"MultiDimensionalRule" Sampling Points and Weights
This subsection gives an example of a calculation of an integral estimate with a fully symmetric multidimensional rule.
Here is the parameter for the number of generators. 
This function takes a generator point and creates its orbit. 
The generators and weights for given number of generators. 
This computes the orbit of each generator. 
This makes graphics primitives for points of the orbits. 
Here is how the different orbits look.
Out[95]=  

Here are all rule points together.
Out[96]=  

"MonteCarloRule"
A Monte Carlo rule estimates an integral by forming a uniformly weighted sum of integrand evaluations over random (quasirandom) sampling points.
Here is an example of using "MonteCarloRule" with 1000 sampling points.
Out[97]=  

  
"Points"  100  number of sampling points 
"PointGenerator"  Random  sampling points coordinates generator 
"AxisSelector"  Automatic  selection algorithm of the splitting axis when global adaptive Monte Carlo integration is used 
"SymbolicProcessing"  Automatic  number of seconds to do symbolic preprocessing 
"MonteCarloRule" options.
In Monte Carlo methods [
KrUeb98], the
ddimensional integral
_{V}f (x)x is interpreted as the following expected (mean) value
where
E (f) is the mean value of the function
f interpreted as a random variable, with respect to the uniform distribution on
V, that is, the distribution with probability density
vol (V)^{1}Boole (xV).
Boole (xV) denotes the characteristic function of the region
V;
vol (V) denotes the volume of
V.
The crude Monte Carlo estimate of the expected value
E (f) is obtained by taking
n independent random vectors
x_{1}, x_{2}, ..., x_{n}^{d} with density
vol (V)^{1}Boole (xV) (that is, the vectors are uniformly distributed on
V), and making the estimate
Remark: The function
vol (V)^{1}Boole (xV) is a valid probability density function because it is nonnegative on the whole of
^{d} and
_{d } vol (V)^{1}Boole (xV)x=1.
According to the strong law of large numbers, the convergence
happens with probability
1. The strong law of large numbers does not provide information for the error
MC (f, n)_{V}f (x)x, so a probabilistic estimate is used.
Let
be defined as
Formula (
5) is an unbiased estimator of
(that is, the expectation of
MC (f, n) for various sets of
is
) and its variance is
where
Var (f) denotes the variance of
f, The standard error of
MC (f, n) is thus
.
In practice the
Var (f) is not known, so it is estimated with the formula
The standard error of
MC (f, n) is then
The result of the Monte Carlo estimation can be written as
MC (f, n)±SD (f, n).
It can be seen from Equation (
6) that the convergence rate of the crude Monte Carlo estimation does not depend on the dimension
d of the integral, and if
n sampling points are used then the convergence rate is
.
The
NIntegrate integration rule
"MonteCarloRule" calculates the estimates
MC (f, n) and
SD (f, n).
The estimates can be improved incrementally. That is, if we have the estimates
MC (f, n_{0}) and
SD (f, n_{0}), and a new additional set of sample function values
{f_{1}, f_{2}, ..., f_{n1}}, then using (
7) and (
8) we have
To compute the estimates
MC (f, n_{0}+n_{1}) and
SD (f, n_{0}+n_{1}), it is not necessary to know the random points used to compute the estimates
MC (f, n_{0}) and
SD (f, n_{0}).
"AxisSelector"
When used for multidimensional global adaptive integration,
"MonteCarloRule" chooses the splitting axis of an integration subregion it is applied to in two ways: (i) by random selection or (ii) by minimizing the sum of the variances of the integral estimates of each half of the subregion, if the subregion is divided along that axis. The splitting axis is selected after the integral estimation.
The random axis selection is done in the following way.
"MonteCarloRule" keeps a set of axes for selection,
A. Initially
A contains all axes. An element of
A is randomly selected. The selected axis is excluded from
A. After the next integral estimation, an axis is selected from
A and excluded from it, and so forth. If
A is empty, it is filled up with all axes.
The minimization of variance axis selection is done in the following way. During the integration over the region, a subset of the sampling points and their integrand values is stored. Then for each axis, the variances of the two subregions that the splitting along this axis will produce are estimated using the stored sampling point and corresponding integrand values. The axis for which the sum of these variances is minimal is chosen to be the splitting axis, since this would mean that if the region is split on that axis, the new integration error estimate will be minimal. If it happens that for some axis all stored points are clustered in one of the halfregions, then that axis is selected for splitting.
 
Random  random splitting axis election 
MinVariance{MinVariance, "SubsampleFraction">frac}  splitting axis selection that minimizes the sum of variances of the new regions 
"AxisSelector" options.
  
"SubsampleFraction"  1/10  fraction of the sampling points used to determine the splitting axis 
MinVariance options.
This is an example of using "MonteCarloRule"'s option "AxisSelector".
Out[98]=  

In the examples below the two axis selection algorithms are compared. In general, the minimization of variance selection uses less number of sampling points. Nevertheless, using the minimization of variance axis selection slows down the application of
"MonteCarloRule". So for integrals for which both axis selection methods would result in the same number of sampling points, it is faster to use random axis selection. Also, using larger fraction sampling points to determine the splitting axis in minimization of variance selection makes the integration slower.
Consider the following function.
Out[3]=  

These are the adaptive Monte Carlo integration sampling points for the function above with random choice of splitting axis.
Out[5]=  

These are the sampling points with choice of splitting axes that minimize the variance. Compared to the previous Monte Carlo integration, the sampling points of this one are more concentrated around the circle (x1/2)^{2}+ (y1/2)^{2}=1/6, and their number is nearly twice as small.
Out[6]=  

Here is an adaptive Monte Carlo integration that uses random axis selection.
Out[104]=  

Here is an adaptive Monte Carlo integration for the preceding integral that uses the minimization of variance axis selection and is slower than using random axis selection.
Out[105]=  

Using a larger fraction of stored points for the minimization of variance axis choice slows down the integration.
Out[106]=  

Comparisons of the Rules
All integration rules, except
"MonteCarloRule", are to be used by
adaptive strategies in
NIntegrate. Changing the type and the number of points of the integration rule component for an integration strategy will make a different integration algorithm. In general these different integration algorithms will perform differently for different integrals. Naturally the following questions arise.
1. Is there a type of rule that is better than other types for any integral or for integrals of a certain type?
2. Given an integration strategy, what rules perform better with it? For what integrals?
3. Given an integral, an integration strategy, and an integration rule, what number of points in the rule will minimize the total number of sampling points used to reach an integral estimate that satisfies the precision goal?
For a given integral and integration strategy the integration rule which achieves a result that satisfies the precision goal with the smallest number of sampling points is called the best integration rule. There are several factors that determine the best integration rule.
4. In general the higher the degree of the rule the faster the integration will be for smooth integrands and for higherprecision goals. On the other hand, the rule degree might be too high for the integrand and hence too many sampling points might be used when the adaptive strategies work around, for example, the integrand's discontinuities.
5. The error estimation functional of a rule influences significantly the total amount of work by the integration strategy. Rules with a smaller number of points might lead (i) to a wrong result because of underestimation of the integral, or (ii) to applying too many sampling points because of overestimation of the integrand. (See "
Examples of Pathological Behavior".) Further, the error estimation functional might be computed with one or several embedded
null rules. In general, the larger the number of the null rules the better the error estimation—fewer
phase errors are expected. The number of the null rules and the weights assigned to them in the sum that computes the error estimate determines the sets of pathological integrals and integrals hard to compute for that rule. (Some of
NIntegrate's multidimensional rules use several embedded null rules to compute the error estimate. All of
NIntegrate's onedimensional integration rules use only one null rule.)
6. Local adaptive strategies are more effective with closed rules that have their sampling points more uniformly distributed (for example,
"ClenshawCurtisRule") than with open rules (for example,
GaussKronrodRule) and closed rules that have sampling points distributed in a nonuniform way (for example,
"LobattoKronrodRule").
7. The percent of points reused by the strategy might greatly determine what is the best rule. For onedimensional integrals,
"LocalAdaptive" reuses all points of the closed rules.
"GlobalAdaptive" throws away almost all points of the regions that need improvement of their error estimate.
Number of points in a rule
This subsection demonstrates with examples that the higher the degree of the rule the faster the integration will be for smooth integrands and for higherprecision goals. It also shows examples in which the degree of the rule is too high for the integrand and hence too many sampling points are used when the adaptive strategies work around the integrand's discontinuities. All examples use Gaussian rules with BerntsenEspelid error estimate.
Here is the error of a Gaussian rule in the interval
[a, b].
Here is a function that calculates the error of a rule for the integral , using the exact value computed by Integrate for comparison. 
This defines a list of functions. 
Here are plots of the functions in the interval [0, 1].
Out[109]=  

Here is the computation of the errors of "GaussBerntsenEspelidRule" for , , , and for a range of points. 
Here are plots of how the logarithm of the error decreases for each of the functions. It can be seen that the integral estimates of discontinuous functions and functions with discontinuous derivatives improve slowly when the number of points is increased.
Out[117]=  

Minimal Number of Sampling Points
Here is a function that finds the number of sampling points used in an integration. 
This finds the number of sampling points used for a range of precision goals and a range of integration rule coarse points. 
This finds the for each precision the minimum total number of sampling points. This way the number of coarse integration rule points used is also found. 
This is a plot of the precision goal and the number of integration rule points with which the minimum number of total sampling points was used.
Out[128]=  

Rule Comparison
Here is a function that calculates the error of a rule for the integral , using the exact value computed by Integrate for comparison. 
This defines a list of functions. 
Here are plots of the functions in the interval [0, 1].
Out[131]=  

This is the computation of the errors of "GaussKronrodRule", "LobattoKronrodRule", "TrapezoidalRule", and "ClenshawCurtisRule" for each of the integrals , , , and for a range of points. 
Here are plots of how the logarithms of the errors decrease for each rule and each function.
Out[136]=  

Examples of Pathological Behavior
Tricking the Error Estimator
In this subsection an integral will be discussed which
NIntegrate underestimates with its default settings since it fails to detect part of the integrand. The part is detected if the precision goal is increased.
The Wrong Estimation
Consider the following function. 
Here is its exact integral over [0, 1].
Out[138]=  

This is too inaccurate when compared to the exact value.
Out[140]=  

Here is the plot of the function, which is also wrong.
Out[141]=  

Better Results
This is a table that finds the precision goal for which no good results are computed.
Out[18]=  

If the plot points are increased, the plot of the function looks different.
Out[144]=  

Here is the zoomed plot of the spike that Plot is missing with the default options.
Out[145]=  

If this part of the function is integrated, the result fits the quantity that is "lost" (or "missed") by NIntegrate with the default option settings.
Out[146]=  
Out[147]=  

Why the Estimator Is Misled
These are the abscissas and weights of a GaussKronrod rule used by default by NIntegrate. 
This defines a function for application of the rule. 
This finds the points at which the adaptive strategy samples the integrand. 
This is a plot of the sampling points. The vertical axis is for the order at which the points have been used to evaluate the integrand.
Out[149]=  

It can be seen on the preceding plot that
NIntegrate does extensive computation around the top of the second spike near
x=0.4.
NIntegrate does not do as much computation around the unintegrated spike near
x=0.6.
These are GaussKronrod and Gauss abscissas in the last set of sampling points, which is over the region [0.5, 0.75].
Out[150]=  
Out[151]=  

Here the integrand is applied over the abscissas. 
Here is a polynomial approximation of the integrand over the abscissas. 
These plots show that the two polynomial approximations almost coincide over x=0.6.
Out[156]=  
Out[158]=  

If the polynomials are integrated over the region where 0.6 is placed, the difference between them, which NIntegrate uses as an error estimate, is really small. 
Since the difference is the error estimate assigned for the region
[0.5, 0.75], with the default precision goal
NIntegrate never picks it up for further integration refinement.
Phase Errors
In this subsection are discussed causes why integration rules might seriously underestimate or overestimate the actual error of their integral estimates. Similar discussion is given in [
LynKag76].
Consider the numerical and symbolic evaluations of the integral of f[x, 0.415, 1.25] over the region [1, 1].
Out[163]=  
Out[164]=  

They differ significantly. The precision goal requested is 2, but relative error is much higher than 10^{2}.
Out[165]=  

(Note that
NIntegrate gives correct results for higherprecision goals.)
Below is an explanation why this happens.
Let the integration rule
R_{2} be embedded in the rule
R_{1}. Accidentally, the error estimate
R_{1}^{V}[f]R_{2}^{V}[f] of the integral estimate
R_{1}^{V}[f], where
V=[1, 1], can be too small compared with the exact error
R_{2}^{V}[f]_{V}f (x)x.
To demonstrate this, consider the GaussKronrod rule GK[f, 5] with 11 sampling points that has an embedded Gauss rule G[f, 5] with 5 sampling points. (This is the rule used in the two integrations above.) 
This defines a function that applies the rule. 
This is the integral of f[x, , ] previously defined. 
We can plot a graph with the estimated error of GK (f, 5) and the real error for different values of in [1, 1]. That is, we plot GK (f, 5)G (f, 5) and .
Out[169]=  

In the plot above the blue graph is for the estimated error,
GK (f, 5)G (f, 5). The graph of the actual error
is red.
We can see that the value
0.415 of the parameter
is very close to one of the
GK (f, 5)G (f, 5) local minimals.
A onedimensional quadrature rule can be seen as the result of the integration of a polynomial that is fitted through the rule's abscissas and the integrand values over them. We can further try to see the actual fitting polynomials for the integration of
f[x, , ].
In the plots below the function
f[x, , ] is plotted in red, the Gauss polynomial is plotted in blue, the GaussKronrod polynomial is plotted in violet, the Gauss sampling points are in black, and the GaussKronrod sampling points are in red.
We can see that since the peak of
f[x, 0.415, 1.25] falls approximately halfway between two abscissas, its approximation is an underestimate.
Out[172]=  
Conversely, we can see that since the peak of
f[x, 0.53, 1.25] falls approximately on one of the abscissas, its approximation is an overestimate.
Out[173]=  
Index of Technical Terms