NAG Library Chapter Introduction
D01 – Quadrature
1 Scope of the Chapter
This chapter provides routines for the numerical evaluation of definite integrals in one or more dimensions and for evaluating weights and abscissae of integration rules.
2 Background to the Problems
The routines in this chapter are designed to estimate:
(a) 
the value of a onedimensional definite integral of the form
where $f\left(x\right)$ is defined by you, either at a set of points $\left({x}_{\mathit{i}},f\left({x}_{\mathit{i}}\right)\right)$, for $\mathit{i}=1,2,\dots ,n$, where $a={x}_{1}<{x}_{2}<\cdots <{x}_{n}=b$, or in the form of a function; and the limits of integration $a,b$ may be finite or infinite.
Some methods are specially designed for integrands of the form
which contain a factor $w\left(x\right)$, called the weightfunction, of a specific form. These methods take full account of any peculiar behaviour attributable to the $w\left(x\right)$ factor. 
(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
where $f\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$ is a function defined by you and ${R}_{n}$ is some region of $n$dimensional space.
The simplest form of ${R}_{n}$ is the $n$rectangle defined by
where ${a}_{i}$ and ${b}_{i}$ are constants. When ${a}_{i}$ and ${b}_{i}$ are functions of ${x}_{j}$ ( $j<i$), the region can easily be transformed to the rectangular form (see page 266 of Davis and Rabinowitz (1975)). Some of the methods described incorporate the transformation procedure. 
2.1 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.
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.
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:
(i) 
whole interval procedures (nonadaptive)
A series of rules using increasing values of $N$ are successively applied over the whole interval $\left[a,b\right]$. It is clearly more economical if abscissae already used for a lower value of $N$ can be used again as part of a higherorder formula. This principle is known as optimal extension. There is no overlap between the abscissae used in Gaussian formulae of different orders. However, the Kronrod formulae are designed to give an optimal $\left(2N+1\right)$point formula by adding $\left(N+1\right)$ points to an $N$point Gauss formula. Further extensions have been developed by Patterson. 
(ii) 
adaptive procedures
The interval $\left[a,b\right]$ is repeatedly divided into a number of subintervals, and integration rules are applied separately to each subinterval. Typically, the subdivision process will be carried further in the neighbourhood of a sharp peak in the integrand than where the curve is smooth. Thus, the distribution of abscissae is adapted to the shape of the integrand.
Subdivision raises the problem of what constitutes an acceptable accuracy in each subinterval. The usual global acceptability criterion demands that the sum of the absolute values of the error estimates in the subintervals should meet the conditions required of the error over the whole interval. Automatic extrapolation over several levels of subdivision may eliminate the effects of some types of singularities. 

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 routines in this chapter cannot be assumed to be $100\%$ reliable. In general, however, the reliability is very high.
2.2 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
where $\left({w}_{i},{x}_{i}\right)$ and $\left({v}_{i},{y}_{i}\right)$ are the weights and abscissae of the rules used in the respective dimensions.
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 routines 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 $d$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, which is nonetheless 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({2}^{\mathit{ld}}\right)$ function evaluations, whereas the sparse grid of level $\ell $ will require $\mathit{O}\left({2}^{\ell}{d}^{\ell 1}\right)$. Hence a sparse grid approach is computationally feasible even for integrals over $d\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. 
3 Recommendations on Choice and Use of Available Routines
This section is divided into five subsections. The first subsection illustrates the difference between direct and reverse communication routines. The second subsection highlights the different levels of vectorization provided by different interfaces.
Sections 3.3,
3.4 and
3.5 consider in turn routines for: onedimensional integrals over a finite interval, and over a semiinfinite or an infinite interval; and multidimensional integrals. Within each subsection, routines are classified by the type of method, which ranges from simple rule evaluation to automatic adaptive algorithms. The recommendations apply particularly when the primary objective is simply to compute the value of one or more integrals, and in these cases the automatic adaptive routines are generally the most convenient and reliable, although also the most expensive in computing time.
Note however that in some circumstances it may be counterproductive to use an automatic routine. If the results of the quadrature are to be used in turn as input to a further computation (e.g., an ‘outer’ quadrature or an optimization problem), then this further computation may be adversely affected by the ‘jagged performance profile’ of an automatic routine; a simple ruleevaluation routine may provide much better overall performance. For further guidance, the article by
Lyness (1983) is recommended.
3.1 Direct and Reverse Communication
Routines in this chapter which evaluate an integral value may be classified as either direct communication or reverse communication. See
Section 3.3.3 in How to Use the NAG Library and its Documentation for a description of these terms.
Currently in this chapter the only routine explicitly using reverse communication is
D01RAF.
3.2 Choice of Interface
This section concerns the design of the interface for the provision of abscissae, and the subsequent collection of calculated information, typically integrand evaluations. Vectorized interfaces typically allow for more efficient operation.
(a) 
Single abscissa interfaces
The algorithm will provide a single abscissa at which information is required. These are typically the most simple to use, although they may be significantly less efficient than a vectorized equivalent. Most of the algorithms in this chapter are of this type.
Examples of this include
D01AJF and D01FBF.

(b) 
Vectorized abscissae interfaces
The algorithm will return a set of abscissae, at all of which information is required. While these are more complicated to use, they are typically more efficient than a nonvectorized equivalent. They reduce the overhead of function calls, allow the avoidance of repetition of computations common to each of the integrand evaluations, and offer greater scope for vectorization and parallelization of your code.
Examples include
D01RGF, D01UAF, and the routines D01ATF and D01AUF, which are vectorized equivalents of D01AJF and D01AKF.

(c) 
Multiple integral interfaces
These are routines which allow for multiple integrals to be estimated simultaneously. As with (b) above, these are more complicated to use than single integral routines, however they can provide higher efficiency, particularly if several integrals require the same subcalculations at the same abscissae. They are most efficient if integrals which are supplied together are expected to have similar behaviour over the domain, particularly when the algorithm is adaptive.
Examples include D01EAF and D01RAF.

3.3 Onedimensional Integrals over a Finite Interval
(a) 
Integrand defined at a set of points
If $f\left(x\right)$ is defined numerically at four or more points, then the Gill–Miller finite difference method ( D01GAF) should be used. The interval of integration is taken to coincide with the range of $x$ values of the points supplied. It is in the nature of this problem that any routine may be unreliable. In order to check results independently and so as to provide an alternative technique you may fit the integrand by Chebyshev series using E02ADF and then use routine E02AJF to evaluate its integral (which need not be restricted to the range of the integration points, as is the case for D01GAF). A further alternative is to fit a cubic spline to the data using E02BAF and then to evaluate its integral using E02BDF. 
(b) 
Integrand defined as a function
If the functional form of $f\left(x\right)$ is known, then one of the following approaches should be taken. They are arranged in the order from most specific to most general, hence the first applicable procedure in the list will be the most efficient.
However, if you do not wish to make any assumptions about the integrand, the most reliable routines to use will be
D01ATF (or D01AJF), D01AUF (or D01AKF), D01ALF, D01RGF or D01RAF, although these will in general be less efficient for simple integrals.
(i) 
Ruleevaluation routines
If $f\left(x\right)$ is known to be sufficiently well behaved (more precisely, can be closely approximated by a polynomial of moderate degree), a Gaussian routine with a suitable number of abscissae may be used.
D01BCF or D01TBF
with D01FBF may be used if it is required to examine the weights and abscissae.
D01TBF
is faster and more accurate, whereas
D01BCF
is more general. D01UAF uses the same quadrature rules as D01TBF, and may be used if you do not explicitly require the weights and abscissae.
If $f\left(x\right)$ is well behaved, apart from a weightfunction of the form
D01BCF
with D01FBF may be used.
D01BCF and D01TBF
generate weights and abscissae for specific Gauss rules. Weights and abscissae for other quadrature formulae may be computed using routines D01TDF or D01TEF. Wherever possible use D01TDF in preference to D01TEF. The former however requires information that may not be readily available. 
(ii) 
Automatic wholeinterval routines
If $f\left(x\right)$ is reasonably smooth, and the required accuracy is not too high, the automatic whole interval
routines D01ARF and D01BDF
may be used. Additionally, D01ESF with $d=1$ may be used with an appropriate transformation from the unit interval.
D01BDF uses the Gauss $10$point rule, with the $21$ point Kronrod extension, and the subsequent $43$ and $87$ point Patterson extensions if required.
D01ESF supports multiple simultaneous integrals, and has a vectorized interface. Either high order Gauss–Patterson rules (of size ${2}^{\ell}1$, for $\ell =1,\dots ,9$), or high order ClenshawCurtis rules (of size ${2}^{\ell 1}+1$, for $\ell =2,\dots ,12$). Gauss–Patterson rules possess greater polynomial accuracy, whereas Clenshaw–Curtis rules are often well suited to oscillatory integrals.
D01ARF incorporates the same high order Gauss–Patterson rules as D01ESF, and is the only routine that may be used for indefinite integration. 
(iii) 
Automatic adaptive routines
Firstly, several routines are available for integrands of the form $w\left(x\right)g\left(x\right)$ where $g\left(x\right)$ is a ‘smooth’ function (i.e., has no singularities, sharp peaks or violent oscillations in the interval of integration) and $w\left(x\right)$ is a weight function of one of the following forms.
1. 
if $w\left(x\right)={\left(bx\right)}^{\alpha}{\left(xa\right)}^{\beta}{\left(\mathrm{log}\left(bx\right)\right)}^{k}{\left(\mathrm{log}\left(xa\right)\right)}^{l}$, where $k,l=0\text{ or}1$, $\alpha ,\beta >1$: use
D01APF; 
2. 
if $w\left(x\right)=\frac{1}{xc}$: use
D01AQF
(this integral is called the Hilbert transform of $g$); 
3. 
if $w\left(x\right)=\mathrm{cos}\left(\omega x\right)$ or $\mathrm{sin}\left(\omega x\right)$: use
D01ANF
(this routine can also handle certain types of singularities in $g\left(x\right)$). 
Secondly, there are multiple routines for general $f\left(x\right)$, using different strategies.
D01ATF (and D01AJF), and D01AUF (and D01AKF)
use the strategy of Piessens et al. (1983), using repeated bisection of the interval, and in the first case the $\epsilon $algorithm ( Wynn (1956)), to improve the integral estimate. This can cope with singularities away from the end points, provided singular points do not occur as abscissae,
D01AUF tends to perform better than D01ATF
on more oscillatory integrals.
D01ALF
uses the same subdivision strategy as
D01ATF
over a set of initial interval segments determined by supplied breakpoints. It is hence suitable for integrals with discontinuities (including switches in definition) or sharp peaks occuring at known points. Such integrals may also be approximated using other routines which do not allow breakpoints, although such integrals should be evaluated over each of the subintervals seperately.
D01RAF again uses the strategy of Piessens et al. (1983), and provides the functionality of
D01ALF, D01ATF and D01AUF
in a reverse communication framework. It also supports multiple integrals and uses a vectorized interface for the abscissae. Hence it is likely to be more efficient if several similar integrals are required to be evaluated over the same domain. Furthermore, its behaviour can be tailored through the use of optional parameters.
D01AHF uses the strategy of Patterson (1968) and the $\epsilon $algorithm to adaptively evaluate the integral in question. It tends to be more efficient than the bisection based algorithms, although these tend to be more robust when singularities occur away from the end points.
D01RGF uses another adaptive scheme due to Gonnet (2010). This attempts to match the quadrature rule to the underlying integrand as well as subdividing the domain. Further, it can explicitly deal with singular points at abscissae, should NaN's or ∞ be returned by the usersupplied (sub)routine, provided the generation of these does not cause the program to halt (see Chapter X07). 

3.4 Onedimensional Integrals over a Semiinfinite or Infinite Interval
(a) 
Integrand defined at a set of points
If $f\left(x\right)$ is defined numerically at four or more points, and the portion of the integral lying outside the range of the points supplied may be neglected, then the Gill–Miller finite difference method, D01GAF, should be used. 
(b) 
Integrand defined as a function
(i) 
Rule evaluation routines
If $f\left(x\right)$ behaves approximately like a polynomial in $x$, apart from a weight function of the form:
1. 
${e}^{\beta x},\beta >0$ (semiinfinite interval, lower limit finite); or 
2. 
${e}^{\beta x},\beta <0$ (semiinfinite interval, upper limit finite); or 
3. 
${e}^{\beta {\left(x\alpha \right)}^{2}},\beta >0$ (infinite interval), 
or if $f\left(x\right)$ behaves approximately like a polynomial in ${\left(x+b\right)}^{1}$ (semiinfinite range), then the Gaussian routines may be used.
D01UAF
may be used if it is not required to examine the weights and abscissae.
D01BCF or D01TBF
with D01FBF may be used if it is required to examine the weights and abscissae.
D01TBF
is faster and more accurate, whereas
D01BCF
is more general.
D01UBF returns an approximation to the specific problem $\underset{0}{\overset{\infty}{\int}}}\mathrm{exp}\left({x}^{2}\right)\text{}g\left(x\right)dx$. 
(ii) 
Automatic adaptive routines
D01AMF
may be used, except for integrands which decay slowly towards an infinite end point, and oscillate in sign over the entire range. For this class, it may be possible to calculate the integral by integrating between the zeros and invoking some extrapolation process (see C06BAF).
D01ASF
may be used for integrals involving weight functions of the form $\mathrm{cos}\left(\omega x\right)$ and $\mathrm{sin}\left(\omega x\right)$ over a semiinfinite interval (lower limit finite).
The following alternative procedures are mentioned for completeness, though their use will rarely be necessary.
1. 
If the integrand decays rapidly towards an infinite end point, a finite cutoff may be chosen, and the finite range methods applied. 
2. 
If the only irregularities occur in the finite part (apart from a singularity at the finite limit, with which
D01AMF
can cope), the range may be divided, with
D01AMF used on the infinite part. 
3. 
A transformation to finite range may be employed, e.g.,
will transform $\left(0,\infty \right)$ to $\left(1,0\right)$ while for infinite ranges we have
If the integrand behaves badly on $\left(\infty ,0\right)$ and well on $\left(0,\infty \right)$ or vice versa it is better to compute it as $\underset{\infty}{\overset{0}{\int}}}f\left(x\right)dx+{\displaystyle \underset{0}{\overset{\infty}{\int}}}f\left(x\right)dx$. This saves computing unnecessary function values in the semiinfinite range where the function is well behaved. 


