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 normally 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 at a set of points, called sampling points (or abscissas).

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 effort (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 with the integrand and regions for which with the integrand .
In[31]:=
Click for copyable input
Out[31]=
Here is a plot of all sampling points used in the integration. The integrand is sampled at the coordinates in the order of the coordinates (in the plot). It can be seen that the sampling points are concentrated near the singularity point . 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.
In[72]:=
Click for copyable input
Out[73]=

"Adaptive Strategies" gives a general description of the adaptive strategies. The default (main) strategy of NIntegrate is , which is explained in "Global Adaptive Strategy". Complementary to it is the local adaptive strategy, which is explained in "Local Adaptive Strategy". Both adaptive strategies use singularity handling mechanisms, which are explained in "Singularity Handling".

The Monte Carlo strategies are explained in "Crude Monte Carlo and Quasi Monte Carlo Strategies" and "Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies".

"GlobalAdaptive"any integrand, adaptive sampling, rule-based
"LocalAdaptive"any integrand, adaptive sampling, rule-based
"DoubleExponential"any integrand, uniform sampling
"Trapezoidal"any integrand, uniform sampling
"MultiPeriodic"multidimensional integrand, uniform sampling
"MonteCarlo"any integrand, uniform random sampling
"QuasiMonteCarlo"any integrand, uniform quasi-random sampling
"AdaptiveMonteCarlo"any integrand, adaptive random sampling
"AdaptiveQuasiMonteCarlo"any integrand, adaptive quasi-random sampling
"DoubleExponentialOscillatory"one-dimensional infinite-range oscillatory integrand
"ExtrapolatingOscillatory"one-dimensional infinite-range oscillatory integrand

NIntegrate integration strategies.

NIntegrate uses certain "preprocessor" strategies for special types of integrals (or integrands). These are explained in "Duffy's Coordinates Strategy", "Oscillatory Strategies", and "Cauchy Principal Value Integration". Preprocessor strategies also handle symbolic preprocessing of the integrand.

"SymbolicPreprocessing"overall preprocessor controller
"EvenOddSubdivision"simplify even and odd integrands
"InterpolationPointsSubdivision"subdivide integrands containing interpolating functions
"OscillatorySelection"detect oscillatory integrands and select suitable methods
"SymbolicPiecewiseSubdivision"subdivide integrands containing piecewise functions
"UnitCubeRescaling"rescale multidimensional integrand to unit cube
"DuffyCoordinates"multidimensional singularity-removing transformation
"PrincipalValue"numerical integral equivalent to Cauchy principal value

NIntegrate preprocessor strategies.

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 and integrand , if an integration region is partitioned into, say, two disjoint subregions and , , , then the sum of the integral estimates of over and is closer to the actual integral . In other words,

and (1) will imply that the sum of the error estimates for and is smaller than the error estimate of .

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 .
In[32]:=
Click for copyable input
Out[32]=
option name
default value
MethodAutomaticintegration rule used to compute integral and error estimates over each subregion
"SingularityDepth"Automaticnumber of recursive bisections before applying a singularity handler
"SingularityHandler"Automaticsingularity handler
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

options.

is the default integration strategy of NIntegrate. It is used for both one-dimensional and multidimensional integration. works with both Cartesian product rules and fully symmetric multidimensional rules.

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".

stops if the following expression is true:

where and 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 and MaxRecursion.

If for any subregion the number of bisections in any of the dimensions is greater than MaxRecursion then the integration by stops.

Setting 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 .

"MaxErrorIncreases"

Since (2) is expected to hold in , 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.

Zero integral.
In[3]:=
Click for copyable input
Out[3]=

(ii) The specified working precision is not dense enough for the specified precision goal.

The working precision is not dense enough.
In[33]:=
Click for copyable input
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 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 option , the integration stops with a message (NIntegrate::eincr).

The default value of is 400 for one-dimensional integrals and 2000 for multidimensional integrals.

The following integration invokes the message NIntegrate::eincr, with the default value of .
In[1]:=
Click for copyable input
Out[1]=
Increasing silences the NIntegrate::eincr message.
In[2]:=
Click for copyable input
Out[2]=
The result compares well with the exact value.
In[3]:=
Click for copyable input
Out[4]=

Example Implementation of a Global Adaptive Strategy

This computes Gauss-Kronrod abscissas, weights, and error weights.
In[15]:=
Click for copyable input
This is a definition of a function that applies the integration rule with abscissas and weights computed to the function f over the interval .
In[16]:=
Click for copyable input
This is a definition of a simple global adaptive algorithm that finds the integral of the function f over the interval with relative error tol.
In[17]:=
Click for copyable input
This defines an integrand.
In[18]:=
Click for copyable input
The global adaptive strategy defined earlier gives the following result.
In[19]:=
Click for copyable input
Out[19]=
Here is the exact result.
In[20]:=
Click for copyable input
Out[20]=
The relative error is within the prescribed tolerance.
In[21]:=
Click for copyable input
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 .
In[5]:=
Click for copyable input
Out[5]=
option name
default value
MethodAutomaticintegration rule used to compute integral and error estimates over the subregions
"SingularityDepth"Automaticnumber of recursive bisections before applying a singularity handler
"SingularityHandler"Automaticsingularity handler
"Partitioning"Automatichow to partition the regions in order to improve their integral estimate
"InitialEstimateRelaxation"Trueattempt to adjust the magnitude of the initial integral estimate in order to avoid unnecessary computation
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

options.

Like , can be used for both one-dimensional and multidimensional integration. works with both Cartesian product rules and fully symmetric multidimensional rules.

The 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, , 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 is the smallest number such that at the working precision, and and are the user-specified precision and accuracy goals.

The recursive routine of 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 criterion (4) is true.

MinRecursion and MaxRecursion

The options and MaxRecursion for have the same meaning and functionality as they do for . See "MinRecursion and MaxRecursion".

"InitialEstimateRelaxation"

After the first recursion is finished a better integral estimate, , will be available. That better estimate is compared to the two integral estimates, and , that the integration rule has used to give the integral estimate () and the error estimate () for the initial step. If

then the integral estimate in (5) can be increased—that is, the condition (6) is relaxed—with the formula

since means that the rule's integral estimate is more accurate than what the rule's error estimate predicts.

"Partitioning"

has the option to specify how to partition the regions that do not satisfy (7). For one-dimensional integrals, if is set to Automatic, 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 is given a list of integers with length that equals the number of integral variables, each dimension of the integration region is divided into equal parts. If is given an integer , all dimensions are divided into equal parts.

Consider the following function.
In[4]:=
Click for copyable input
Out[4]=
These are the sampling points used by 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.
In[1]:=
Click for copyable input
Out[2]=
These are the sampling points used by 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.
In[5]:=
Click for copyable input
Out[6]=
Multidimensional example of using the option. To make the plot, the sampling points of the first region to be integrated, , are removed.
In[7]:=
Click for copyable input
Out[10]=

Reuse of Integrand Values

With its default partitioning settings for one-dimensional integrals reuses the integrand values at the endpoints of the subintervals that have integral and error estimates that do not satisfy (8).

Sampling points of the integration of by . The variable determines the number of points in the integration rule used by .
In[13]:=
Click for copyable input
Out[15]=
Out[16]=
The percent of reused points in the integration.
In[17]:=
Click for copyable input
Out[19]=

Example Implementation of a Local Adaptive Strategy

This computes Clenshaw-Curtis abscissas, weights, and error weights.
In[33]:=
Click for copyable input
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 .
In[34]:=
Click for copyable input
This defines a simple local adaptive algorithm that finds the integral of the function f over the interval with relative error tol.
In[35]:=
Click for copyable input
This defines a function.
In[37]:=
Click for copyable input
The local adaptive strategy gives the result.
In[38]:=
Click for copyable input
Out[38]=
This is the exact result.
In[39]:=
Click for copyable input
Out[39]=
The relative error is within the prescribed tolerance.
In[40]:=
Click for copyable input
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 and :

1. stops when the sum of the errors of all regions satisfies the precision goal, while stops when the error of each region is small enough compared to an estimate of the integral.

2. To improve the integral estimate bisects the region with largest error, while partitions all regions for which the error is not small enough.

For multidimensional integrals is much faster because 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 is faster and performs better than , it is because the precision-goal-stopping criteria and partitioning strategy of are more suited for the integrand's nature. Another factor is the ability of to reuse the integral values of all points already sampled. 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 and .

"GlobalAdaptive" Is Generally Better than "LocalAdaptive"

The table that follows, with timing ratios and numbers of integrand evaluations, demonstrates that is better than for the most common cases. All integrals considered are one dimensional over , because: (1) for multidimensional integrals the performance differences should be expected to deepen, since partitions the regions along each axis; and (2) one-dimensional integrals over different ranges can be rescaled to be over .

Here are the definitions of some functions, precision goals, number of integrations, and the integration rule. The variable can be changed in order to compare the profiling runs with the same integration rule. The last function is derived from by the variable change that transforms into .
In[70]:=
Click for copyable input
In[74]:=
Click for copyable input
Out[74]=
Exact integral values. All integrals are over .
In[75]:=
Click for copyable input
timings.
In[76]:=
Click for copyable input
timings.
In[77]:=
Click for copyable input
function evaluations.
In[78]:=
Click for copyable input
function evaluations.
In[79]:=
Click for copyable input
Table with the timing ratios and integrand evaluations.
In[80]:=
Click for copyable input
Out[80]=
Table with the errors of the integrations. Both and reach the required precision goals.
In[81]:=
Click for copyable input
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 .
In[58]:=
Click for copyable input
Out[58]=
Here is an example of a two-dimensional integral with a singular point at .
In[59]:=
Click for copyable input
Out[59]=
Here is an example of an integral that has two singular points at and specified with the Exclusions option.
In[60]:=
Click for copyable input
Out[60]=
Here is an example of a two-dimensional integral with a singular point at specified with the Exclusions option.
In[61]:=
Click for copyable input
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 .
In[62]:=
Click for copyable input
Out[62]=
Using Exclusions the integral is quickly calculated.
In[12]:=
Click for copyable input
Out[12]=
NIntegrate will reach convergence much more slowly if no singularity specification is given.
In[35]:=
Click for copyable input
Out[35]=
Here is an example of a case in which a singular curve can be specified with the variable ranges. If and this would not be possible—see the following example.
In[10]:=
Click for copyable input
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 .
In[66]:=
Click for copyable input
Out[66]=
Using Boole the integral is calculated quickly.
In[9]:=
Click for copyable input
Out[9]=
This two-dimensional function has singular points along the curve .
In[68]:=
Click for copyable input
Out[68]=
Again, using Boole the integral is calculated quickly.
In[8]:=
Click for copyable input
Out[8]=
This is how the sampling points of the integration look.
In[6]:=
Click for copyable input
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.
In[1]:=
Click for copyable input
This defines a three-dimensional function.
In[2]:=
Click for copyable input
Here is the integral of a three-dimensional function with singular points along the surface .
In[3]:=
Click for copyable input
Out[3]=
These are the sampling points of the integration.
In[4]:=
Click for copyable input
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 , it is decided that subregion has a singularity. Then the integration over that subregion is done with the singularity handler specified by .

option name
default value
"SingularityDepth"Automaticnumber of recursive bisections before applying a singularity handler
"SingularityHandler"Automaticsingularity handler

and 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", ) 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 endpoints 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.
In[13]:=
Click for copyable input
Out[13]=
option name
default value
"TuningParameters"10a pair of numbers that are the tuning parameters in the IMT transformation formula ; if only a number a is given, it is interpreted as

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 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 , is defined.
In[14]:=
Click for copyable input

The parameters a and p are called tuning parameters [MurIri82].

The limit of the derivative of the IMT transformation is 0.
In[16]:=
Click for copyable input
Out[16]=
Here is the plot of the IMT transformation.
In[17]:=
Click for copyable input
Out[17]=

From the graph above, it 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.
In[18]:=
Click for copyable input
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 that has a singularity at 0.
In[29]:=
Click for copyable input
In[30]:=
Click for copyable input
Out[30]=
Assume the integration is done with , with singularity handler IMT and singularity depth 4. After four bisections will have a region with boundaries that contains the singular endpoint. For that region the IMT variable transformation will change its boundaries to and its integrand to the following.
In[31]:=
Click for copyable input
Out[32]=
Here is the plot of the new integrand.
In[33]:=
Click for copyable input
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".

For example, consider the sampling points and weights of the Gauss-Kronrod rule.
In[35]:=
Click for copyable input
The Gauss-Kronrod sampling points for the region and the derivatives of the rescaling follow.
In[36]:=
Click for copyable input
Out[36]=
In[37]:=
Click for copyable input
Out[37]=
Here is the integral estimate.
In[38]:=
Click for copyable input
Out[38]=
With the IMT transformation, these are the sampling points and derivatives.
In[39]:=
Click for copyable input
Out[39]=
In[40]:=
Click for copyable input
Out[40]=
Here is the integral estimate with the IMT transformation.
In[41]:=
Click for copyable input
Out[41]=
The estimate calculated with the IMT variable transformation is much closer to the exact value.
In[42]:=
Click for copyable input
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.
In[103]:=
Click for copyable input
Out[103]=

IMT versus "DoubleExponential" versus No Singularity Handling for One-Dimensional Integrals

Both singularity handlers ( and ) are applied to regions that are obtained through too many bisections (as specified by ).

The main difference between them is that does not change the integration rule used to compute integral estimates on the region it is applied to— is only a variable transformation. On the other hand, uses both variable transformation and a different integration rule—the trapezoidal rule—to compute integral estimates on the region it is applied to. In other words, the singularity handler delegates the integration to the double exponential quadrature as described in "Double Exponential Strategy".

As a consequence, a region to which the singularity handler is applied is still going to be subject to bisection by the adaptive integration strategy. Therefore, until the precision goal is reached the integrand evaluations done before the last bisection will be thrown away. On the other hand, a region to which the singularity handler is applied will not be bisected. The trapezoidal rule quadrature used by will compute integral estimates over the region with an increasing number of sampling points at each step, completely reusing the integrand evaluations of the sampling points from the previous steps.

So, if the integrand is "very" analytic (i.e. no rapid or sudden changes of the integrand and its derivatives with respect to the integration variable) over the regions with endpoint singularity, the singularity handler is going to be much faster than the singularity handler. In the cases where the integrand is not analytic in the region given to the singularity handler, or the double transformation of the integrand converges too slowly, it is better to switch to the singularity handler. This is done if the option is set to Automatic.

Following are tables that compare the , , and Automatic singularity handlers applied at different depths of bisection.

This loads a package that defines the profiling function that gives the number of sampling points and the time needed by a numerical integration command.
In[17]:=
Click for copyable input
Table for a "very" analytical integrand that the singularity handler easily computes.
In[34]:=
Click for copyable input
Out[36]//TableForm=
Table for an integrand that does not have a singularity and has a nearly discontinuous derivative (i.e. it is not "very" analytical). The Automatic singularity handler starts with and then switches to .
In[37]:=
Click for copyable input
Out[40]//TableForm=
A table for an integrand for which the Automatic singularity handler starts with and then switches to .
In[41]:=
Click for copyable input
Out[44]//TableForm=

IMT Multidimensional Singularity Handling

When used for multidimensional integrals, the IMT singularity handler speeds up the integration process only when the singularity is along one of the axes. When the singularity is at a corner of the integration region, using IMT is counterproductive. The function defined earlier is used in the following examples.

The number of integrand evaluations and timings for an integrand that has a singularity only along the axis. The default (automatic) singularity handler chooses to apply IMT to regions obtained after the default (four) bisections.
In[19]:=
Click for copyable input
Out[19]=
The number of integrand evaluations and timings for an integrand that has a singularity only along the axis with no singularity handler application.
In[20]:=
Click for copyable input
Out[20]=
The number of integrand evaluations and timings for an integrand that has a singularity at a corner of the integration region. The default (automatic) singularity handler chooses to apply the singularity handler DuffyCoordinates to regions obtained after the default (four) bisections.
In[21]:=
Click for copyable input
Out[21]=
The number of integrand evaluations and timings for an integrand that has a singularity at a corner of the integration region. IMT is applied to regions obtained after the default (four) bisections.
In[22]:=
Click for copyable input
Out[22]=
The number of integrand evaluations and timings for an integrand that has a singularity at a corner of the integration region with no singularity handler application.
In[23]:=
Click for copyable input
Out[23]=

Duffy's Coordinates for Multidimensional Singularity Handling

Duffy's coordinates is a technique that transforms an integrand over a square, cube, or hypercube with a singular point in one of the corners into an integrand with a singularity over a line, which might be easier to integrate.

The following integration uses Duffy's coordinates.
In[63]:=
Click for copyable input
Out[63]=
The following integration does not use Duffy's coordinates.
In[62]:=
Click for copyable input
Out[62]=

The NIntegrate strategies and apply the Duffy's coordinates technique only at the corners of the integration region.

When the singularity of a multidimensional integral occurs at a point, the coupling of the variables will make the singularity variable transformations used in one-dimensional integration counterproductive. A variable transformation that has a geometrical nature, proposed by Duffy in [Duffy82], makes a change of variables that replaces a point singularity at a corner of the integration region with a "softer" one on a plane.

If is the dimension of integration and , then Duffy's coordinates is a suitable technique for singularities of the following type (see again [Duffy82]):

1.  , , ;

2.  , , ;

3.  , , , .

For example, consider the integral

If the integration region is changed to with the rule , the Jacobian of which is , the integral becomes

The last integral has no singularities at all!

Now consider the integral

which is equivalent to the sum

The first integral of that sum is transformed as in (9); for the second one, though, the change of into by has the Jacobian , which will not bring the desired cancellation of terms. Fortunately, a change of the order of integration:

makes the second integral amenable for the transformation in (10):

(In the second integral in equation (3) the variables were permuted, which is not necessary to prove the mathematical equivalence, but it is faster when computing the integrals.)

So the integral (11) can be rewritten as an integral with no singularities:

If the integration variables were not permuted in (12), the integral (13) is going to be rewritten as

That is a more complicated integral, as its integrand is not simple along both axes. Subsequently it is harder to compute than the former one.

Here is the number of sampling points for the simpler integral.
In[58]:=
Click for copyable input
Out[58]=
Here is the number of sampling points for the more complicated integral.
In[59]:=
Click for copyable input
Out[59]=

NIntegrate uses a generalization to arbitrary dimension of the technique in the example above. (In [Duffy82] only the third dimension is described.) An example implementation together with the generalization description are given below.

Here is a table that compares the different singularity handling methods for . (The profiling function defined earlier is used.)
In[74]:=
Click for copyable input
Out[75]=

Duffy's Coordinates Strategy

When Duffy's coordinates are applicable, a numerical integration result is obtained faster if Duffy's coordinate change is made before the actual integration begins. Making the transformation beforehand, though, requires knowing at which corner(s) of the integration region the singularities occur. The strategy in NIntegrate facilitates such pre-integration transformation.

Here is an example with an integrand that has singularities at two different corners of its integration region.
In[84]:=
Click for copyable input
Out[84]=
option name
default value
Method{"GlobalAdaptive","SingularityDepth"->∞}the strategy with which the integration will be made after applying Duffy's coordinates transformation
"Corners"Alla vector or a list of vectors that specifies the corner(s) to apply the Duffy's coordinates tranformation to; the elements of the vectors are either 0 or 1; each vector length equals the dimension of the integral

options.

The first thing does is to rescale the integral into one that is over the unit hypercube (or square, or cube). If only one corner is specified, applies Duffy's coordinates transformation as described earlier. If more than one corner is specified, the unit hypercube of the previous step is partitioned into disjoint cubes with side length of one-half. Consider the integrals over these disjoint cubes. Duffy's coordinates transformation is applied to the ones that have a vertex that is specified to be singular. The rest are transformed into integrals over the unit cube. Since all integrals at this point have an integration region that is the unit cube, they are summated, and that sum is given to NIntegrate with a Method option that is the same as the one given to .

The actual integrand used by can be obtained through , which has the same arguments as NIntegrate.

Here is an example for the integrand of a three-dimensional function that is singular at one of the corners of the integration region.
In[78]:=
Click for copyable input
Out[78]=
Here is an example for the integrand for a two-dimensional function that is singular at two of the corners of the integration region.
In[79]:=
Click for copyable input
Out[79]=

might considerably improve speed for the types of integrands described in "Duffy's Coordinates for Multidimensional Singularity Handling".

Integration with .
In[80]:=
Click for copyable input
Out[80]=
Integration with the default NIntegrate options settings, which is much slower than the previous one.
In[81]:=
Click for copyable input
Out[81]=
Here is another example of a speedup by .
In[82]:=
Click for copyable input
Out[82]=
Integration with the default NIntegrate options settings, which is much slower than the previous one.
In[83]:=
Click for copyable input
Out[83]=

Duffy's Coordinates Generalization and Example Implementation

See "Duffy's Coordinates for Multidimensional Singularity Handling" for the theory of Duffy's coordinates.

The implementation is based on the following two theorems.

Theorem 1: A -dimensional cube can be divided into disjoint geometrically equivalent -dimensional pyramids (with bases -dimensional cubes) and with coinciding apexes.

Proof: Assume the side length of the cube is 1, the cube has a vertex at the origin, and the coordinates of all other vertexes are or . Consider the -dimensional cube walls , where . Their number is exactly , and the origin does not belong to them. Each of the walls can form a pyramid with the origin. This proves the theorem.

Here is a plot that illustrates the theorem in 3D.
In[89]:=
Click for copyable input
Out[92]=

If the axes are denoted the pyramid formed with the wall can be described as . Let denote the permutation derived after rotating cyclically times to the left (i.e. applying times RotateLeft to ). Then the following theorem holds:

Theorem 2: For any integral over the unit cube the following equalities hold:

Proof: The first equality follows from Theorem 1. The second equality is just a change of variables that transforms a pyramid to a cube.

Here is a function that gives the rules and the Jacobian for the transformation of a hypercube with a specified side into a region.
In[93]:=
Click for copyable input
Here is an example of unit-square to infinite-region rescaling.
In[96]:=
Click for copyable input
Out[96]=
Here is a function that computes the integrals obtained by the Duffy's coordinates technique when the singularity is at the origin.
In[97]:=
Click for copyable input
Here is a function that computes the integrals obtained by the Duffy's coordinates technique for a specified corner of the hypercube where the singularity occurs.
In[99]:=
Click for copyable input
Here is a symbolic example.
In[101]:=
Click for copyable input
Out[101]=
Here is another symbolic example.
In[102]:=
Click for copyable input
Out[102]=
Here is a computational example.
In[103]:=
Click for copyable input
Out[103]=
Using Duffy's coordinates is much faster than using no singularity handling (see the next example).
In[108]:=
Click for copyable input
Out[110]=
Integration using no singularity handling.
In[111]:=
Click for copyable input
Out[111]=
Of course, the internal implementation of NIntegrate gives similar performance results.
In[107]:=
Click for copyable input
Out[107]=

Ignoring the Singularity

Another way of handling a singularity is to ignore it. NIntegrate ignores sampling points that coincide with a singular point.

Consider the following integral that has a singular point at 1.

The integrand goes to when the integration variable is close to 1.

Here is a plot of the integrand.
In[118]:=
Click for copyable input
Out[118]=
NIntegrate gives a result that is close to the exact one.
In[114]:=
Click for copyable input
Out[115]=
Convergence is achieved when MaxRecursion is increased.
In[45]:=
Click for copyable input
Out[45]=

With its default options NIntegrate has a sampling point at 1, as can be seen from the following.

Check that NIntegrate has 1 as a sampling point.
In[119]:=
Click for copyable input
Out[119]=

But for NIntegrate[Log[(1-x)2], {x, 0, 2}] the evaluation monitor has not picked a sampling point that is 1.

Sampling points that belong to the interval .
In[120]:=
Click for copyable input
Out[120]=

In other words, the singularity at 1 is ignored. Ignoring the singularity is equivalent to having an integrand that is zero at the singular sampling point.

Note that the integral is easily integrated if the singular point is specified in the variable range. Following are the numbers of sampling points and timings for NIntegrate with the singular and nonsingular range specifications.

Integration with the singular point specified.
In[123]:=
Click for copyable input
Out[123]=
Integration by ignoring the singularity.
In[122]:=
Click for copyable input
Out[122]=

A more interesting example of ignoring the singularity is using Bessel functions in the denominator of the integrand.

Integral with several (five) integrable singularities.
In[124]:=
Click for copyable input
Out[124]//InputForm=

The result can be checked using NIntegrate with singular range specification with the zeros of BesselJ[2, x] (see BesselJZero).

Integration with the Bessel zeros specified as singular points.
In[125]:=
Click for copyable input
Out[125]//InputForm=

Needless to say, the last integration required the calculation of the BesselJ zeros. The former one "just integrates" without any integrand analysis.

Ignoring the singularity may not work with oscillating integrands.

For example, these two integrals are equivalent.
In[126]:=
Click for copyable input
Out[126]=
NIntegrate can do the first one.
In[127]:=
Click for copyable input
Out[127]=
NIntegrate cannot do the second one.
In[128]:=
Click for copyable input
Out[128]=

However, if the integrand is monotonic in a neighborhood of its singularity, or more precisely, if it can be majorized by a monotonic integrable function, it can be shown that by ignoring the singularity, convergence will be reached.

For theoretical justification and practical recommendations of ignoring the singularity see [DavRab65IS] and [DavRab84].

Automatic Singularity Handling

One-Dimensional Integration

When the option is set to Automatic for a one-dimensional integral, is used on regions that are obtained by number of partitionings. As explained earlier, this region will not be partitioned further as long as the singularity handler works over it. If the error estimate computed by does not evolve in a way predicted by the theory of the double exponential quadrature, the singularity handling for this region is switched to .

As explained in "Convergence Rate", the following dependency of the error is expected with respect to the number of double exponential sampling points:

where is a positive constant. Consider the relative errors and of two consecutive double exponential quadrature calculations, made with and number of sampling points respectively, for which . Assuming , , and it should be expected that

The switch from to happens when:

(i) the region error estimate is larger than the absolute value of the region integral estimate (hence the relative error is not smaller than 1);

(ii) the inequality (2) is not true in two different instances;

(iii) the integrand values calculated with the double exponential transformation do not decay fast enough.

Here is an example of a switch from to singularity handling. On the plot the integrand is sampled at the coordinates in the order of the coordinates. The patterns of the sampling points over show the change from Gaussian quadrature () to double exponential quadrature (), which later is replaced by Gaussian quadrature using the IMT variable transformation ().
In[143]:=
Click for copyable input
Out[145]=

Multidimensional Integration

When the option is set to Automatic for a multidimensional integral, both and are used.

A region needs to meet the following conditions in order for to be applied:

the region is obtained by number of bisections (or partitionings) along each axis;

the region is a corner of one of the initial integration regions (the specified integration region can be partitioned into integration regions by piecewise handling or by user-specified singularities).

A region needs to meet the following conditions in order for to be applied:

the region is obtained with number of bisections (or partitionings) along predominantly one axis;

the region is not a corner region and it is on a side of one of the initial integration regions.

In other words, is applied to regions that are derived through number of partitionings but that do not satisfy the conditions of the automatic application.

is effective if the singularity is along one of the axes. Using for point singularities can be counterproductive.

Sampling points of two-dimensional integration, , with Automatic (left) and (right) singularity handling. It can be seen that the automatic singularity handling uses almost two times more points than . To illustrate the effect of the singularity handlers they are applied after two bisections.
In[133]:=
Click for copyable input
Out[134]=
Timings for the integral with singularity handlers Automatic, , and , and with no singularity handling. The integral has a point singularity at .
In[47]:=
Click for copyable input
Out[47]//TableForm=
Timings for the integral , singular along the axis with singularity handlers Automatic, , and , and with no singularity handling.
In[46]:=
Click for copyable input
Out[46]//TableForm=

Cauchy Principal Value Integration

To evaluate the Cauchy principal value of an integral, NIntegrate uses the strategy PrincipalValue.

Cauchy principal value integration with singular point at 2.
In[153]:=
Click for copyable input
Out[153]=

In NIntegrate, PrincipalValue uses the strategy specified by its Method option to work directly on those regions where there is no difficulty and by pairing values symmetrically about the specified singularities in order to take advantage of the cancellation of the positive and negative values.

option name
default value
MethodAutomaticmethod specification used to compute estimates over subregions
SingularPointIntegrationRadiusAutomatica number or a list of numbers that correspond to the singular points , , ..., in the range specification; with each pair an integral of the form is formed

options.

Thus the specification

is evaluated as

where each of the integrals is evaluated using NIntegrate with Method->methodspec. If is not given explicitly, a value is chosen based upon the differences and . The option can take a list of numbers that equals the number of singular points. For the derivation of the formula see [DavRab84].

This finds the Cauchy principal value of .
In[14]:=
Click for copyable input
Out[14]=
Here is the Cauchy principal value of . Note that there are two singularities that need to be specified.
In[114]:=
Click for copyable input
Out[114]=
The singular points can be specified using the Exclusions option.
In[30]:=
Click for copyable input
Out[30]=
This checks the value. The result would be 0 if everything were done exactly.
In[31]:=
Click for copyable input
Out[31]=

It should be noted that the singularities must be located exactly. Since the algorithm pairs together the points on both sides of the singularity, if the singularity is slightly mislocated the cancellation will not be sufficiently good near the pole and the result can be significantly in error if NIntegrate converges at all.

Sampling Points Visualization

Consider the calculation of the principal value of

The following examples show two ways of visualizing the sampling points. The first shows the sampling points used. Since the integrand is modified in order to do the principal value integration, it might be desired to see the points at which the original integrand is evaluated. This is shown on the second example.

Here are sampling points used by NIntegrate. There are no points over the interval because of the PrincipalValue integration , and there are sampling points over .
In[154]:=
Click for copyable input
Out[156]=
This defines a function that accumulates the argument values given to the integrand.
In[1]:=
Click for copyable input
Here are the points at which the integrand has been evaluated. Note the symmetric pattern over the interval .
In[166]:=
Click for copyable input
Out[168]=

Double Exponential Strategy

The double exponential quadrature consists of applying the trapezoidal rule after a variable transformation. The double exponential quadrature was proposed by Mori and Takahasi in 1974 and it was inspired by the so-called IMT rule and TANH rule. The transformation is given the name "double exponential" since its derivative decreases double exponentially when the integrand's variable reaches the ends of the integration region.

The double exponential algorithm for NIntegrate is specified with the Method option value .
In[169]:=
Click for copyable input
Out[169]=
option name
default value
"ExtraPrecision"50maximum extra precision to be used internally
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

options.

The double exponential strategy can be used for one-dimensional and multidimensional integration. When applied to multidimensional integrals it uses the Cartesian product of the trapezoidal rule.

A double exponential transformation transforms the integral

into

where can be finite, half-infinite (), or infinite (). The integrand must be analytic in and might have singularity at one or both of the endpoints.

The transformed integrand decreases double exponentially, that is, as .

The function is analytic in . It is known that for an integral like (14) of an analytic integrand the trapezoidal rule is an optimal rule [Mori74].

The transformations used for the different types of integration regions are:

where and are finite numbers.

The trapezoidal rule is applied to (15):

The terms in (16) decay double exponentially for large enough . Therefore the summation in (17) is cut off at the terms that are too small to contribute to the total sum. (A criterion similar to (18) for the local adaptive strategy is used. See also the following double exponential example implementation.)

The strategy employs the double exponential quadrature.

The strategy works best for analytic integrands; see "Comparison of Double Exponential and Gaussian Quadrature".

uses the Cartesian product of double exponential quadratures for multidimensional integrals.

Cartesian double exponential quadrature.
In[48]:=
Click for copyable input
Out[48]=

As with the other Cartesian product rules, if is used for dimensions higher than three, it might be very slow due to combinatorial explosion.

The following plot illustrates the Cartesian product character of the multidimensional integration.
In[49]:=
Click for copyable input
Out[50]=

Double exponential quadrature can be used for singularity handling in adaptive strategies; see "Singularity Handling".

MinRecursion and MaxRecursion

The option has the same meaning and functionality as it does for and described in "MinRecursion and MaxRecursion". MaxRecursion for restricts how many times the trapezoidal quadrature estimates are improved; see "Example Implementation of Double Exponential Quadrature".

Comparison of Double Exponential and Gaussian Quadrature

The strategy works best for analytic integrands. For example, the following integral is done by three times faster than the Gaussian quadrature (using a global adaptive algorithm).

Integration with .
In[215]:=
Click for copyable input
Out[215]=
Integration with Gauss quadrature. (The default strategy of NIntegrate, , uses by default a Gauss-Kronrod integration rule with 5 Gaussian points and 6 Kronrod points.)
In[51]:=
Click for copyable input
Out[51]=

Since converges double exponentially with respect to the number of evaluation points, increasing the precision goal slightly increases the work done by . This is illustrated for two integrals, and . Each table entry shows the error and number of evaluations.

Double exponential quadrature and Gaussian quadrature for . Increasing the precision goal does not change the number of sampling points used by .
In[217]:=
Click for copyable input
Out[219]=
Double exponential quadrature and Gaussian quadrature for . Increasing the precision goal does not change the number of sampling points used by . (The integrations are done without symbolic preprocessing.)
In[220]:=
Click for copyable input
Out[222]=

On the other hand, for nonanalytic integrands is quite slow, and a global adaptive algorithm using Gaussian quadrature can resolve the singularities easily.

needs more than 10000 integrand evaluations to compute this integral with a nonanalytic integrand.
In[52]:=
Click for copyable input
Out[53]=
Gaussian quadrature is much faster for the integral.
In[54]:=
Click for copyable input
Out[54]=

Further, might be slowed down by integrands that have nearly discontinuous derivatives, that is, integrands that are not "very" analytical.

Here is an example with a not "very" analytical integrand.
In[226]:=
Click for copyable input
Out[226]=
Again, Gaussian quadrature is much faster.
In[227]:=
Click for copyable input
Out[227]=
Here are the plots of the integrand and its derivative.
In[228]:=
Click for copyable input
Out[228]=

Convergence Rate

This section demonstrates that the asymptotic error of the double exponential quadrature in terms of the number of evaluation points used is

where is a positive constant.

This defines a double exponential integration function that returns an integral estimate and the number of points used.
In[229]:=
Click for copyable input
This defines a function.
In[230]:=
Click for copyable input
This is the exact integral.
In[231]:=
Click for copyable input
Out[231]=
This finds the errors and number of evaluation points for a range of step sizes of the trapezoidal rule.
In[232]:=
Click for copyable input
This fits through the logarithms of the errors; see (19).
In[239]:=
Click for copyable input
Here is the fitted function. The summation term is just a translation parameter.
In[240]:=
Click for copyable input
Out[240]=
You see that the errors fit the model (20).
In[241]:=
Click for copyable input
Out[241]=

Example Implementation of Double Exponential Quadrature

Following is an example implementation of the double exponential quadrature with the finite-region variable transformation (transformation (21) earlier).

This is a definition of a function that applies the trapezoidal rule to a transformed integrand. The function uses (22) and it is made to reuse integral estimates computed with a twice-larger step.
In[173]:=
Click for copyable input
This is a definition of a simple double exponential strategy, which finds the integral of the function f over the finite interval with relative error tol.
In[189]:=
Click for copyable input
This defines a function that is singular at 0.
In[194]:=
Click for copyable input
Here is the integral estimate from the double exponential strategy.
In[195]:=
Click for copyable input
Out[195]//InputForm=
Here is the exact result.
In[177]:=
Click for copyable input
Out[177]//InputForm=

The two results are the same.

This defines an oscillating function.
In[178]:=
Click for copyable input
Here is the integral estimate given by the double exponential strategy.
In[179]:=
Click for copyable input
Out[179]//InputForm=
Here is the exact result.
In[180]:=
Click for copyable input
Out[180]=
Here is the exact result in machine precision.
In[181]:=
Click for copyable input
Out[181]//InputForm=
The relative error is within the prescribed tolerance.
In[182]:=
Click for copyable input
Out[182]=

"Trapezoidal" Strategy

The strategy gives optimal convergence for analytic periodic integrands when the integration interval is exactly one period.

option name
default value
"ExtraPrecision"50maximum extra precision to be used internally
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

options.

takes the same options as . If the integration ranges are infinite or semi-infinite, becomes .

For theoretical background, examples, and explanations of periodic functions integration (with trapezoidal quadrature) see [Weideman2002].

In[109]:=
Click for copyable input
Out[109]=
Here is a table that shows the number of sampling points for different values of the parameter used by and respectively for the integral .
In[33]:=
Click for copyable input
Out[35]//TableForm=

Example Implementation

This function makes a trapezoidal quadrature integral estimate with specified points.
In[242]:=
Click for copyable input
This function improves a trapezoidal quadrature integral estimate using sampling points between the old ones.
In[257]:=
Click for copyable input
This function is an interface to the preceding one.
In[272]:=
Click for copyable input
Here is a definition of a (Bessel) function.
In[269]:=
Click for copyable input
Here is the trapezoidal quadrature estimate.
In[274]:=
Click for copyable input
Out[274]=
Here is the exact value.
In[278]:=
Click for copyable input
Out[278]=
The relative error is within the prescribed tolerance.
In[279]:=
Click for copyable input
Out[279]=

Oscillatory Strategies

NIntegrate includes several methods for oscillatory integration. The integration rule "LevinRule" applies to a large class of one-dimensional and multidimensional oscillatory integrals. It is used with adaptive rule-based strategies such as .

NIntegrate also includes several integration strategies specialized to certain types of one-dimensional oscillatory integrals. Generally these algorithms are for either finite-region integrals or infinite-region integrals only. For oscillatory integrals on a finite region, NIntegrate uses Chebyshev expansions of the integrand and the global adaptive integration strategy. For oscillatory integrals on an infinite region, NIntegrate uses either a modification of the double exponential algorithm or sequence summation acceleration over the sequence of integrals over regions between the zeros of the integrand.

The following example uses all the different specialized oscillatory strategies on different subregions.
In[77]:=
Click for copyable input
Out[77]=

NIntegrate automatically detects oscillatory (one-dimensional) integrands, and selects which strategy or rule to use according to the integrand's region and specific oscillatory form.

The integrals handled by the specialized strategies described here are of the form

where the oscillating kernel is of the form:

1.  , , for finite;

2.  , , , , , , , , or for infinite or semi-infinite.

In these oscillating kernel forms, , , and are real constants, and is a positive integer.

See "LevinRule" for a description of the more general forms of oscillatory integrand handled by that integration rule.

Finite Region Oscillatory Integration

Modified Clenshaw-Curtis quadrature ([PiesBrand75][PiesBrand84]) is for finite-region one-dimensional integrals of the form

where , , , , are finite real numbers.

The modified Clenshaw-Curtis quadrature rule approximates with a single polynomial through Chebyshev polynomials expansion. This leads to simplified computations because of the orthogonality of the Chebyshev polynomials with sine and cosine functions. The modified Clenshaw-Curtis quadrature rule is used with the strategy . For smooth the modified Clenshaw-Curtis quadrature is usually superior [KrUeb98] to other approaches for oscillatory integration (as Filon's quadrature and multi-panel integration between the zeros of the integrand).

Modified Clenshaw-Curtis quadrature is quite good for highly oscillating integrals of the form (23). For example, modified Clenshaw-Curtis quadrature uses less than a hundred integrand evaluations for both and .

Number of integrand evaluations for modified Clenshaw-Curtis quadrature for slowly oscillating kernel.
In[1]:=
Click for copyable input
Out[1]=
Timing and integral estimates for modified Clenshaw-Curtis quadrature for slowly oscillating kernel.
In[3]:=
Click for copyable input
Out[3]=
Number of integrand evaluations for modified Clenshaw-Curtis quadrature for highly oscillating kernel.
In[5]:=
Click for copyable input
Out[5]=
Timing and integral estimates for modified Clenshaw-Curtis quadrature for highly oscillating kernel.
In[6]:=
Click for copyable input
Out[6]=

On the other hand, without symbolic preprocessing, the default NIntegrate method— strategy with a Gauss-Kronrod rule—uses thousands of evaluations for , and it cannot integrate .

Number of integrand evaluations for Gaussian quadrature for slowly oscillating kernel.
In[7]:=
Click for copyable input
Out[7]=
Timing and integral estimates for Gaussian quadrature for slowly oscillating kernel.
In[8]:=
Click for copyable input
Out[8]=
Number of integrand evaluations for Gaussian quadrature for highly oscillating kernel.
In[9]:=
Click for copyable input
Out[9]=
Timing and integral estimates for Gaussian quadrature for highly oscillating kernel.
In[10]:=
Click for copyable input
Out[10]=

Extrapolating Oscillatory Strategy

The NIntegrate strategy is for oscillating integrals in infinite one-dimensional regions. The strategy uses sequence convergence acceleration for the sum of the sequence that consists of each of the integrals with regions between two consecutive zeros of the integrand.

Here is an example of an integration using .
In[294]:=
Click for copyable input
Out[294]=
option name
default value
MethodGlobalAdaptiveintegration strategy used to integrate between the zeros and which will be used if fails
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

Consider the integral

where the function is the oscillating kernel and the function is smooth. Let be the zeros of enumerated from the lower (finite) integration bound, that is, the inequality holds. If the integral (24) converges then the sequence

converges too. The elements of the sequence (25) are the partial sums of the sequence

Often a good estimate of the limit of the sequence (26) can be computed with relatively few elements of it through some convergence acceleration technique.

The strategy uses NSum with Wynn's extrapolation method for the integrals in (27). Each integral in (28) is calculated by NIntegrate without oscillatory methods.

The strategy applies its algorithm to oscillating kernels in (29) that are of the form , , , , , , , or , where , , , and are real constants.

Example Implementation

The following example implementation illustrates how the strategy works.

Here is a definition of an oscillation function that will be integrated in the interval . The zeros of the oscillating function are .
In[1]:=
Click for copyable input
Here is a plot of the oscillatory function in the interval .
In[89]:=
Click for copyable input
Out[89]=
This is a definition of a function that integrates between two consequent zeros. The zeros of the oscillating function are .
In[5]:=
Click for copyable input
Here is the integral estimate computed by sequence convergence acceleration (extrapolation).
In[6]:=
Click for copyable input
Out[6]=
Here is the exact integral value.
In[7]:=
Click for copyable input
Out[7]=
The integral estimate is very close to the exact value.
In[8]:=
Click for copyable input
Out[8]=
Here is another check using the strategy.
In[94]:=
Click for copyable input
Out[94]=
The integral estimate by is very close to the exact value.
In[95]:=
Click for copyable input
Out[95]=

Double Exponential Oscillatory Integration

The strategy is for slowly decaying oscillatory integrals over one-dimensional infinite regions that have integrands of the form , , or , where is the integration variable and , , are constants.

Integration with .
In[2]:=
Click for copyable input
Out[2]=
option name
default value
MethodNoneintegration strategy used if fails
"TuningParameters"Automatictuning parameters of the error estimation
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

options.

is based on the strategy "DoubleExponential", but instead of using a transformation that reaches double exponentially the ends of the integration interval, uses a transformation that reaches double exponentially the zeros of and . The theoretical foundations and properties of the algorithm are explained in [OouraMori91], [OouraMori99], [MoriOoura93]. The implementation of uses the formulas and the integrator design in [OouraMori99].

The algorithm of will be explained using the sine integral

Consider the following transformation

where and are constants satisfying

The parameters and are chosen to satisfy

(taken from [OouraMori99]).

Transformation (30) is applied to (31) to obtain

Note that disappeared in the sine term. The trapezoidal formula with equal mesh size applied to (32) gives

which is approximated with the truncated series sum

and are chosen to satisfy

The integrand decays double exponentially at large negative , as can be seen from (33). While the double exponential transformation, (34) in "Double Exponential Strategy", also makes the integrand decay double exponentially at large positive , the transformation (35) does not decay the integrand at large positive . Instead it makes the sampling points approach double exponentially to the zeros of at large positive . Moreover

As is explained in [OouraMori99], since is linear near any of its zeros, the integrand decreases double exponentially as approaches a zero of . This is the sense in which (36) is considered a double exponential formula.

The relative error is assumed to satisfy

In [OouraMori99] the suggested value for is 5.

Since the formulas cannot be made progressive, (as proposed in [OouraMori99]) does between two and four integration estimates with different . If the desired relative error is the integration steps are the following:

1. Choose such that

and compute (37) with . Let the result be .

2. Next, set , and compute (38) with . Let the result be . The relative error of the first integration step is assumed to be . From (39) , and therefore, if

is satisfied, where is a robustness factor (by default 10) exits with result .

3. If (40) does not hold, compute

and compute (41) with . If

exits with result .

4. If (42) does not hold, compute

and compute (43) with . Let the result be . If

does not hold, issues the message NIntegrate::deoncon. If the value of the method option is None, then is returned. Otherwise will return the result of NIntegrate called with the method option.

For the cosine integral

the transformation corresponding to (44) is

Generalized Integrals

Here is the symbolic computation of the regularized divergent integral .
In[110]:=
Click for copyable input
Out[110]=
computes the nonregularized integral above in a generalized sense.
In[111]:=
Click for copyable input
Out[111]=

More about the properties of for divergent Fourier-type integrals can found in [MoriOoura93].

Nonalgebraic Multiplicand

Symbolic integration of an oscillatory integral.
In[116]:=
Click for copyable input
Out[116]=
If the oscillatory kernel is multiplied by a nonalgebraic function, still gives a good result.
In[117]:=
Click for copyable input
Out[117]=
Plots of the integrand and its oscillatory kernel.
In[119]:=
Click for copyable input
Out[119]=

Crude Monte Carlo and Quasi Monte Carlo Strategies

The crude Monte Carlo algorithm estimates a given integral by averaging integrand values over uniformly distributed random points in the integral's region. The number of points is incremented until the estimated standard deviation is small enough to satisfy the specified precision or accuracy goals. A Monte Carlo algorithm is called a quasi Monte Carlo algorithm if it uses equidistributed, deterministically generated sequences of points instead of uniformly distributed random points.

Here is a crude Monte Carlo integration.
In[3]:=
Click for copyable input
Out[3]=
Here is a crude quasi Monte Carlo integration.
In[4]:=
Click for copyable input
Out[4]=
option name
default value
Method"MonteCarloRule"Monte Carlo rule specification
MaxPoints50000maximum number of sampling points
"RandomSeed"Automatica seed to reset the random generator
"Partitioning"1partitioning of the integration region along each axis
"SymbolicProcessing"0number of seconds to do symbolic preprocessing

options.

option name
default value
MaxPoints50000maximum number of sampling points
"Partitioning"1partitioning of the integration region along each axis
"SymbolicProcessing"0number of seconds to do symbolic preprocessing

options.

In Monte Carlo methods [KrUeb98] the -dimensional integral is interpreted as the following expected (mean) value:

where is the mean value (the expectation) of the function interpreted as a random variable, with respect to the uniform distribution on , that is, the distribution with probability density vol(V)-1Boole(xV). Boole(xV) denotes the characteristic function of the region , while denotes the volume of .

The crude Monte Carlo estimate is made with the integration rule . The formulas for the integral and error estimation are given in "MonteCarloRule" in "NIntegrate Integration Rules".

Consider the integral

If the original integration region is partitioned into the set of disjoint subregions , , then the integral estimate is

and integration error is

The number of sampling points used on each subregion generally can be different, but in the Monte Carlo algorithms all are equal ().

The partitioning is called stratification, and each is called strata. Stratification can be used to improve crude Monte Carlo estimations. (The adaptive Monte Carlo algorithm uses recursive stratification.)

AccuracyGoal and PrecisionGoal

The default values for AccuracyGoal and PrecisionGoal are Infinity and 2, respectively, when the Monte Carlo algorithms of NIntegrate are used.

MaxPoints

The option specifies the maximum number of (pseudo) random sampling points to be used to compute the Monte Carlo estimate of an integral.

Here is an example in which the maximum number of sampling points is reached and NIntegrate stops with a message.
In[261]:=
Click for copyable input
Out[261]=

"RandomSeed"

The value of the option is used to seed the random generator used to make the sampling integration points. In that respect the use of in the Monte Carlo method is similar to the use of SeedRandom and RandomReal.

By using the results of a Monte Carlo integration can be reproduced. The results of the following two runs are identical.

Here is a Monte Carlo integration that uses .
In[56]:=
Click for copyable input
Out[56]//InputForm=
This Monte Carlo integration gives the same number.
In[57]:=
Click for copyable input
Out[57]//InputForm=
The following shows the first 20 points used in the Monte Carlo integrations.
In[65]:=
Click for copyable input
Out[66]=
The points coincide with the points made using SeedRandom and Random.
In[67]:=
Click for copyable input
Out[67]=

Stratified Crude Monte Carlo Integration

In stratified sampling Monte Carlo integration you break the region into several subregions and apply the crude Monte Carlo estimate on each subregion separately.

From the expected (mean) value formula, equation (45) at the beginning of "Crude Monte Carlo and Quasi Monte Carlo Strategies", you have

Let the region be bisected into two half-regions, and . is the expectation of on , and is the variance of on . From the theorem [PrFlTeuk92]

you can see that the stratified sampling gives a variance that is never larger than the crude Monte Carlo sampling variance.

There are two ways to specify strata for the strategy. One is to specify "singular" points in the variable range specifications, the other is to use the method suboption .

Stratified crude Monte Carlo integration using variable ranges specifications.
In[124]:=
Click for copyable input
Out[124]=
Stratified crude Monte Carlo integration using the suboption .
In[123]:=
Click for copyable input
Out[123]=

If is given a list of integers, with length that equals the number of integral variables, each dimension of the integration region is divided into equal parts. If is given an integer , all dimensions are divided into equal parts.

This graph demonstrates the stratified sampling specified with . Each cell contains 3 points, as specified by the option .
In[95]:=
Click for copyable input
Out[100]=

Stratified Monte Carlo sampling can be specified if the integration variable ranges are given with intermediate singular points.

Stratified Monte Carlo sampling through specification of intermediate singular points.
In[18]:=
Click for copyable input
Out[23]=

Stratified sampling improves the efficiency of the crude Monte Carlo estimation: if the number of strata is , the standard deviation of the stratified Monte Carlo estimation is times less of the standard deviation of the crude Monte Carlo estimation. (See the following example.)

The following benchmark shows that stratification speeds up the convergence.
In[120]:=
Click for copyable input
Out[121]=

Convergence Speedup of the Stratified Monte Carlo Integration

The following example confirms that if the number of strata is , the standard deviation of the stratified Monte Carlo estimation is times less than the standard deviation of the crude Monte Carlo estimation.

Here is a stratified integration definition based on the expected (mean) value formula (46).
In[122]:=
Click for copyable input
Consider this integral.
In[123]:=
Click for copyable input
Out[124]=
Here the integral above is approximated with 1000 points for the number of strata running from 1 to 40.
In[125]:=
Click for copyable input
These are the ratios between the standard deviations and the nonstratified crude Monte Carlo estimation.
In[126]:=
Click for copyable input
Note that is the ratio for the Monte Carlo estimation with i number of strata. This allows you to try a least-squares fit of the function to ratios.
In[127]:=
Click for copyable input
Out[128]=
The fitting of shows a coefficient very close to 1, which is a confirmation of the rule of thumb that number of strata gives -times faster convergence. This is the plot of the ratios and the least-squares fit.
In[130]:=
Click for copyable input
Out[130]=

Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies

The global adaptive Monte Carlo and quasi Monte Carlo strategies recursively bisect the subregion with the largest variance estimate into two halves, and compute integral and variance estimates for each half.

Here is an example of adaptive Monte Carlo integration.
In[1]:=
Click for copyable input
Out[1]=
option name
default value
MethodMonteCarloRule specification
"Partitioning"Automaticinitial partitioning of the integration region along each axis
"BisectionDithering"0offset from the middle of the region side that is parallel to the bisection axis
"MaxPoints"Automaticmaximum number of (pseudo) random sampling points to be used
"RandomSeed"Automaticrandom seed used to generate the (pseudo) random sampling points

Adaptive (quasi) Monte Carlo uses crude (quasi) Monte Carlo estimation rule on each subregion.

The process of subregion bisection and subsequent bi-integration is expected to reduce the global variance, and it is referred to as recursive stratified sampling. It is motivated by a theorem that states that if a region is partitioned into disjoint subregions the random variable variance over the region is greater than or equal to the sum of the random variable variances over each subregion. (See "Stratified Monte Carlo Integration" in "Crude Monte Carlo and Quasi Monte Carlo Strategies".)

The global adaptive Monte Carlo strategy is similar to . There are some important differences though.

1.   does not use singularity flattening, and does not have detectors for slow convergence and noisy integration.

2.   chooses randomly the bisection dimension. To avoid irregular separation of different coordinates a dimension recurs only if other dimensions have been chosen for bisection.

3.   can be tuned to bisect the subregions away from the middle. More at "BisectionDithering".

MinRecursion and MaxRecursion

The options and MaxRecursion for the adaptive Monte Carlo methods have the same meaning and functionality as they do for . See "MinRecursion and MaxRecursion".

"Partitioning"

The option of provides initial stratification of the integration. It has the same meaning and functionality as of the strategy .

"BisectionDithering"

When the integrand has some special symmetry that puts significant parts of it in the middle of the region, it is better if the bisection is done slightly away from the middle. The value of the option "BisectionDithering"->dith specifies that the splitting fraction of the region's splitting dimension side should be at instead of . The sign of dith is changed in a random manner. The default value given to is . The allowed values for dith are reals in the interval .

Consider the function.
In[195]:=
Click for copyable input
Out[196]=
Consider the integral.
In[197]:=
Click for copyable input
Out[197]=
Out[198]=
The integral is seriously underestimated if no bisection dithering is used, i.e. is given 0.
In[199]:=
Click for copyable input
Out[199]=
The following picture shows why the integral is underestimated. The black points are the integration sampling points. It can be seen that half of the peak of the integrand is undersampled.
In[204]:=
Click for copyable input
Out[207]=
Specifying bisection dithering of 10 percent gives a satisfactory estimation.
In[212]:=
Click for copyable input
Out[212]=
It can be seen on this plot that the peak of the integrand is sampled better.
In[213]:=
Click for copyable input
Out[216]=

Choice of Bisection Axis

For multidimensional integrals the adaptive Monte Carlo algorithm chooses the splitting axis of an integration region in two ways: (1) by random selection; or (2) by minimizing the variance of the integral estimates of each half. The axis selection is a responsibility of the "MonteCarloRule".

Example: Comparison with Crude Monte Carlo

Generally, the strategy has greater performance than . This is demonstrated with the examples in this subsection.

Consider the function.
In[217]:=
Click for copyable input
This is a plot of the function.
In[218]:=
Click for copyable input
Out[218]=

It can be seen from the following profiling that uses nearly three times fewer sampling points than the crude strategy.

These are the sampling points and timing for .
In[219]:=
Click for copyable input
Out[219]=
These are the sampling points and timing for .
In[220]:=
Click for copyable input
Out[220]=
This is the exact result.
In[221]:=
Click for copyable input
Out[221]=
Here is the timing for 100 integrations with .
In[222]:=
Click for copyable input
Out[222]=
The integration compares well with the exact result. The numbers below show the error of the mean of the integral estimates, the mean of the relative errors of the integral estimates, and the variance of the integral estimates.
In[223]:=
Click for copyable input
Out[223]=
Here is the timing for 100 integrations with , which is several times faster than integrations.
In[233]:=
Click for copyable input
Out[233]=
The integration result compares well with the exact result. The numbers below show the error of the mean of the integral estimates, the mean of the relative errors of the integral estimates, and the variance of the integral estimates.
In[234]:=
Click for copyable input
Out[234]=

"MultiPeriodic"

The strategy transforms all integrals into integrals over the unit cube and periodizes the integrands to be one periodic with respect to each integration variable. Different periodizing functions (or none) can be applied to different variables. works for integrals with dimension less than or equal to 12. If is given for integrals with higher dimension, the strategy is used.

Here is an example of integration with .
In[2]:=
Click for copyable input
Out[2]=
option name
default value
"Transformation"SidiSinperiodizing transformation applied to the integrand
"MinPoints"0minimal number of sampling points
"MaxPoints"105maximum number of sampling points
"SymbolicProcessing"Automaticnumber of seconds to be used for symbolic preprocessing

can be seen as a multidimensional generalization of the strategy . It can also be seen as a quasi Monte Carlo method.

uses lattice integration rules; see [SloanJoe94] [KrUeb98].

Here integration lattice in , , is understood to be a discrete subset of that is closed under addition and subtraction, and that contains . A lattice integration rule [SloanJoe94] is a rule of the form

where are all the points of an integration lattice contained in .

If is called on, a -dimensional integral option takes a list of one-argument functions that is used to transform the corresponding variables. If is given a list with length smaller than , then the last function, , is used for the last integration variables. If is given a function, that function will be used to transform all the variables.

Let be the dimension of the integral. If , calls after applying the periodizing transformation. For dimensions higher than 12, is called without applying periodizing transformations. uses the so-called copy rules for . For each , has a set of copy rules that are used to compute a sequence of integral estimates. The rules with a smaller number of points are used first. If the error estimate of a rule satisfies the precision goal, or if the difference of two integral estimates in the sequence satisfies the precision goal, the integration stops.

Number of points for the copy rules in the rule sets for different dimensions.
In[3]:=
Click for copyable input
Out[7]=

Comparison with "MultidimensionalRule"

Generally is slower than using . For computing high-dimension integrals with lower precision, might give results faster.

This defines the function of eight arguments.
In[8]:=
Click for copyable input
Timing in seconds for computing using and with .
In[11]:=
Click for copyable input
Out[12]//TableForm=
Number of integrand evaluations for computing using and with .
In[13]:=
Click for copyable input
Out[13]//TableForm=

Preprocessors

The capabilities of all strategies are extended through symbolic preprocessing of the integrals. The preprocessors can be seen as strategies that delegate integration to other strategies (preprocessors included).

"SymbolicPiecewiseSubdivision"

is a preprocessor that divides an integral with a piecewise integrand into integrals with disjoint integration regions on each of which the integrand is not piecewise.

option name
default value
MethodAutomaticintegration strategy or preprocessor to which the integration will be passed
"ExpandSpecialPiecewise"Truewhich piecewise functions should be expanded
TimeConstraint5the maximum number of seconds for which the piecewise subdivision will be attempted
"MaxPiecewiseCases"100the maximum number of subregions the piecewise preprocessor can return
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

Options of .

As was mentioned at the beginning of the tutorial, NIntegrate is able to integrate simultaneously integrals with disjoint domains, each having a different integrand. Hence, after the preprocessing with the integration continues in the same way as if, say, NIntegrate were given ranges with singularity specifications (which can be seen as specifying integrals with disjoint domains with the same integrand). For example, the strategy tries to improve the integral estimate of the region with the largest error through bisection, and will choose that largest error region regardless of which integrand it corresponds to.

Below are the sampling points for the numerical estimation of the integral . On the plot, the integrand is sampled at the coordinates in the order of the coordinates. It can be seen that alternates sampling for the piece with sampling for the piece .
In[12]:=
Click for copyable input
Out[13]=
Here are the sampling points for the numerical estimation of the integral . The integrand is plotted on the left, and the sampling points are plotted on the right. The integral has been partitioned into + + + , which is why the sampling points form a different pattern for .
In[14]:=
Click for copyable input
Out[17]=

"ExpandSpecialPiecewise"

In some cases it is preferable to do piecewise expansion only over certain piecewise functions. In these cases the option can be given a list of functions with which to do the piecewise expansion.

This Monte Carlo integral is done faster with piecewise expansion only over Boole.
In[18]:=
Click for copyable input
Out[19]=
Here is a Monte Carlo integration with piecewise expansion over both Boole and Abs.
In[20]:=
Click for copyable input
Out[20]=

"EvenOddSubdivision"

is a preprocessor that reduces the integration region if the region is symmetric around the origin and the integrand is determined to be even or odd. The convergence of odd integrals is verified by default.

option name
default value
MethodAutomaticintegration strategy or preprocessor to which the integration will be passed
VerifyConvergenceAutomaticwhether the convergence should be verified if an odd integral is detected
"SymbolicProcessing"Automaticnumber of seconds to do symbolic preprocessing

Options of .

When the integrand is an even function and the integration region is symmetric around the origin, the integral can be computed by integrating only on some part of the integration region and multiplying with a corresponding factor.

Here is a plot of an even function and the sampling points without any preprocessing.
In[21]:=
Click for copyable input
Out[24]=
These are the sampling points used by NIntegrate after has been applied. Note that the sampling points are only in the region .
In[25]:=
Click for copyable input
Out[25]=

Transformation Theorem

The preprocessor is based on the following theorem.

Theorem: Given the -dimensional integral

assume that for some these equalities hold:

a) ,

b) for all :

