Introduction to Numerical Integration in the Wolfram Language

OverviewStrategies, Rules, and Preprocessors

Overview

The Wolfram Language function NIntegrate is a general numerical integrator. It can handle a wide range of one-dimensional and multidimensional integrals.

NIntegrate[f[x1,x2,…,xn],{x1,a1,b1},{x2,a2,b2},…,{xn,an,bn}]
find a numerical integral for the function f over the region [a1,b1]×[a2,b2]×…×[an,bn]

Finding a numerical integral of a function over a region.

In general, NIntegrate estimates the integral through sampling of the integrand value over the integration region. The various numerical integration methods prescribe the initial sampling steps and how the sampling evolves.

This numerically computes the integral :
This plots the sequence of values of at which the integrand is evaluated. The samples are more dense near , where the integrand varies more rapidly:

NIntegrate uses algorithms called "integration strategies" that attempt to compute integral estimates that satisfy user-specified precision or accuracy goals. The integration strategies use "integration rules" that compute a single integral estimate from a set of integrand values, often using a weighted sum.

This computes a numerical integral using a specific integration strategy and rule:

Symbolic Preprocessing

NIntegrate uses symbolic preprocessing to simplify integrals with special structure and to automatically select integration methods. Integrands that are even or odd functions or that contain piecewise functions may lead to the integration region being transformed or separated into multiple distinct integration regions.

This integrates a piecewise function over the interval . The integration region is automatically separated at the singularity :

Highly oscillatory integrands are identified and specialized integration rules are applied.

This integrates a highly oscillatory function over the interval . A specialized method suited to Fourier sine integrals is automatically selected:
This is a plot of the previous oscillatory integrand over of the integration region:

Symbolic preprocessing allows the automatic computation of a wide variety of integrals containing discontinuities and regions of extremely rapid variation.

This integrates a piecewise function that is highly oscillatory in one region. The integration region is automatically separated at and and a specialized Fourier method is applied to the oscillatory region :

Quadrature Rules

NIntegrate includes most classical one-dimensional quadrature rules.

Here is a one-dimensional integral performed with several different standard quadrature rules:

The classical quadrature rules work by forming a linear combination of the sampled integrand values. The vector of weights in the linear combination is fixed for each quadrature rule.

These are the nodes and weights for a three-point Newton–Cotes quadrature rule on the unit interval :

For multidimensional integrals, NIntegrate includes a class of rules based on sparse grids and also allows rules formed from the Cartesian product of one-dimensional rules.

The default method "MultidimensionalRule" is a sparse grid rule:
Cartesian product rules can be specified as a list of one-dimensional rules, one for each dimension:
This shows the sampling points used in the first step of each of the above two multidimensional rules:

Boole can be used to specify more complicated multidimensional regions. Regions specified this way may also be further simplified during symbolic preprocessing.

Here is a two-dimensional function: a cone with base inscribed in the square ×:
Here is the numerical integral of the cone function:
Here are the sampling points used by NIntegrate. Note that the sampling points are only in one quadrant of the integration region, and all sampling points lie inside the region :
Here are the sampling points used by NIntegrate without symbolic preprocessing. The adaptive algorithm must sample the entire integration region. Also note the dense sampling around the boundary of the region . Without symbolic preprocessing, NIntegrate must effectively numerically locate this boundary:

Oscillatory Integration

NIntegrate contains general oscillatory integration methods applicable to a very wide range of integrands, over finite or infinite regions, and in either one dimension or multiple dimensions. Additionally, NIntegrate contains several methods that are specifically suited to one-dimensional integrals of functions of particular forms involving Exp, trigonometric functions such as Sin and Cos, and certain other special functions such as BesselJ.

This computes a numerical integral of an irregularly oscillatory function:
This is a plot of the above irregularly oscillatory function over a finite range:

Automatic Singularity Handling

NIntegrate has several ways to deal with singular integrands. The deterministic adaptive strategies "GlobalAdaptive" and "LocalAdaptive" use singularity handling techniques (based on variable transformations) to speed up the convergence of the integration process. The strategy "DoubleExponential" employs trapezoidal quadrature with a special variable transformation on the integrand. This rule-transformation combination achieves optimal convergence for integrands analytic on an open set in the complex plane containing the interval of integration.

Here is a one-dimensional integration with singularity handling:
You can perform the integration 100 times to get more precise timing information:
Without singularity handling the integral is computed more slowly:

