NIntegrate Integration Strategies

Introduction

In the following integration, the symbolic piecewise preprocessor in NIntegrate recognizes the integrand as a piecewise function, and the integration is done over regions for which x≥1 with the integrand and regions for which x≤1 with the integrand . In[31]:= Out[31]=

Here is a plot of all sampling points used in the integration. The integrand is sampled at the x coordinates in the order of the y coordinates (in the plot). It can be seen that the sampling points are concentrated near the singularity point 1. The patterns formed by the sampling points at the upper part of the plot differ from the patterns of the lower part of the plot because a singularity handler is applied. In[10]:= Out[11]=

Adaptive Strategies

Global Adaptive Strategy

The global adaptive algorithm for NIntegrate is specified with the Method option value "GlobalAdaptive". In[32]:= Out[32]=

"GlobalAdaptive" options.

MinRecursion and MaxRecursion

"MaxErrorIncreases"

Zero integral. In[3]:= Out[3]=

The working precision is not dense enough. In[33]:= Out[33]//InputForm=

The following integration invokes the message NIntegrate::eincr, with the default value of "MaxErrorIncreases". In[1]:= Out[1]=

Increasing "MaxErrorIncreases" silences the NIntegrate::eincr message. In[2]:= Out[2]=

The result compares well with the exact value. In[3]:= Out[4]=

Example Implementation of a Global Adaptive Strategy

This computes Gauss-Kronrod abscissas, weights, and error weights. In[15]:=

This is a definition of a function that applies the integration rule with abscissas and weights computed to the function f over the interval {a, b}. In[16]:=

This is a definition of a simple global adaptive algorithm that finds the integral of the function f over the interval {aArg, bArg} with relative error tol. In[17]:=

This defines an integrand. In[18]:=

The global adaptive strategy defined earlier gives the following result. In[19]:= Out[19]=

Here is the exact result. In[20]:= Out[20]=

The relative error is within the prescribed tolerance. In[21]:= Out[21]=

Local Adaptive Strategy

The local adaptive algorithm for NIntegrate is specified with the Method option value "LocalAdaptive". In[5]:= Out[5]=

"LocalAdaptive" options.

MinRecursion and MaxRecursion

"InitialEstimateRelaxation"

"Partitioning"

Consider the following function. In[4]:= Out[4]=

These are the sampling points used by "LocalAdaptive" with its automatic region partitioning. It can be seen that the sampling points of each recursion level are between the sampling points of the previous recursion level. In[1]:= Out[2]=

These are the sampling points used by "LocalAdaptive" integration which partitions the regions with large error into three subregions. The patterns formed clearly show the three next recursion level subregions of each region of the first and second recursion levels. In[5]:= Out[6]=

Multidimensional example of using the "Partitioning" option. To make the plot, the sampling points of the first region to be integrated, [0, 1]×[0, 1], are removed. In[7]:= Out[10]=

Reuse of Integrand Values

Sampling points of the integration of by "LocalAdaptive". The variable rulePoints determines the number of points in the integration rule used by "LocalAdaptive". In[13]:= Out[15]= Out[16]=

The percent of reused points in the integration. In[17]:= Out[19]=

Example Implementation of a Local Adaptive Strategy

This computes Clenshaw-Curtis abscissas, weights, and error weights. In[33]:=

This is a definition of a function that applies the integration rule, with the abscissas and weights computed in the previous example, to the function f over the interval {a, b}. In[34]:=

This defines a simple local adaptive algorithm that finds the integral of the function f over the interval {aArg, bArg} with relative error tol. In[35]:=

This defines a function. In[37]:=

The local adaptive strategy gives the result. In[38]:= Out[38]=

This is the exact result. In[39]:= Out[39]=

The relative error is within the prescribed tolerance. In[40]:= Out[40]=

"GlobalAdaptive" versus "LocalAdaptive"

"GlobalAdaptive" Is Generally Better than "LocalAdaptive"

Here are the definitions of some functions, precision goals, number of integrations, and the integration rule. The variable integrationRule can be changed in order to compare the profiling runs with the same integration rule. The last function is derived from -xsin (x) by the variable change that transforms [0, 1) into [0, ). In[70]:= In[74]:= Out[74]=