3.5 Multidimensional Integrals
A number of techniques are available in this area and the choice depends to a large extent on the dimension and the required accuracy. It can be advantageous to use more than one technique as a confirmation of accuracy, particularly for highdimensional integrations. Several routines include a transformation procedure, using a usersupplied subroutine, which allows general product regions to be easily dealt with in terms of conversion to the standard $n$cube region.
(a) 
Products of onedimensional rules (suitable for up to about $5$ dimensions)
If $f\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$ is known to be a sufficiently well behaved function of each variable ${x}_{i}$, apart possibly from weight functions of the types provided, a product of Gaussian rules may be used. These are provided by
D01BCF or D01TBF
with D01FBF. Rules for finite, semiinfinite and infinite ranges are included.
For twodimensional integrals only, unless the integrand is very badly behaved, the automatic wholeinterval product procedure of D01DAF may be used. The limits of the inner integral may be userspecified functions of the outer variable. Infinite limits may be handled by transformation (see Section 3.4); end point singularities introduced by transformation should not be troublesome, as the integrand value will not be required on the boundary of the region.
If none of these routines proves suitable and convenient, the onedimensional routines may be used recursively. For example, the twodimensional integral
may be expressed as
The usersupplied code to evaluate $F\left(x\right)$ will call the integration routine for the $y$integration, which will call more usersupplied code for $f\left(x,y\right)$ as a function of $y$ ( $x$ being effectively a constant).
From Mark 24 onwards, all direct communication routines may be called recursively. As such, you may use any routine, including the same routine, for each dimension. Note however, in previous releases,
direct communication routines were not defined as recursive, and thus a different library integration routine must be used for each dimension if you are using an older product. Apart from this restriction, the following combinations were not permitted:
D01AJF and D01ALF, D01ANF and D01APF, D01APF and D01AQF, D01ANF and D01AQF, D01ANF and D01ASF, D01AMF and D01ASF, D01ATF and D01AUF.
Otherwise the full range of onedimensional routines are available, for finite/infinite intervals, constant/variable limits, rule evaluation/automatic strategies etc.
The reverse communication routine D01RAF may be used by itself in a pseudorecursive manner, in that it may be called to evaluate an inner integral for the integrand value of an outer integral also being calculated by D01RAF. 
(b) 
Sag–Szekeres method
Two routines are based on this method.
D01FDF is particularly suitable for integrals of very large dimension although the accuracy is generally not high. It allows integration over either the general product region (with builtin transformation to the $n$cube) or the $n$sphere. Although no error estimate is provided, two adjustable arguments may be varied for checking purposes or may be used to tune the algorithm to particular integrals.
D01JAF is also based on the Sag–Szekeres method and integrates over the $n$sphere. It uses improved transformations which may be varied according to the behaviour of the integrand. Although it can yield very accurate results it can only practically be employed for dimensions not exceeding $4$. 
(c) 
Number Theoretic method
Two subroutines are based on this method, D01GCF and a vectorized equivalent D01GDF.
Algorithms of this type carry out multidimensional integration using the Korobov–Conroy method over a product region with builtin transformation to the $n$cube. A stochastic modification of this method is incorporated into the routines in this library, hybridising the technique with the Monte–Carlo procedure. An error estimate is provided in terms of the statistical standard error. A number of precomputed optimal coefficient rules for up to $20$ dimensions are provided; others can be computed using D01GYF and D01GZF. Like the Sag–Szekeres method it is suitable for large dimensional integrals although the accuracy is not high.
D01GCF requires a function to be provided to evaluate the value of the integrand at a single abscissa, and a subroutine to return the upper and lower limits of integration in a given dimension.
D01GDF has a vectorized interface which can result in faster execution, especially on vectorprocessing machines. You are required to provide two subroutines, the first to return an array of values of the integrand at each of an array of points, and the second to evaluate the limits of integration at each of an array of points. This reduces the overhead of function calls, avoids repetitions of computations common to each of the evaluations of the integral and limits of integration, and offers greater scope for vectorization of your code. 
(d) 
A combinatorial extrapolation method
D01PAF computes a sequence of approximations and an error estimate to the integral of a function over a multidimensional simplex using a combinatorial method with extrapolation. 
(e) 
Sparse Grid method
D01ESF implements a sparse grid quadrature scheme for the integration of a vector of multidimensional integrals over the unit hypercube,
The routine uses a vectorized interface, which returns a set of points at which the integrands must be evaluated in a sparse storage format for efficiency.
Other domains can be readily integrated over by using an appropriate mapping inside the provided subroutine for evaluating the integrands. It is suitable for $d$ up to $\mathit{O}\left(100\right)$, although no upper bound on the number of dimensions is enforced. It will also evaluate onedimensional integrals, although in this case the sparse grid used is in fact the full grid.
The routine uses optional parameters, set and queried using the routines D01ZKF and D01ZLF respectively. Amongst other options, these allow the parallelization of the routine to be controlled. 
(f) 
Automatic routines
(D01FCF and D01GBF)
Both routines are for integrals of the form
D01GBF
is an adaptive Monte–Carlo routine. This routine is usually slow and not recommended for highaccuracy work. It is a robust routine that can often be used for lowaccuracy results with highly irregular integrands or when $n$ is large.
D01FCF
is an adaptive deterministic routine. Convergence is fast for well behaved integrands. Highly accurate results can often be obtained for $n$ between $2$ and $5$, using significantly fewer integrand evaluations than would be required by
D01GBF.
The routine will usually work when the integrand is mildly singular and for $n\le 10$ should be used before
D01GBF.
If it is known in advance that the integrand is highly irregular, it is best to compare results from at least two different routines.
There are many problems for which one or both of the routines will require large amounts of computing time to obtain even moderately accurate results. The amount of computing time is controlled by the number of integrand evaluations you have allowed, and you should set this argument carefully, with reference to the time available and the accuracy desired.
D01EAF extends the technique of D01FCF to integrate adaptively more than one integrand, that is to calculate the set of integrals
for a set of similar integrands ${f}_{1},{f}_{2},\dots ,{f}_{m}$ where ${f}_{i}={f}_{i}\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$. 
4 Decision Trees
Tree 1: Onedimensional integrals over a finite interval
Is the functional form of the integrand known? 