Special Strategies

The strategy "DuffyCoordinates" simplifies or eliminates certain types of singularities in multidimensional integrals. Integrals with certain spherical symmetry can converge very quickly.

Here is an integration using "DuffyCoordinates":
Here is the same integral using default settings for NIntegrate. The correct result is obtained, but the computation time is greater:

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

Here is a calculation of an integral computed with the trapezoidal strategy. The result is compared with the exact value. The result computed with "Trapezoidal" is obtained faster and it is more precise than the one with the default NIntegrate settings:
Here is a (slower) computation of the same integral but with the default Method settings for NIntegrate:

For high-dimensional integrals, or in cases when only a rough integral estimate is needed, Monte Carlo methods are useful. NIntegrate has both crude and adaptive Monte Carlo and quasi Monte Carlo strategies.

Here is a high-dimensional integral computed quickly using an adaptive Monte Carlo algorithm:

Design

The principal features of the NIntegrate framework are:

  • Object orientation (method property specification and communication)
  • Separation of method initialization phase and runtime computation
  • Hierarchical and reentrant numerical methods
  • Type- and precision-dynamic methods
  • User extensibility and prototyping through plugin capabilities
  • Specialized data structures

Strategies, Rules, and Preprocessors

NIntegrate integration strategies can be classified according to how they sample the integration region, the class of integrands to which they can be applied, and whether they are "rule-based" 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.

Adaptive sampling strategies try to improve the integral estimate by sampling more often in subregions with a larger error estimate, typically by subdividing those subregions. Uniform sampling strategies try to improve the integral estimate by uniformly increasing the density of sampling throughout the whole integration region.

Rule-based strategies apply a given integration rule to each subregion to obtain integral and error estimates for that region. The integration rule can be specified with the setting Method->{"strategy",Method->"rule"}.

Here is how to specify an integration using globally adaptive subdivision with Clenshaw–Curtis quadrature on each subregion:

NIntegrate integration rules can be classified according to whether they apply to one-dimensional or multidimensional regions, and according to the type of integration rule.

"BooleRule"one-dimensional, weighted sum
"ClenshawCurtisRule"one-dimensional, weighted sum
"GaussBerntsenEspelidRule"one-dimensional, weighted sum
"GaussKronrodRule"one-dimensional, weighted sum
"LobattoKronrodRule"one-dimensional, weighted sum
"LobattoPeanoRule"one-dimensional, weighted sum
"NewtonCotesRule"one-dimensional, weighted sum
"TrapezoidalRule"one-dimensional, weighted sum
"ClenshawCurtisOscillatoryRule"one-dimensional, specialized oscillatory rule
"LevinRule"one- or multidimensional, general oscillatory rule
"MultipanelRule"one-dimensional, weighted sum, combination of 1D rules
"CartesianRule"multidimensional, weighted sum, product of 1D rules
"MultidimensionalRule"multidimensional, weighted sum

Integration rules that can be used with the rule-based strategies "GlobalAdaptive" and "LocalAdaptive".

Classical "weighted sum"-type rules estimate the integral as a predetermined linear combination of the function values at a set of points. Oscillatory rules estimate the integral using quadrature weights that depend on the particular oscillatory "kernel" of the integrand.

Combination rules construct a quadrature rule from one or more subrules. They are specified with the setting Method->{"rule",Method->{"subrule1",…}}.

Here is how to specify the Cartesian product multidimensional rule of two one-dimensional subrules:

The capabilities of all strategies are extended through symbolic preprocessing of the integrand. Preprocessing is controlled by preprocessor strategies that first transform or analyze the integral, then delegate integration to another strategy (often another preprocessor strategy).

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

Preprocessor strategies are specified with the setting Method->{"preprocessor",Method->m}, where m is the strategy or rule to which the integration is delegated after preprocessing is complete.

Here is how to explicitly apply a preprocessor strategy to an integrand that is even with respect to both of its variables. After preprocessing, the integration is delegated to the "LocalAdaptive" strategy:

Preprocessor strategies often reduce the amount of work required by the final integration strategy.

Here are the sampling points of the previous integration. Without preprocessing, four times as many sampling points would be required, covering the entire integration region:

User Extensibility

Built-in methods can be used as building blocks for the efficient construction of special-purpose integrators. User-defined integration rules, integration strategies, and preprocessor strategies can also be added.