In other words the range of is symmetric around the origin, and the boundaries of the variables are even functions with respect to .

Then:

a) the integral is equivalent to

if the integrand is even with respect to , that is,

b) the integral is equivalent to 0 if the integrand is odd with respect to , that is,

Note that the theorem above can be applied several times over an integral.

To illustrate the theorem consider the integral . It is symmetric along , and the integrand and the bounds of are even functions with respect to .

Here is a plot of the sampling points without the application of (black) and with applied (red).
In[26]:=
Click for copyable input
Out[28]=

If the bounds of are not even functions with respect to then the symmetry along is broken. For example, the integral has no symmetry NIntegrate can exploit.

Here is a plot of the sampling points with applied (red). The region has no symmetry along .
In[29]:=
Click for copyable input
Out[30]=

"VerifyConvergence"

Consider the following divergent integral . NIntegrate detects it as an odd function over a symmetric domain and tries to integrate (that is, check the convergence of ). Since no convergence was reached as is indicated by the message, NIntegrate gives the message that the integral might be divergent.
In[31]:=
Click for copyable input
Out[31]=
If the option VerifyConvergence is set to False no convergence verification—and hence no integrand evaluation—will be done after the integral is found to be odd.
In[32]:=
Click for copyable input
Out[32]=

"OscillatorySelection"