Exact integral values. All integrals are over [0, 1]. In[75]:=

"GlobalAdaptive" timings. In[76]:=

"LocalAdaptive" timings. In[77]:=

"GlobalAdaptive" function evaluations. In[78]:=

"LocalAdaptive" function evaluations. In[79]:=

Table with the timing ratios and integrand evaluations. In[80]:= Out[80]=

Table with the errors of the integrations. Both "GlobalAdaptive" and "LocalAdaptive" reach the required precision goals. In[81]:= Out[81]=

Singularity Handling

User-Specified Singularities

Point Singularities

Here is an example of an integral that has two singular points at and . In[58]:= Out[58]=

Here is an example of a two-dimensional integral with a singular point at (1, 1). In[59]:= Out[59]=

Here is an example of an integral that has two singular points at and specified with the Exclusions option. In[60]:= Out[60]=

Here is an example of a two-dimensional integral with a singular point at (1, 1) specified with the Exclusions option. In[61]:= Out[61]=

Curve, Surface, and Hypersurface Singularities

This two-dimensional function is singular along the curve x2+y2=1. In[62]:= Out[62]=

Using Exclusions the integral is quickly calculated. In[12]:= Out[12]=

NIntegrate will reach convergence much more slowly if no singularity specification is given. In[35]:= Out[35]=

Here is an example of a case in which a singular curve can be specified with the variable ranges. If x[0, 2] and y[0, 2] this would not be possible—see the following example. In[10]:= Out[10]=

Example Implementation of Curve, Surface, and Hypersurface Singularity Handling

This two-dimensional function has singular points along the curve x2+y2=1. In[66]:= Out[66]=

Using Boole the integral is calculated quickly. In[9]:= Out[9]=

This two-dimensional function has singular points along the curve x+ (1-y)2=1. In[68]:= Out[68]=

Again, using Boole the integral is calculated quickly. In[8]:= Out[8]=

This is how the sampling points of the integration look. In[6]:= 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]:=

This defines a three-dimensional function. In[2]:=

Here is the integral of a three-dimensional function with singular points along the surface x+ (1-y)2+ (1-z)2=1. In[3]:= Out[3]=

These are the sampling points of the integration. In[4]:= Out[5]=

"SingularityHandler" and "SingularityDepth"

"GlobalAdaptive" and "LocalAdaptive" singularity handling options.

Use of the IMT Variable Transformation

Here is a numerical integration that uses the IMT variable transformation for singularity handling. In[13]:= Out[13]=

IMT singularity handler options.

The IMT transformation a, p (t): (0, 1]→ (0, 1], a>0, p>0, is defined. In[14]:=

The limit of the derivative of the IMT transformation is 0. In[16]:= Out[16]=

Here is the plot of the IMT transformation. In[17]:= Out[17]=

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]:= Out[19]=

IMT Transformation by Example

Consider the function over (0, 1] that has a singularity at 0. In[29]:= In[30]:= Out[30]=

Assume the integration is done with "GlobalAdaptive", with singularity handler IMT and singularity depth 4. After four bisections "GlobalAdaptive" will have a region with boundaries {0, 1/16} that contains the singular end point. For that region the IMT variable transformation will change its boundaries to {0, 1} and its integrand to the following. In[31]:= Out[32]=

Here is the plot of the new integrand. In[33]:= Out[34]=

For example, consider the sampling points and weights of the Gauss-Kronrod rule. In[35]:=

The Gauss-Kronrod sampling points for the region {0, 1/16} and the derivatives of the rescaling follow. In[36]:= Out[36]= In[37]:= Out[37]=

Here is the integral estimate. In[38]:= Out[38]=

With the IMT transformation, these are the sampling points and derivatives. In[39]:= Out[39]= In[40]:= Out[40]=

Here is the integral estimate with the IMT transformation. In[41]:= Out[41]=

The estimate calculated with the IMT variable transformation is much closer to the exact value. In[42]:= Out[42]=

Use of Double-Exponential Quadrature

Here is a numerical integration that uses double-exponential quadrature for singularity handling. In[103]:= Out[103]=

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

This loads a package that defines the profiling function NIntegrateProfile that gives the number of sampling points and the time needed by a numerical integration command. In[17]:=