Is indefinite integration required? 

D01ARF 
yes  yes 
 no    no  

Do you require reverse communication? 

D01RAF 
 yes 
  no  

Are you concerned with efficiency for simple integrals? 

Is the integrand smooth (polynomiallike) apart from weight function ${\leftx\left(a+b\right)/2\right}^{c}$ or ${\left(bx\right)}^{c}{\left(xa\right)}^{d}$? 

D01ARF, D01UAF, D01TBF or
D01BCF and D01FBF,
or D01GCF 
 yes  yes 
  no    no  


Is the integrand reasonably smooth and the required accuracy not too great? 

D01ARF, D01BDF, D01ESF or D01UAF 
  yes 
   no  


Are multiple integrands to be integrated simultaneously? 

D01ESF or D01RAF 
  yes 
   no  


Has the integrand discontinuities, sharp peaks or singularities at known points other than the end points? 

Split the range and begin again; or use D01ALF or D01RGF 
  yes 
   no  


Is the integrand free of singularities, sharp peaks and violent oscillations apart from weight function ${\left(bx\right)}^{\alpha}{\left(xa\right)}^{\beta}\phantom{\rule{0ex}{0ex}}{\left(\mathrm{log}\left(bx\right)\right)}^{k}{\left(\mathrm{log}\left(xa\right)\right)}^{l}$?