is a preprocessor that selects specialized algorithms for efficient evaluation of oscillatory integrals.

On each integration region, detects whether the integrand is of Fourier finite-range type, Fourier infinite-range type (see "DoubleExponentialOscillatory"), Bessel infinite-range type (see "ExtrapolatingOscillatory"), Levin type (see "LevinRule"), or no special oscillatory type. Options control which method is selected for each type of integrand.

option name
default value
"BesselInfiniteRangeMethod"Automaticmethod for infinite-range Bessel integrals
"FourierFiniteRangeMethod"Automaticmethod for finite-range Fourier integrals
"FourierInfiniteRangeMethod"Automaticmethod for infinite-range Fourier integrals
"LevinMethod"Automaticmethod for Levin-type oscillatory integrals
Method"GlobalAdaptive"method for non-oscillatory integrals
"TermwiseOscillatory"Falsewhether to separately process terms in a sum
"SymbolicProcessing"5number of seconds to do symbolic processing

method options.

The preprocessor is used by default in NIntegrate.
In[98]:=
Click for copyable input
Out[98]=
Without the preprocessor NIntegrate does not reach convergence with its default option settings.
In[97]:=
Click for copyable input
Out[97]=

The preprocessor is designed to work with the internal output of the preprocessor. itself partitions oscillatory integrals that include the origin or have oscillatory kernels that need to be expanded or transformed into forms for which the oscillatory algorithms are designed.