Table for a "very" analytical integrand that the "DoubleExponential" singularity handler easily computes. In[34]:= 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 "DoubleExponential" and then switches to IMT. In[37]:= Out[40]//TableForm=

A table for an integrand, , for which the Automatic singularity handler starts with "DoubleExponential" and then switches to IMT. In[41]:= Out[44]//TableForm=

IMT Multidimensional Singularity Handling

The number of integrand evaluations and timings for an integrand that has a singularity only along the x axis. The default (automatic) singularity handler chooses to apply IMT to regions obtained after the default (four) bisections. In[19]:= Out[19]=

The number of integrand evaluations and timings for an integrand that has a singularity only along the x axis with no singularity handler application. In[20]:= 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]:= 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]:= 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]:= Out[23]=

Duffy's Coordinates for Multidimensional Singularity Handling

The following integration uses Duffy's coordinates. In[63]:= Out[63]=

The following integration does not use Duffy's coordinates. In[62]:= Out[62]=

Here is the number of sampling points for the simpler integral. In[58]:= Out[58]=

Here is the number of sampling points for the more complicated integral. In[59]:= Out[59]=

Here is a table that compares the different singularity handlings for . (The profiling function NIntegrateProfile defined earlier is used.) In[74]:= Out[75]=

Duffy's Coordinates Strategy

Here is an example with an integrand that has singularities at two different corners of its integration region. In[84]:= Out[84]=

"DuffyCoordinates" options.

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

Here is an example for the "DuffyCoordinates" integrand for a two-dimensional function that is singular at two of the corners of the integration region. In[79]:= Out[79]=

Integration with "DuffyCoordinates". In[80]:= Out[80]=

Integration with the default NIntegrate options settings which is much slower than the previous one. In[81]:= Out[81]=

Here is another example of a speedup by "DuffyCoordinates". In[82]:= Out[82]=

Integration with the default NIntegrate options settings which is much slower than the previous one. In[83]:= Out[83]=

Duffy's Coordinates Generalization and Example Implementation

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

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]:=

Here is an example of unit-square to infinite region rescaling. In[96]:= 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]:=

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]:=

Here is a symbolic example. In[101]:= Out[101]=

Here is another symbolic example. In[102]:= Out[102]=

Here is a computational example. In[103]:= Out[103]=

Using Duffy's coordinates is much faster than using no singularity handling (see the next example). In[108]:= Out[110]=

Integration using no singularity handling. In[111]:= Out[111]=

Of course, the internal implementation of NIntegrate gives similar performance results. In[107]:= Out[107]=

Ignoring the Singularity

Here is a plot of the integrand. In[118]:= Out[118]=

NIntegrate gives a result that is close to the exact one. In[114]:= Out[115]=

Convergence is achieved when MaxRecursion is increased. In[45]:= Out[45]=

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

Sampling points that belong to the interval [1-106, 1+106]. In[120]:= Out[120]=

Integration with the singular point specified. In[123]:= Out[123]=

Integration by ignoring the singularity. In[122]:= Out[122]=

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

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

For example, these two integrals are equivalent. In[126]:= Out[126]=

NIntegrate can do the first one. In[127]:= Out[127]=

NIntegrate cannot do the second one. In[128]:= Out[128]=

Automatic Singularity Handling

One-Dimensional Integration

Here is an example of a switch from "DoubleExponential" to IMT singularity handling. On the plot the integrand is sampled at the x coordinates in the order of the y coordinates. The patterns of the sampling points over show the change from Gaussian quadrature (y[0, 97]) to double-exponential quadrature (y[98, 160]), which later is replaced by Gaussian quadrature using the IMT variable transformation (y[160, 400]). In[143]:= Out[145]=

Multidimensional Integration

Sampling points of two-dimensional integration, , with Automatic (left) and "DuffyCoordinates" (right) singularity handling. It can be seen that the automatic singularity handling uses almost two times more points than "DuffyCoordinates". To illustrate the effect of the singularity handlers they are applied after two bisections. In[133]:= Out[134]=

Timings for the integral, , with singularity handlers Automatic, "DuffyCoordinates", and IMT and with no singularity handling. The integral has a point singularity at {0, 0}. In[47]:= Out[47]//TableForm=