D01APF 
  yes 
   no  


Is the integrand free of singularities, sharp peaks and violent oscillations apart from weight function ${\left(xc\right)}^{1}$? 

D01AQF 
  yes 
   no  


Is the integrand free of violent oscillations apart from weight function $\mathrm{cos}\left(\omega x\right)$ or $\mathrm{sin}\left(\omega x\right)$? 

D01ANF 
  yes 
   no  


Is the integrand free of singularities? 

D01AJF, D01AKF, D01AUF or D01ESF 
  yes 
   no  


Is the integrand free of discontinuities and of singularities except possibly at the end points? 

D01AHF 
  yes 
   no  


D01AJF, D01ATF, D01RAF or D01RGF 
 
 

D01AHF, D01AJF, D01ATF, D01RAF or D01RGF 


D01GAF 

Note: D01ATF,
D01AUF,
D01RAF and
D01RGF are likely to be more efficient due to their vectorized interfaces than
D01AJF and
D01AKF, which use a more conventional userinterface, consistent with other routines in the chapter.
Tree 2: Onedimensional integrals over a semiinfinite or infinite interval
Is the functional form of the integrand known? 

Are you concerned with efficiency for simple integrands? 

Is the integrand smooth (polynomiallike) with no exceptions? 

D01UAF, D01BDF, D01ARF or D01ESF with transformation. See Section 3.4 (b)(ii). 
yes  yes  yes 
 no    no    no  