Here is a piecewise function integration in which all of the specialized integration strategies for one-dimensional integrals are automatically selected by . For this integral the preprocessor first divides the integral into four different integrals; for each of these integrals selects an appropriate specialized algorithm.
In[96]:=
Click for copyable input
Out[96]=
The following table shows the names of the options used to specify the algorithms for each subinterval in the integral above.
x(-∞,0]"BesselInfiniteRangeMethod"
x[0,20]"FourierFiniteRangeMethod"
x[30,∞)"FourierInfiniteRangeMethod"
x[20,30]Method
In this example is called twice. is a special algorithm for Fourier integrals, and the formula makes the integrand a sum of two Fourier integrands.
In[99]:=
Click for copyable input
Out[99]//InputForm=
To demonstrate that has used the formula , here is the integral above split "by hand". The result is identical with the last result.
In[100]:=
Click for copyable input
Out[100]//InputForm=

The value Automatic for the option means that if the integration strategy specified with the option Method is one of or then that strategy will be used for the finite-range Fourier integration, otherwise will be used.

Here is a piecewise function integration that uses strategy for the non-oscillatory integral and for the finite-range oscillatory integral.
In[37]:=
Click for copyable input
Out[37]=
These are the sampling points of the preceding integration and integral but with default option settings. The pattern between on the left plot is typical for the local adaptive quadrature—the recursive partitioning into three parts can be seen (because of the option given to ). The pattern over on the right plot comes from . The pattern between on the first plot is typical for the double exponential quadrature. The same pattern can be seen on the second plot between since uses by default the singularity handler.
In[106]:=
Click for copyable input
Out[110]=

