This chapter provides methods for the numerical evaluation of definite integrals in one or more dimensions and for evaluating weights and abscissae of integration rules.
Syntax
C# 

public static class D01 
Visual Basic 

Public NotInheritable Class D01 
Visual C++ 

public ref class D01 abstract sealed 
F# 

[<AbstractClassAttribute>] [<SealedAttribute>] type D01 = class end 
Background to the Problems
The methods in this chapter are designed to estimate:
(a)  the value of a onedimensional definite integral of the form
Some methods are specially designed for integrands of the form


(b)  the values of the onedimensional indefinite integrals arising from (1) where the ranges of integration are interior to the interval $\left[a,b\right]$.  
(c)  the value of a multidimensional definite integral of the form
The simplest form of ${R}_{n}$ is the $n$rectangle defined by

Onedimensional Integrals
To estimate the value of a onedimensional integral, a quadrature rule uses an approximation in the form of a weighted sum of integrand values, i.e.,
The points ${x}_{i}$ within the interval $\left[a,b\right]$ are known as the abscissae, and the ${w}_{i}$ are known as the weights.
$$\underset{a}{\overset{b}{\int}}f\left(x\right)dx\simeq \sum _{i=1}^{N}{w}_{i}f\left({x}_{i}\right)\text{.}$$  (5) 
More generally, if the integrand has the form (2), the corresponding formula is
If the integrand is known only at a fixed set of points, these points must be used as the abscissae, and the weighted sum is calculated using finite difference methods. However, if the functional form of the integrand is known, so that its value at any abscissa is easily obtained, then a wide variety of quadrature rules are available, each characterised by its choice of abscissae and the corresponding weights.
$$\underset{a}{\overset{b}{\int}}w\left(x\right)g\left(x\right)dx\simeq \sum _{i=1}^{N}{w}_{i}g\left({x}_{i}\right)\text{.}$$  (6) 
The appropriate rule to use will depend on the interval $\left[a,b\right]$ – whether finite or otherwise – and on the form of any $w\left(x\right)$ factor in the integrand. A suitable value of $N$ depends on the general behaviour of $f\left(x\right)$; or of $g\left(x\right)$, if there is a $w\left(x\right)$ factor present.
Among possible rules, we mention particularly the Gaussian formulae, which employ a distribution of abscissae which is optimal for $f\left(x\right)$ or $g\left(x\right)$ of polynomial form.
The choice of basic rules constitutes one of the principles on which methods for onedimensional integrals may be classified. The other major basis of classification is the implementation strategy, of which some types are now presented.
(a)  Single rule evaluation procedures
A fixed number of abscissae, $N$, is used. This number and the particular rule chosen uniquely determine the weights and abscissae. No estimate is made of the accuracy of the result. 

(b)  Automatic procedures
The number of abscissae, $N$, within $\left[a,b\right]$ is gradually increased until consistency is achieved to within a level of accuracy (absolute or relative) you requested. There are essentially two ways of doing this; hybrid forms of these two methods are also possible:

An ideal generalpurpose method would be an automatic method which could be used for a wide variety of integrands, was efficient (i.e., required the use of as few abscissae as possible), and was reliable (i.e., always gave results to within the requested accuracy). Complete reliability is unobtainable, and generally higher reliability is obtained at the expense of efficiency, and vice versa. It must therefore be emphasized that the automatic methods in this chapter cannot be assumed to be $100\%$ reliable. In general, however, the reliability is very high.
Multidimensional Integrals
A distinction must be made between cases of moderately low dimensionality (say, up to $4$ or $5$ dimensions), and those of higher dimensionality. Where the number of dimensions is limited, a onedimensional method may be applied to each dimension, according to some suitable strategy, and high accuracy may be obtainable (using product rules). However, the number of integrand evaluations rises very rapidly with the number of dimensions, so that the accuracy obtainable with an acceptable amount of computational labour is limited; for example a product of $3$point rules in $20$ dimensions would require more than ${10}^{9}$ integrand evaluations. Special techniques such as the Monte–Carlo methods can be used to deal with high dimensions.
(a)  Products of onedimensional rules
Using a twodimensional integral as an example, we have
A different onedimensional rule may be used for each dimension, as appropriate to the range and any weight function present, and a different strategy may be used, as appropriate to the integrand behaviour as a function of each independent variable.
For a ruleevaluation strategy in all dimensions, the formula (8) is applied in a straightforward manner. For automatic strategies (i.e., attempting to attain a requested accuracy), there is a problem in deciding what accuracy must be requested in the inner integral(s). Reference to formula (7) shows that the presence of a limited but random error in the $y$integration for different values of ${x}_{i}$ can produce a ‘jagged’ function of $x$, which may be difficult to integrate to the desired accuracy and for this reason products of automatic onedimensional methods should be used with caution (see Lyness (1983)). 

