NIntegrate Integration Strategies
Introduction
An integration strategy is an algorithm that attempts to compute integral estimates that satisfy user-specified precision or accuracy goals.
An integration strategy prescribes how to manage and create new elements of a set of disjoint subregions of the initial integral region. Each subregion might have its own integrand and integration rule associated with it. The integral estimate is the sum of the integral estimates of all subregions. Integration strategies use integration rules to compute the subregion integral estimates. An integration rule samples the integrand with a set of points, called abscissas (or sampling points).
To improve an integral estimate the integrand should be sampled at additional points. There are two principal approaches: (i) adaptive strategies try to identify the problematic integration areas and concentrate the computational efforts (i.e., sampling points) on them; (ii) non-adaptive strategies increase the number of sampling points over the whole region in order to compute a higher degree integration rule estimate that reuses the integrand evaluations of the former integral estimate.
Both approaches can use symbolic preprocessing and variable transformation or sequence summation acceleration to achieve faster convergence.
In the following integration, the symbolic piecewise preprocessor in NIntegrate recognizes the integrand as a piecewise function, and the integration is done over regions for which x≥1 with the integrand  and regions for which x≤1 with the integrand  .
| Out[31]= |  |
|
Here is a plot of all sampling points used in the integration. The integrand is sampled at the x coordinates in the order of the y coordinates (in the plot). It can be seen that the sampling points are concentrated near the singularity point 1. The patterns formed by the sampling points at the upper part of the plot differ from the patterns of the lower part of the plot because a singularity handler is applied.
| Out[11]= |  |
|
The section
"Adaptive Strategies" gives a general description of the adaptive strategies. The default (main) strategy of
NIntegrate is global adaptive, which is explained in the section
"Global Adaptive Strategy". Complementary to it is the local adaptive strategy, which is explained in the section
"Local Adaptive Strategy". Both adaptive strategies use singularity handling mechanisms, which are explained in the section
"Singularity Handling".
The strategies
NIntegrate uses for special types of integrals (or integrands) are explained in the corresponding sections:
"Duffy's coordinates strategy",
"Oscillatory strategies", and
"Cauchy principal value integration".
Here is a table with links to descriptions of built-in integration strategies of
NIntegrate.
Adaptive Strategies
Adaptive strategies try to concentrate computational efforts where the integrand is discontinuous or has some other kind of singularity. Adaptive strategies differ by the way they partition the integration region into disjoint subregions. The integral estimates of each subregion contribute to the total integral estimate.
The basic assumption for the adaptive strategies is that for given integration rule
R and integrand
f, if an integration region
V is partitioned into, say, two disjoint subregions
V1 and
V2 ,
V=V1
V2,
V1
V2=0, then the sum of the integral estimates of
R over
V1 and
V2 is closer to the actual integral
Vf
x. In other words,
and (
1) will imply that the sum of the error estimates for
RV1 (f) and
RV2 (f) is smaller than the error estimate of
RV (f).
Hence an adaptive strategy has these components [
MalcSimp75]:
(i) an integration rule to compute the integral and error estimates over a region;
(ii) a method for deciding which elements of a set of regions
to partition/subdivide;
(iii) stopping criteria for deciding when to terminate the adaptive strategy algorithm.
Global Adaptive Strategy
A global adaptive strategy reaches the required precision and accuracy goals of the integral estimate by recursive bisection of the subregion with the largest error estimate into two halves, and computes integral and error estimates for each half.
The global adaptive algorithm for NIntegrate is specified with the Method option value "GlobalAdaptive".
| Out[32]= |  |
|
| | |
| Method | Automatic | integration rule used to compute integral and error estimates over each subregion |
| "SingularityDepth" | Automatic | number of recursive bisections before applying a singularity handler |
| "SingularityHandler" | Automatic | singularity handler |
| "SymbolicProcessing" | Automatic | number of seconds to do symbolic preprocessing |
"GlobalAdaptive" options.
"GlobalAdaptive" is the default integration strategy of
NIntegrate. It is used for both one-dimensional and multidimensional integration.
"GlobalAdaptive" works with both
Cartesian product rules and
fully symmetric multidimensional
rules.
"GlobalAdaptive" uses a data structure called a "heap" to keep the set of regions partially sorted, with the largest error region being at the top of the heap. In the main loop of the algorithm the largest error region is bisected in the dimension that is estimated to be responsible for most of its error.
It can be said that the algorithm produces the leaves of a binary tree, the nodes of which are the regions. The children of a node/region are its subregions obtained after bisection.
After a bisection of a region and the subsequent integration over the new (sub)regions, new global integral and global error estimates are computed, which are sums of the integral and error estimates of all regions that are leaves of the binary tree.
Each region has a record of how many bisections are made per dimension in order to produce it. When a region has been produced through too many bisections a singularity flattening algorithm is applied to it; see
Singularity Handling.
"GlobalAdaptive" stops if the following expression is true:
where
pg and
ag are precision and accuracy goals.
The strategy also stops when the number of recursive bisections of a region exceeds a certain number (see
MinRecursion and MaxRecursion), or when the global integration error oscillates too much (see
"MaxErrorIncreases").
Theoretical and practical evidence show that the global adaptive strategies have in general better performance than the local adaptive strategies [
MalcSimp75][
KrUeb98].
MinRecursion and MaxRecursion
The minimal and maximal depths of the recursive bisections are given by the values of the options
MinRecursion and
MaxRecursion.
If for any subregion the number of bisections in any of the dimensions is greater than
MaxRecursion then the integration by
"GlobalAdaptive" stops.
Setting
MinRecursion to a positive integer forces recursive bisection of the integration regions before the integrand is ever evaluated. This can be done to ensure that a narrow spike in the integrand is not missed. (See
Tricking the error estimator.)
For multidimensional integration an effort is made to bisect in each dimension for each level of recursion in
MinRecursion.
"MaxErrorIncreases"
Since (
2) is expected to hold in
"GlobalAdaptive" the global error is expected to decrease after the bisection of the largest error region and the integration over its new parts. In other words the global error is expected to be more or less monotonically decreasing with respect to the number of integration steps.
The global error might oscillate due to
phase errors of the integration rules. Still, the global error is assumed at some point to start decreasing monotonically.
Below are listed cases in which this assumption might become false.
(i) The actual integral is zero.
| Out[3]= |  |
|
(ii) The specified working precision is not dense enough for the specified precision goal.
The working precision is not dense enough.
Out[33]//InputForm= |
| |  |
|
(iii) The integration is badly conditioned [
KrUeb98]. For example, the reason might be that the integrand is defined by complicated expressions or in terms of approximate solutions of mathematical problems (such as differential equations or nonlinear algebraic equations).
The strategy
"GlobalAdaptive" keeps track of the number of times the total error estimate has not decreased after the bisection of the region with the largest error estimate. When that number becomes bigger than the value of the
"GlobalAdaptive" option
"MaxErrorIncreases", the integration stops with a message (
NIntegrate::eincr).
The default value of
"MaxErrorIncreases" is 400 for one-dimensional integrals and 2000 for multidimensional integrals.
The following integration invokes the message NIntegrate::eincr, with the default value of "MaxErrorIncreases".
| Out[1]= |  |
|
Increasing "MaxErrorIncreases" silences the NIntegrate::eincr message.
| Out[2]= |  |
|
The result compares well with the exact value.
| Out[4]= |  |
|
Example Implementation of a Global Adaptive Strategy
This computes Gauss-Kronrod abscissas, weights, and error weights. |
This is a definition of a function that applies the integration rule with abscissas and weights computed to the function f over the interval {a, b}. |
This is a definition of a simple global adaptive algorithm that finds the integral of the function f over the interval {aArg, bArg} with relative error tol. |
This defines an integrand. |
The global adaptive strategy defined earlier gives the following result.
| Out[19]= |  |
|
Here is the exact result.
| Out[20]= |  |
|
The relative error is within the prescribed tolerance.
| Out[21]= |  |
|
Local Adaptive Strategy
In order to reach the required precision and accuracy goals of the integral estimate, a local adaptive strategy recursively partitions the subregion into smaller disjoint subregions and computes integral and error estimates for each of them.
The local adaptive algorithm for NIntegrate is specified with the Method option value "LocalAdaptive".
| Out[5]= |  |
|
| | |
| Method | Automatic | integration rule used to compute integral and error estimates over the subregions |
| "SingularityDepth" | Automatic | number of recursive bisections before applying a singularity handler |
| "SingularityHandler" | Automatic | singularity handler |
| "Partitioning" | Automatic | how to partition the regions in order to improve their integral estimate |
| "InitialEstimateRelaxation" | True | attempt to adjust the magnitude of the initial integral estimate in order to avoid unnecessary computation |
| "SymbolicProcessing" | Automatic | number of seconds to do symbolic preprocessing |
"LocalAdaptive" options.
Like
"GlobalAdaptive",
"LocalAdaptive" can be used for both one-dimensional and multidimensional integration.
"LocalAdaptive" works with both
Cartesian product rules and
fully symmetric multidimensional
rules.
The
"LocalAdaptive" strategy has an initialization routine and a Recursive Routine (RR). RR produces the leaves of a tree, the nodes of which are regions. The children of a node/region are subregions obtained by its partition. RR takes a region as an argument and returns an integral estimate for it.
RR uses an integration rule to compute integral and error estimates of the region argument. If the error estimate is too big, RR calls itself on the region's disjoint subregions obtained by partition. The sum of the integral estimates returned from these recursive calls becomes the region's integral estimate.
RR makes the decision to continue the recursion knowing only the integral and error estimates of the region at which it is executed. (This is why the strategy is called "local adaptive.")
The initialization routine computes an initial estimation of the integral over the initial regions. This initial integral estimate is used in the stopping criteria of RR: if the error of a region is significant compared to the initial integral estimate then that region is partitioned into disjoint regions and RR is called on them; if the error is insignificant the recursion stops.
The error estimate of a region,
regionError, is considered insignificant if
The stopping criteria (
3) will compute the integral to the working precision. Since you want to compute the integral estimate to user-specified precision and accuracy goals, the following stopping criteria is used:
where
eps is the smallest number such that
1+eps≠1 at the working precision, and
pg and
ag are the user-specified precision and accuracy goals.
The recursive routine of
"LocalAdaptive" stops the recursion if:
1. there are no numbers of the specified working precision between region's boundaries;
2. the maximum recursion level is reached;
3. the error of the region is insignificant, i.e., the criteria (
4) is true.
MinRecursion and MaxRecursion
The options
MinRecursion and
MaxRecursion for
"LocalAdaptive" have the same meaning and functionality as they do for
"GlobalAdaptive". See
MinRecursion and MaxRecursion.
"InitialEstimateRelaxation"
After the first recursion is finished a better integral estimate,
I2, will be available. That better estimate is compared to the two integral estimates,
I1 and
I1e that the integration rule has used to give the integral estimate (
I1) and the error estimate (
|I1-I1e|) for the initial step. If
then the integral estimate
integralEst in (
5) can be increased—that is, the condition (
6) is relaxed—with the formula
since
<1 means that the rule's integral estimate is more accurate than what the rule's error estimate predicts.
"Partitioning"
"LocalAdaptive" has the option
"Partitioning" to specify how to partition the regions that do not satisfy (
7). For one-dimensional integrals, if
"Partitioning" is set to
Automatic,
"LocalAdaptive" partitions a region between the sampling points of the (rescaled) integration rule. In this way, if the integration rule is of
closed type, every integration value can be reused. If
"Partitioning" is given a list of integers
{p1, p2, ..., pn} with length
n that equals the number of integral variables, each dimension
i of the integration region is divided into
pi equal parts. If
"Partitioning" is given an integer
p, all dimensions are divided into
p equal parts.
Consider the following function.
| Out[4]= |  |
|
These are the sampling points used by "LocalAdaptive" with its automatic region partitioning. It can be seen that the sampling points of each recursion level are between the sampling points of the previous recursion level.
| Out[2]= |  |
|
These are the sampling points used by "LocalAdaptive" integration which partitions the regions with large error into three subregions. The patterns formed clearly show the three next recursion level subregions of each region of the first and second recursion levels.
| Out[6]= |  |
|
Multidimensional example of using the "Partitioning" option. To make the plot, the sampling points of the first region to be integrated, [0, 1]×[0, 1], are removed.
| Out[10]= |  |
|
Reuse of Integrand Values
With its default partitioning settings for one-dimensional integrals
"LocalAdaptive" reuses the integrand values at the end points of the sub-intervals that have integral and error estimates that do not satisfy (
8).
Sampling points of the integration of  by "LocalAdaptive". The variable rulePoints determines the number of points in the integration rule used by "LocalAdaptive".
| Out[15]= |  |
| Out[16]= |  |
|
The percent of reused points in the integration.
| Out[19]= |  |
|
Example Implementation of a Local Adaptive Strategy
This computes Clenshaw-Curtis abscissas, weights, and error weights. |
This is a definition of a function that applies the integration rule, with the abscissas and weights computed in the previous example, to the function f over the interval {a, b}. |
This defines a simple local adaptive algorithm that finds the integral of the function f over the interval {aArg, bArg} with relative error tol. |
The local adaptive strategy gives the result.
| Out[38]= |  |
|
This is the exact result.
| Out[39]= |  |
|
The relative error is within the prescribed tolerance.
| Out[40]= |  |
|
"GlobalAdaptive" versus "LocalAdaptive"
In general the global adaptive strategy has better performance than the local adaptive one. In some cases though the local adaptive strategy is more robust and/or gives better performance.
There are two main differences between
"GlobalAdaptive" and
"LocalAdaptive":
1.
"GlobalAdaptive" stops when the sum of the errors of all regions satisfies the precision goal, while
"LocalAdaptive" stops when the error of each region is small enough compared to an estimate of the integral.
2. To improve the integral estimate
"GlobalAdaptive" bisects the region with largest error, while
"LocalAdaptive" partitions all regions the error for which is not small enough.
For multidimensional integrals
"GlobalAdaptive" is much faster because
"LocalAdaptive" does partitioning along each axis, so the number of regions can explode combinatorically.
Why and how global adaptive strategy is faster for one-dimensional smooth integrands is proved (and explained) in [
MalcSimp75].
When
"LocalAdaptive" is faster and performs better than
"GlobalAdaptive", it is because the precision-goal-stopping criteria and partitioning strategy of
"LocalAdaptive" are more suited for the integrand's nature. Another factor is the ability of
"LocalAdaptive" to reuse the integral values of all points already sampled.
"GlobalAdaptive" has the ability to reuse very few integral values (at most 3 per rule application, 0 for the default one-dimensional rule, the Gauss-Kronrod rule).
The following subsection demonstrates the performance differences between
"GlobalAdaptive" and
"LocalAdaptive".
"GlobalAdaptive" Is Generally Better than "LocalAdaptive"
The table that follows, with timing ratios and numbers of integrand evaluations, demonstrates that
"GlobalAdaptive" is better than
"LocalAdaptive" for the most common cases. All integrals considered are one-dimensional over
[0, 1], because (i) for multidimensional integrals the performance differences should be expected to deepen, since
"LocalAdaptive" partitions the regions along each axis, and (ii) one-dimensional integrals over different ranges can be rescaled to be over
[0, 1].
Here are the definitions of some functions, precision goals, number of integrations, and the integration rule. The variable integrationRule can be changed in order to compare the profiling runs with the same integration rule. The last function is derived from  -xsin (x) by the variable change  that transforms [0, 1) into [0,  ).
| Out[74]= |  |
|
Exact integral values. All integrals are over [0, 1]. |
"GlobalAdaptive" timings. |
"GlobalAdaptive" function evaluations. |
"LocalAdaptive" function evaluations. |
Table with the timing ratios and integrand evaluations.
| Out[80]= |  |
|
Table with the errors of the integrations. Both "GlobalAdaptive" and "LocalAdaptive" reach the required precision goals.
| Out[81]= |  |
|
Singularity Handling
The adaptive strategies of
NIntegrate speed up their convergence through variable transformations at the integration region boundaries and user-specified singular points or manifolds. The adaptive strategies also ignore the integrand evaluation results at singular points.
Singularity specification is discussed in
"User-specified Singularities".
Multidimensional singularity handling with variable transformations should be used with caution; see
"IMT Multidimensional Singularity Handling". Coordinate change for a multidimensional integral can simplify or eliminate singularities; see
"Duffy's Coordinates for Multidimensional Singularity Handling".
For details about how
NIntegrate ignores singularities see
"Ignoring the Singularity".
The computation of Cauchy principal value integrals is described in
"Cauchy Principal Value Integration".
User-Specified Singularities
Point Singularities
If it is known where the singularities occur, they can be specified in the ranges of integration, or through the option
Exclusions.
Here is an example of an integral that has two singular points at  and  .
| Out[58]= |  |
|
Here is an example of a two-dimensional integral with a singular point at (1, 1).
| Out[59]= |  |
|
Here is an example of an integral that has two singular points at  and  specified with the Exclusions option.
| Out[60]= |  |
|
Here is an example of a two-dimensional integral with a singular point at (1, 1) specified with the Exclusions option.
| Out[61]= |  |
|
Curve, Surface, and Hypersurface Singularities
Singularities over curves, surfaces, or hypersurfaces in general can be specified through the option
Exclusions with their equations. Such singularities, generally, cannot be specified using variable ranges.
This two-dimensional function is singular along the curve x2+y2=1.
| Out[62]= |  |
|
Using Exclusions the integral is quickly calculated.
| Out[12]= |  |
|
NIntegrate will reach convergence much more slowly if no singularity specification is given.
| Out[35]= |  |
|
Here is an example of a case in which a singular curve can be specified with the variable ranges. If x [0, 2] and y[0, 2] this would not be possible—see the following example.
| Out[10]= |  |
|
Example Implementation of Curve, Surface, and Hypersurface Singularity Handling
If the curve, surface, or hypersurface on which the singularities occur is known in implicit form (i.e., it can be described as a single equation) the function
Boole can be used to form integration regions that have the singular curves, surfaces, or hypersurfaces as boundaries.
This two-dimensional function has singular points along the curve x2+y2=1.
| Out[66]= |  |
|
Using Boole the integral is calculated quickly.
| Out[9]= |  |
|
This two-dimensional function has singular points along the curve x+ (1-y)2=1.
| Out[68]= |  |
|
Again, using Boole the integral is calculated quickly.
| Out[8]= |  |
|
This is how the sampling points of the integration look.
| Out[7]= |  |
|
Here is a function that takes a singular curve, surface, or hypersurface specification and uses the function Boole to make integration regions that have the singularities on their boundaries. |
This defines a three-dimensional function. |
Here is the integral of a three-dimensional function with singular points along the surface x+ (1-y)2+ (1-z)2=1.
| Out[3]= |  |
|
These are the sampling points of the integration.
| Out[5]= |  |
|
"SingularityHandler" and "SingularityDepth"
Adaptive strategies improve the integral estimate by region bisection. If an adaptive strategy subregion is obtained by the number of bisections specified by the option
"SingularityDepth", it is decided that subregion has a singularity. Then the integration over that subregion is done with the singularity handler specified by
"SingularityHandler".
| | |
| "SingularityDepth" | Automatic | number of recursive bisections before applying a singularity handler |
| "SingularityHandler" | Automatic | singularity handler |
"GlobalAdaptive" and "LocalAdaptive" singularity handling options.
If there is an integrable singularity at the boundary of a given region of integration, bisection could easily recur to
MaxRecursion before convergence occurs. To deal with these situations the adaptive strategies of
NIntegrate use variable transformations (
IMT,
"DoubleExponential",
SidiSin) to speed up the integration convergence, or a region transformation (
Duffy's coordinates) that relaxes the order of the singularity. The theoretical background of the variable transformation singularity handlers is given by the Euler-Maclaurin formula [
DavRab84].
Use of the IMT Variable Transformation
The IMT variable transformation is the variable transformation in a quadrature method proposed by Iri, Moriguti, and Takasawa, called in the literature the IMT rule [
DavRab84][
IriMorTak70]. The IMT rule
is based upon the idea of transforming the independent variable in such a way that all derivatives of the new integrand vanish at the end points of the integration interval. A trapezoidal rule is then applied to the new integrand, and under proper conditions high accuracy of the result might be attained [
IriMorTak70][
Mori74].
Here is a numerical integration that uses the IMT variable transformation for singularity handling.
| Out[13]= |  |
|
| | |
| "TuningParameters" | 10 | a pair of numbers {a, p} that are the tuning parameters in the IMT transformation formula  ; if only a number a is given, it is interpreted as {a, 1} |
IMT singularity handler options.
Adaptive strategies of
NIntegrate employ only the transformation of
the
IMT rule. With the decision that a region might have a singularity, the IMT transformation is applied to its integrand. The integration continues, though not with a trapezoidal rule, but with the same integration rule used before the transformation. (Singularity handling with
"DoubleExponential" switches to a trapezoidal integration rule.)
Also, adaptive strategies of
NIntegrate use a variant of the original IMT transformation, with the transformed integrand vanishing only at one of the ends.
The IMT transformation  a, p (t): (0, 1]→ (0, 1], a>0, p>0, is defined. |
The parameters
a and
p are called tuning parameters [
MurIri82].
The limit of the derivative of the IMT transformation is 0.
| Out[16]= |  |
|
Here is the plot of the IMT transformation.
| Out[17]= |  |
|
From the graph above follows that the transformed sampling points are much denser around 0. This means that if the integrand is singular at 0 it will be sampled more effectively, since a larger part of the integration rule sampling points will contribute large integrand values to the integration rule's integral estimate.
Since for any given working precision the numbers around 0 are much denser than the numbers around 1, after a region bisection the adaptive strategies of NIntegrate reverse the bisection variable of the subregion that has the right end of the bisected interval. This can be seen from the following plot.
| Out[19]= |  |
|
No other singularity handler is applied to the subregions of a region to which the IMT variable transformation has been applied.
IMT Transformation by Example
Consider the function  over (0, 1] that has a singularity at 0.
| Out[30]= |  |
|
Assume the integration is done with "GlobalAdaptive", with singularity handler IMT and singularity depth 4. After four bisections "GlobalAdaptive" will have a region with boundaries {0, 1/16} that contains the singular end point. For that region the IMT variable transformation will change its boundaries to {0, 1} and its integrand to the following.
| Out[32]= |  |
|
Here is the plot of the new integrand.
| Out[34]= |  |
|
The singularity is smashed!
Some of the sampling points, though, become too close to the singular end, and therefore special care should be taken for sampling points that coincide with the singular point because of the IMT transformation.
NIntegrate ignores evaluations at singular points; see
"Ignoring the Singularity".
The Gauss-Kronrod sampling points for the region {0, 1/16} and the derivatives of the rescaling follow.
| Out[36]= |  |
| Out[37]= |  |
|
Here is the integral estimate.
| Out[38]= |  |
|
With the IMT transformation, these are the sampling points and derivatives.
| Out[39]= |  |
| Out[40]= |  |
|
Here is the integral estimate with the IMT transformation.
| Out[41]= |  |
|
The estimate calculated with the IMT variable transformation is much closer to the exact value.
| Out[42]= |  |
|
Use of Double-Exponential Quadrature
When adaptive strategies use the IMT variable transformation they do not change the integration rule on the IMT-transformed regions. In contrast to this you can use both a variable transformation and a different integration rule on the regions considered to have singularity. (This is more in the spirit of the
IMT rule [
DavRab84].) This is exactly what happens when
double-exponential quadrature is used—double-exponential quadrature uses the trapezoidal rule.
NIntegrate can use double-exponential quadrature for singularity handling only for one-dimensional integration.
Here is a numerical integration that uses double-exponential quadrature for singularity handling.
| Out[103]= |  |
|
IMT versus "DoubleExponential" versus No Singularity Handling for One-Dimensional Integrals
Both singularity handlers (
IMT and
"DoubleExponential") are applied to regions that are obtained through too many bisections (as specified by
"SingularityDepth").
The main difference between them is that
IMT does not change the integration rule used to compute integral estimates on the region it is applied to—
IMT is only a variable transformation. On the other hand,
"DoubleExponential" uses both variable transformation and a different integration rule&mdas