If the application of a particular oscillatory method is desired for a particular type of oscillatory integrals, either the corresponding options of should be changed, or the Method option in NIntegrate should be used without the preprocessor .

Here is a piecewise function integration that uses for any of the infinite-range oscillatory integrals.
In[111]:=
Click for copyable input
Out[111]=
If is given as the method, uses it for infinite-range oscillatory integration.
In[112]:=
Click for copyable input
Out[112]=
The integration above is faster with the default options of NIntegrate. For this integral , which is applied by default, uses .
In[113]:=
Click for copyable input
Out[113]=

"UnitCubeRescaling"

is a preprocessor that transforms the integration region into a unit cube or hypercube. The variables of the original integrand are replaced and the result is multiplied by the Jacobian of the transformation.

option name
default value
"FunctionalRangesOnly"Truewhat ranges should be transformed to the unit cube
Method"GlobalAdaptive"integration strategy or preprocessor to which the integration will be passed
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

Options of .

This uses unit cube rescaling and it is faster than the computation that follows.
In[10]:=
Click for copyable input
Out[10]=
This integration does not use unit cube rescaling. It is done approximately three times slower than the previous one.
In[11]:=
Click for copyable input
Out[11]=

transforms the integral

into an integral over the hypercube . Assuming that and are finite and are piecewise continuous functions, the transformation used by is