Is the integrand of the form $\underset{0}{\overset{\infty}{\int}}}\mathrm{exp}\left({x}^{2}\right)\text{}g\left(x\right)dx$? 

D01UBF 
  yes 
   no  


Is the integrand smooth (polynomiallike) apart from weight function ${e}^{\beta \left(x\right)}$ (semiinfinite range) or ${e}^{{\beta \left(xa\right)}^{2}}$ (infinite range) or is the integrand polynomiallike in $\frac{1}{x+b}$? (semiinfinite range)? 

D01UAF, or
D01BCF and D01FBF, or,
D01TBF and D01FBF, or
D01TDF and D01FBF
(D01TDF may require use of D01TEF) 
  yes 
   no  


Has integrand discontinuities, sharp peaks or singularities at known points other than a finite limit? 

Split range; begin again using finite or infinite range trees 
  yes 
   no  


Does the integrand oscillate over the entire range? 

Does the integrand decay rapidly towards an infinite limit? 

Use D01AMF; or set cutoff and use finite range tree 
  yes  yes 
   no    no  



Is the integrand free of violent oscillations apart from weight function $\mathrm{cos}\left(\omega x\right)$ or $\mathrm{sin}\left(\omega x\right)$ (semiinfinite range)? 

D01ASF 
   yes 
    no  