Timings for the integral, , singular along y axis with singularity handlers Automatic, "DuffyCoordinates", and IMT and with no singularity handling. In[46]:= Out[46]//TableForm=

Cauchy Principal Value Integration

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

"PrincipalValue" options.

This finds the Cauchy principal value of . In[14]:= Out[14]=

Here is the Cauchy principal value of . Note that there are two singularities that need to be specified. In[114]:= Out[114]=

The singular points can be specified using the Exclusions option. In[30]:= Out[30]=

This checks the value. The result would be 0 if everything were done exactly. In[31]:= Out[31]=

Sampling Points Visualization

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]:= Out[156]=

This defines a function which accumulates the argument values given to the integrand. In[1]:=

Here are the points at which the integrand has been evaluated. Note the symmetric pattern over the interval . In[166]:= Out[168]=

Double-Exponential Strategy

The double-exponential algorithm for NIntegrate is specified with the Method option value "DoubleExponential". In[169]:= Out[169]=

"DoubleExponential" options.

Cartesian double-exponential quadrature. In[48]:= Out[48]=

The following plot illustrates the Cartesian product character of the "DoubleExponential" multidimensional integration. In[49]:= Out[50]=

MinRecursion and MaxRecursion

Comparison of Double-Exponential and Gaussian Quadrature

Integration with "DoubleExponential". In[215]:= Out[215]=

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

Double-exponential quadrature and Gaussian quadrature for . Increasing the precision goal does not change the number of sampling points used by "DoubleExponential". In[217]:= Out[219]=

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

"DoubleExponential" needs more than 10000 integrand evaluations to compute this integral with a non-analytic integrand. In[52]:= Out[53]=

Gaussian quadrature is much faster for the integral. In[54]:= Out[54]=

Here is an example with a not "very" analytical integrand. In[226]:= Out[226]=

Again, Gaussian quadrature is much faster. In[227]:= Out[227]=

Here are the plots of the integrand and its derivative. In[228]:= Out[228]=

Convergence Rate

This defines a double-exponential integration function that an returns integral estimate and the number of points used. In[229]:=

This defines a function. In[230]:=

This is the exact integral. In[231]:= Out[231]=

This finds the errors and number of evaluation points for a range of step sizes of the trapezoidal rule. In[232]:=

This fits through the logarithms of the errors; see (19). In[239]:=

Here is the fitted function. The summation term 30.48 is just a translation parameter. In[240]:= Out[240]=

You see that the errors fit the model (20): In[241]:= Out[241]=

Example Implementation of Double-Exponential Quadrature

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]:=

This is a definition of a simple double-exponential strategy, which finds the integral of the function f over the finite interval {a, b} with relative error tol. In[189]:=

This defines a function that is singular at 0. In[194]:=

Here is the integral estimate from the double-exponential strategy. In[195]:= Out[195]//InputForm=

Here is the exact result. In[177]:= Out[177]//InputForm=

This defines an oscillating function. In[178]:=

Here is the integral estimate given by the double-exponential strategy. In[179]:= Out[179]//InputForm=

Here is the exact result. In[180]:= Out[180]=

Here is the exact result in machine precision. In[181]:= Out[181]//InputForm=

The relative error is within the prescribed tolerance. In[182]:= Out[182]=

"Trapezoidal" Strategy

"Trapezoidal" options.

In[109]:=

Out[109]=

Here is a table that shows the number of sampling points for different values of the parameter t used by "GlobalAdaptive" and "Trapezoidal" respectively for the integral . In[33]:= Out[35]//TableForm=

Example Implementation

This function makes a trapezoidal quadrature integral estimate with specified points. In[242]:=

This function improves a trapezoidal quadrature integral estimate using sampling points between the old ones. In[257]:=

This function is an interface to the preceding one. In[272]:=

Here is a definition of a (Bessel) function. In[269]:=

Here is the trapezoidal quadrature estimate. In[274]:= Out[274]=

Here is the exact value. In[278]:= Out[278]=

The relative error is within the prescribed tolerance. In[279]:= Out[279]=

Oscillatory Strategies

Here is an example that uses both algorithms. In[13]:= Out[13]=

Finite Region Oscillatory Integration

Number of integrand evaluations for modified Clenshaw-Curtis quadrature for slowly oscillating kernel. In[1]:= Out[1]=