The Jacobian of this transformation is

If for the ^(th) axis one or both of and are infinite, then the formula for in (47) is a non-affine transformation that maps into . NIntegrate uses the following transformations:

where .

Applying makes the integrand more complicated if the integration region boundaries are constants (finite or infinite). Since NIntegrate has efficient affine and infinite internal variable transformations the integration process would become slower. If some of the integration region boundaries are functions, applying would make the integration faster since the computations that involve the integration variables are done only when the integrand is evaluated. Because of these performance considerations, has the option . If is set to True the rescaling is applied only to multidimensional functional ranges.

This integration uses unit cube rescaling.
In[12]:=
Click for copyable input
Out[12]=
This integration does not use unit cube rescaling. It is done approximately two times faster than the previous one.
In[13]:=
Click for copyable input
Out[13]=

Example Implementation

The transformation process used by is the same as the following one implemented by the function (also defined in "Duffy's Coordinates Generalization and Example Implementation").

This function provides the transformation (48) and its Jacobian (49) for a list of integration ranges and a list of rectangular parallelepiped sides or a hypercube side.
In[14]:=
Click for copyable input

Each transformation of the transformation (50) can be done with Rescale.

In[17]:=
Click for copyable input
Out[17]=

Note that for given axis the transformation rules already derived for axes need to be applied to the original boundaries before the rescaling of boundaries along the ^(th) axis.