Use finiterange integration between the zeros and extrapolate (see C06BAF) 
  
  


D01AMF 
 
 

D01AMF 


D01GAF (integrates over the range of the points supplied) 

Tree 3: Multidimensional integrals
Is dimension $\text{}=2$ and product region? 

D01DAF 
yes 
 no  
Is dimension $\text{}\le 4$ 

Is region an $n$sphere? 

D01FBF with user transformation or D01JAF 
yes  yes 
 no    no  

Is region a Simplex? 

D01FBF with user transformation or D01PAF 
 yes 
  no  

Is the integrand smooth (polynomiallike) in each dimension apart from weight function? 

D01TBF or D01BCF with D01FBF 
 yes 
  no  

Is integrand free of extremely bad behaviour? 

D01ESF, D01FCF, D01FDF or D01GCF

 yes 
  no  

Is bad behaviour on the boundary? 

D01FCF or D01FDF 
 yes 
  no  

Compare results from at least two of D01FCF, D01FDF, D01GBF and D01GCF, D01ESF and onedimensional recursive application 


Is region an $n$sphere? 

D01FDF 
yes 
 no  
Is region a Simplex? 

D01PAF 
yes 
 no  
Is high accuracy required? 

D01FDF with argument tuning 
yes 
 no  
Is dimension high? 