Timing and integral estimates for modified Clenshaw-Curtis quadrature for slowly oscillating kernel. In[3]:= Out[3]=

Number of integrand evaluations for modified Clenshaw-Curtis quadrature for highly oscillating kernel. In[5]:= Out[5]=

Timing and integral estimates for modified Clenshaw-Curtis quadrature for highly oscillating kernel. In[6]:= Out[6]=

Number of integrand evaluations for Gaussian quadrature for slowly oscillating kernel. In[7]:= Out[7]=

Timing and integral estimates for Gaussian quadrature for slowly oscillating kernel. In[8]:= Out[8]=

Number of integrand evaluations for Gaussian quadrature for highly oscillating kernel. In[9]:= Out[9]=

Timing and integral estimates for Gaussian quadrature for highly oscillating kernel. In[10]:= Out[10]=

Extrapolating Oscillatory Strategy

Here is an example of an integration using "ExtrapolatingOscillatory". In[294]:= Out[294]=

Example Implementation

Here is a definition of an oscillation function that will be integrated in the interval [0, ). The zeros of the oscillating function sin ( x) are . In[1]:=

Here is a plot of the oscillatory function in the interval [0, 10]. In[89]:= Out[89]=

This is a definition of a function that integrates between two consequent zeros. The zeros of the oscillating function k (x)=sin ( x) are . In[5]:=

Here is the integral estimate computed by sequence convergence acceleration (extrapolation). In[6]:= Out[6]=

Here is the exact integral value. In[7]:= Out[7]=

The integral estimate is very close to the exact value. In[8]:= Out[8]=

Here is another check using the "ExtrapolatingOscillatory" strategy. In[94]:= Out[94]=

The integral estimate by "ExtrapolatingOscillatory" is very close to the exact value. In[95]:= Out[95]=

Double-Exponential Oscillatory Integration

Integration with "DoubleExponentialOscillatory". In[2]:= Out[2]=

Options of "DoubleExponentialOscillatory".

Generalized Integrals

Here is the symbolic computation of the regularized divergent integral . In[110]:= Out[110]=

"DoubleExponentialOscillatory" computes the nonregularized integral above in a generalized sense. In[111]:= Out[111]=

Non-Algebraic Multiplicand

Symbolic integration of an oscillatory integral. In[116]:= Out[116]=

If the oscillatory kernel is multiplied by a nonalgebraic function, "DoubleExponentialOscillatory" still gives a good result. In[117]:= Out[117]=

Plots of the integrand and its oscillatory kernel. In[119]:= Out[119]=

Crude Monte Carlo and Quasi Monte Carlo Strategies

Here is a crude Monte Carlo integration. In[3]:= Out[3]=

Here is a crude quasi Monte Carlo integration. In[4]:= Out[4]=

"MonteCarlo" options.

"QuasiMonteCarlo" options.

AccuracyGoal and PrecisionGoal

MaxPoints

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

"RandomSeed"

Here is a Monte Carlo integration that uses "RandomSeed". In[56]:= Out[56]//InputForm=

This Monte Carlo integration gives the same number. In[57]:= Out[57]//InputForm=

The following shows the first 20 points used in the Monte Carlo integrations. In[65]:= Out[66]=

The points coincide with the points made using SeedRandom and Random. In[67]:= Out[67]=

Stratified Crude Monte Carlo Integration

Stratified crude Monte Carlo integration using variable ranges specifications. In[124]:= Out[124]=

Stratified crude Monte Carlo integration using the sub-option "Partitioning". In[123]:= Out[123]=

This graph demonstrates the stratified sampling specified with "Partitioning". Each cell contains 3 points, as specified by the "MonteCarloRule" option "Points". In[95]:= Out[100]=

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

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

Convergence Speedup of the Stratified Monte Carlo Integration

Here is a stratified integration definition based on the expected (mean) value formula (46). In[122]:=

Consider this integral. In[123]:= Out[124]=

Here the integral above is approximated with 1000 points for the number of strata running from 1 to 40. In[125]:=

These are the ratios between the standard deviations and the nonstratified, crude Monte Carlo estimation. In[126]:=

Note that ratiosi 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]:= Out[128]=