(b)  Monte–Carlo methods
These are based on estimating the mean value of the integrand sampled at points chosen from an appropriate statistical distribution function. Usually a variance reducing procedure is incorporated to combat the fundamentally slow rate of convergence of the rudimentary form of the technique. These methods can be effective by comparison with alternative methods when the integrand contains singularities or is erratic in some way, but they are of quite limited accuracy. 

(c)  Number theoretic methods
These are based on the work of Korobov and Conroy and operate by exploiting implicitly the properties of the Fourier expansion of the integrand. Special rules, constructed from socalled optimal coefficients, give a particularly uniform distribution of the points throughout $n$dimensional space and from their number theoretic properties minimize the error on a prescribed class of integrals. The method can be combined with the Monte–Carlo procedure. 

(d)  Sag–Szekeres method
By transformation this method seeks to induce properties into the integrand which make it accurately integrable by the trapezoidal rule. The transformation also allows effective control over the number of integrand evaluations. 

(e)  Sparse grid methods
Given a set of onedimensional quadrature rules of increasing levels of accuracy, the sparse grid method constructs an approximation to a multidimensional integral using ${n}_{\mathit{dim}}$ dimensional tensor products of the differences between rules of adjacent levels. This provides a lower theoretical accuracy than the methods in (a), the full grid approach, however this is still sufficient for various classes of sufficiently smooth integrands. Furthermore, it requries substantially fewer evaluations than the full grid approach. Specifically, if a onedimensional quadrature rule has $N\sim \mathit{O}\left({2}^{\ell}\right)$ points, the full grid will require $\mathit{O}\left({N}^{{n}_{\mathit{dim}}}\right)$ function evaluations, whereas the sparse grid of level $\ell $ will require $\mathit{O}\left({2}^{\ell}{{n}_{\mathit{dim}}}^{\ell 1}\right)$. Hence a sparse grid approach is computationally feasible even for integrals over ${n}^{\mathit{dim}}\sim \mathit{O}\left(100\right)$.
Sparse grid methods are deterministic, and may be viewed as automatic whole domain procedures if their level $\ell $ is allowed to increase. 

(f)  Automatic adaptive procedures
An automatic adaptive strategy in several dimensions normally involves division of the region into subregions, concentrating the divisions in those parts of the region where the integrand is worst behaved. It is difficult to arrange with any generality for variable limits in the inner integral(s). For this reason, some methods use a region where all the limits are constants; this is called a hyperrectangle. Integrals over regions defined by variable or infinite limits may be handled by transformation to a hyperrectangle. Integrals over regions so irregular that such a transformation is not feasible may be handled by surrounding the region by an appropriate hyperrectangle and defining the integrand to be zero outside the desired region. Such a technique should always be followed by a Monte–Carlo method for integration.
The method used locally in each subregion produced by the adaptive subdivision process is usually one of three types: Monte–Carlo, number theoretic or deterministic. Deterministic methods are usually the most rapidly convergent but are often expensive to use for high dimensionality and not as robust as the other techniques. 
References
Davis P J and Rabinowitz P (1975) Methods of Numerical Integration Academic Press
Gonnet P (2010) Increasing the reliability of adaptive quadrature using explicit interpolants ACM Trans. Math. software 37 26
Lyness J N (1983) When not to use an automatic quadrature routine SIAM Rev. 25 63–87
Patterson T N L (1968) The Optimum addition of points to quadrature formulae Math. Comput. 22 847–856
Piessens R, de Doncker–Kapenga E, Überhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration Springer–Verlag
Sobol I M (1974) The Monte Carlo Method The University of Chicago Press
Stroud A H (1971) Approximate Calculation of Multiple Integrals Prentice–Hall
Wynn P (1956) On a device for computing the ${e}_{m}\left({S}_{n}\right)$ transformation Math. Tables Aids Comput. 10 91–96