D01FDF, D01GCF or D01GDF, D01ESF

yes 
 no  
D01FCF 

Note: in the case where there are many integrals to be evaluated
D01EAF should be preferred to
D01FCF.
D01GDF is likely to be more efficient than
D01GCF, which uses a more conventional userinterface, consistent with other routines in the chapter.
5 Functionality Index
when number of points is a product of 2 primes   D01GZF 
when number of points is prime   D01GYF 
Multidimensional quadrature,   
over a finite twodimensional region   D01DAF 
over a general product region,   
Korobov–Conroy numbertheoretic method   D01GCF 
Sag–Szekeres method (also over nsphere)   D01FDF 
variant of D01GCF especially efficient on vector machines   D01GDF 
Gaussian quadrature ruleevaluation   D01FBF 
sparse grid method (with user transformation),   
muliple integrands, vectorized interface   D01ESF 
over an nsphere (n ≤ 4),   
allowing for badly behaved integrands   D01JAF 
Onedimensional quadrature,   
adaptive integration of a function over a finite interval,   
suitable for badly behaved integrals,   
strategy due to Patterson,   
suitable for wellbehaved integrands, except possibly at endpoints   D01AHF 
strategy due to Piessens and de Doncker,   
allowing for singularities at userspecified breakpoints   D01ALF 
suitable for badly behaved integrands,   
single abscissa interface   D01AJF 
suitable for highly oscillatory integrals,   
single abscissa interface   D01AKF 
weight function 1 / (x − c) Cauchy principal value (Hilbert transform)   D01AQF 
weight function cos(ωx) or sin(ωx)   D01ANF 
weight function with endpoint singularities of algebraicologarithmic type   D01APF 
adaptive integration of a function over an infinite interval or semiinfinite interval,   
weight function cos(ωx) or sin(ωx)   D01ASF 
integration of a function defined by data values only,   
nonadaptive integration over a finite, semiinfinite or infinite interval,   
using precomputed weights and abscissae   
specific integral with weight exp( − x^{2}) over semiinfinite interval   D01UBF 
nonadaptive integration over a finite interval   D01BDF 
nonadaptive integration over a finite interval,   
with provision for indefinite integrals also   D01ARF 
adaptive integration over a finite interval,   
efficient on vector machines   D01RAF 
general option setting and initialization   D01ZKF 
Weights and abscissae for Gaussian quadrature rules,   
method of Golub and Welsch,   
calculating the weights and abscissae   D01TDF 
generate recursive coefficients   D01TEF 
more general choice of rule,   
calculating the weights and abscissae   D01BCF 
restricted choice of rule,   
using precomputed weights and abscissae   D01TBF 
6 Auxiliary Routines Associated with Library Routine Arguments
D01FDV  nagf_quad_md_sphere_dummy_region See the description of the argument
REGION in D01FDF. 
D01RBM  nagf_quad_d01rb_dummy See the description of the argument MONIT in D01RBF. 
7 Routines Withdrawn or Scheduled for Withdrawal
The following lists all those routines that have been withdrawn since Mark 19 of the Library or are scheduled for withdrawal at one of the next two marks.
8 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