The fitting of shows a coefficient very close to 1, which is a confirmation of the rule of thumb that s number of strata give s-times faster convergence. This is the plot of the ratios and the least squares fit. In[130]:= Out[130]=

Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies

Here is an example of adaptive Monte Carlo integration. In[1]:= Out[1]=

MinRecursion and MaxRecursion

"Partitioning"

"BisectionDithering"

Consider the function. In[195]:= Out[196]=

Consider the integral. In[197]:= Out[197]= Out[198]=

The integral is seriously underestimated if no bisection dithering is used i.e., "BisectionDithering" is given 0. In[199]:= 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]:= Out[207]=

Specifying bisection dithering of 10% gives a satisfactory estimation. In[212]:= Out[212]=

It can be seen on this plot, that the peak of the integrand is sampled better. In[213]:= Out[216]=

Choice of Bisection Axis

Example: Comparison with Crude Monte Carlo

Consider the function. In[217]:=

This is a plot of the function. In[218]:= Out[218]=

These are the sampling points and timing for "MonteCarlo". In[219]:= Out[219]=

These are the sampling points and timing for "AdaptiveMonteCarlo". In[220]:= Out[220]=

This is the exact result. In[221]:= Out[221]=

Here is the timing for 100 integrations with "MonteCarlo". In[222]:= Out[222]=

The "MonteCarlo" 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]:= Out[223]=

Here is the timing for 100 integrations with "AdaptiveMonteCarlo", which is several times faster than "MonteCarlo" integrations. In[233]:= Out[233]=

The "AdaptiveMonteCarlo" 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]:= Out[234]=

"MultiPeriodic"

Here is an example of integration with "MultiPeriodic". In[2]:= Out[2]=

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

Comparison with "MultiDimensionalRule"

This defines the function of eight arguments. In[8]:=

Timing in seconds for computing using "MultiPeriodic" and "GlobalAdaptive" with "MultiDimensionalRule". In[11]:= Out[12]//TableForm=

Number of integrand evaluations for computing using "MultiPeriodic" and "GlobalAdaptive" with "MultiDimensionalRule". In[13]:= Out[13]//TableForm=

Preprocessors

"SymbolicPiecewiseSubdivision"

Options of "SymbolicPiecewiseSubdivision".

Below are the sampling points for the numerical estimation of the integral . On the plot, the integrand is sampled at the x coordinates in the order of the ord coordinates. It can be seen that "GlobalAdaptive" alternates sampling for the piece with sampling for the piece . In[12]:= Out[13]=

Here are the sampling points for the numerical estimation of the integral . The integrand is plotted on the left, the sampling points are plotted on the right. The integral has been partitioned into + + + , that is why the sampling points form a different pattern for -1≤ x≤1. In[14]:= Out[17]=

"ExpandSpecialPiecewise"

This Monte Carlo integral is done faster with piecewise expansion only over Boole. In[18]:= Out[19]=

Here is a Monte Carlo integration with piecewise expansion over both Boole and Abs. In[20]:= Out[20]=

"EvenOddSubdivision"

Options of "EvenOddSubdivision".

Here is a plot of an even function and the sampling points without any preprocessing. In[21]:= Out[24]=

These are the sampling points used by NIntegrate after "EvenOddSubdivision" has been applied. Note that the sampling points are only in the region [0, ]×[0, ]. In[25]:= Out[25]=

Transformation Theorem

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

Here is a plot of the sampling points with "EvenOddSubdivision" applied (red). The region has no symmetry along y. In[29]:= 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 ncvb message, NIntegrate gives the message oidiv that the integral might be divergent. In[31]:= 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]:= Out[32]=

"OscillatorySelection"

Options of "OscillatorySelection".

"OscillatorySelection" is used by default. In[33]:= Out[33]=

Without the "OscillatorySelection" preprocessor NIntegrate does not reach convergence with its default option settings. In[34]:= Out[34]=

Here is a piecewise function integration in which all methods of "OscillatorySelection" are used. For this integral the preprocessor "SymbolicPiecewiseSubdivision" divides the integral into four different integrals; for each of these integrals "OscillatorySelection" selects an appropriate algorithm. In[1]:= Out[1]=