The transformation rules and the Jacobian for .
In[18]:=
Click for copyable input
Out[19]=
Out[20]=
Application of the transformation to a function.
In[21]:=
Click for copyable input
Out[21]=
The transformation rules and the Jacobian for .
In[22]:=
Click for copyable input
Out[23]=
Out[24]=
The transformation rules and the Jacobian for .
In[25]:=
Click for copyable input
Out[26]=
Out[27]=

"SymbolicPreprocessing"

is a composite preprocessor made to simplify the switching on and off of the other preprocessors.

option name
default value
MethodAutomaticintegration strategy or preprocessor to which the integration will be passed
"SymbolicPiecewiseSubdivision"Truepiecewise subdivision
"EvenOddSubdivision"Trueeven-odd subdivision
"OscillatorySelection"Truedetection of products with an oscillatory function
"UnitCubeRescaling"Automaticrescaling to the unit hypercube
"SymbolicProcessing"Automaticnumber of seconds to do symbolic processing

options.

When is set to Automatic it is applied only to multidimensional functional ranges.

Here is an example of the integration of with different combinations of preprocessor application.
In[30]:=
Click for copyable input
Out[44]=

Examples and Applications

Closed-Contour Integrals

This function calculates the mass of a closed contour given in polar coordinates parametrization.
In[42]:=
Click for copyable input
This is the circumference of the ellipse with radii 2 and 3 using Integrate.
In[43]:=
Click for copyable input
Out[44]=
Here is the circumference approximation of the ellipse with radii 2 and 3 using the same function.
In[45]:=
Click for copyable input
Out[45]=
The approximation compares quite well with the exact value.
In[46]:=
Click for copyable input
Out[46]=

Fourier Series Calculation

This is a Mathematica function that calculates a truncated Fourier series approximation of a function.
In[83]:=
Click for copyable input
Fourier approximation of over using Integrate.
In[84]:=
Click for copyable input
Out[84]=
This a plot of and the Fourier series approximation.
In[85]:=
Click for copyable input
Out[85]=
This calculates a 60-term Fourier approximation of Sin[x3+] over using NIntegrate. If Integrate is used the calculation will be very slow.
In[86]:=
Click for copyable input
Out[86]=
This a plot of Sin[x3+] and the Fourier series approximation.
In[87]:=
Click for copyable input
Out[87]=
New to Mathematica? Find your learning path »
Have a question? Ask support »