In this example "DoubleExponentialOscillatory" is called twice. "DoubleExponentialOscillatory" is a special algorithm for Fourier integrals, and the formula 2 x2=cos (2x2)+ sin (2x2) makes the integrand a sum of two Fourier integrands. In[35]:= Out[35]//InputForm=

To demonstrate that "OscillatorySelection" has used the formula 2 x2=cos (2x2)+ sin (2x2), here is the integral above split "by hand." The result is identical with the last result. In[36]:= Out[36]//InputForm=

Here is a piecewise function integration that uses "DoubleExponential" strategy for the non-oscillatory integral and "LocalAdaptive" for the finite range oscillatory integral. In[37]:= Out[37]=

These are the sampling points of the preceding integration and integral but with default option settings. The pattern between [0, 20] on the left picture is typical for the local adaptive quadrature—the recursive partitioning into three parts can be seen (because of the option "Partitioning"->3 given to "LocalAdaptive"). The pattern over [0, 20] on the right picture comes from "GlobalAdaptive". The pattern between [20, 40] on the first picture is typical for the double-exponential quadrature. The same pattern can be seen on the second picture between [20, 21+1/4] since "GlobalAdaptive" uses by default the "DoubleExponential" singularity handler. In[38]:= Out[42]=

Here is a piecewise function integration that uses "ExtrapolatingOscillatory" for any of the infinite range oscillatory integrals. In[10]:= Out[10]=

If "ExtrapolatingOscillatory" is given as the method, "OscillatorySelection" uses it for infinite range oscillatory integration. In[1]:= Out[1]=

The integration above is faster with the default options of NIntegrate. For this integral "OscillatorySelection", which is applied by default, uses "DoubleExponentialOscillatory". In[2]:= Out[2]=

Working with Sums of Oscillating Terms

Here is an example of integration that applies the specialized oscillatory NIntegrate algorithms to each term of the integrand. In[4]:= Out[4]=

By default this option is set to False, and the integral cannot be computed. In[5]:= Out[5]=

Here is a convergent integral. If it is split into two integrals each will be divergent. In[6]:= Out[6]=

If "TermwiseOscillatory"->True is used the result is some big number (and lots of messages). In[8]:= Out[8]=

If "TermwiseOscillatory"->False is used (with increased MaxRecursion) the result is closer to the exact one. In[7]:= Out[7]=

"UnitCubeRescaling"

Options of "UnitCubeRescaling".

This uses unit cube rescaling and it is faster than the computation that follows. In[10]:= Out[10]=

This integration does not use unit cube rescaling. It is done approximately three times slower than the previous one. In[11]:= Out[11]=

This integration uses unit cube rescaling. In[12]:= Out[12]=

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

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]:=

In[17]:=

Out[17]=

The transformation rules and the Jacobian for [0, 1]×[0, 1]→ [0, 1]×[a (x), b (x)]. In[18]:= Out[19]= Out[20]=

Application of the transformation to a function. In[21]:= Out[21]=

The transformation rules and the Jacobian for [0, 1]×[0, 1]→ [0, ]×[a (x), b (x)]. In[22]:= Out[23]= Out[24]=

The transformation rules and the Jacobian for [0, 1]×[0, 1]→ [a0, b0]×[a1 (x), b1 (x)]. In[25]:= Out[26]= Out[27]=

"SymbolicPreprocessing"

"SymbolicPreprocessing" options.

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

Examples and Applications

Closed-Contour Integrals

This function calculates the mass of a closed contour given in polar coordinates parametrization. In[42]:=

This is circumference of the ellipse with radii 2 and 3 using Integrate. In[43]:= Out[44]=

Here is the circumference approximation of the ellipse with radii 2 and 3 using the same function. In[45]:= Out[45]=

The approximation compares quite well with the exact value. In[46]:= Out[46]=

Fourier Series Calculation

This is a Mathematica function that calculates a truncated Fourier series approximation of a function. In[83]:=

Fourier approximation of x3+x2 over [-2, 2] using Integrate. In[84]:= Out[84]=

This a plot of x3+x2 and the Fourier series approximation. In[85]:= Out[85]=

This calculates a 60-term Fourier approximation of over [-4, 4] using NIntegrate. If Integrate is used the calculation will be very slow. In[86]:= Out[86]=

This a plot of and the Fourier series approximation. In[87]:= Out[87]=