NAG Library Chapter Introduction
d01 – Quadrature
1 Scope of the Chapter
This chapter provides functions for the numerical evaluation of definite integrals in one or more dimensions.
2 Background to the Problems
The functions 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) requested by you. 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 functions 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 functions 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) 
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 Functions
The following three subsections consider in turn functions for: onedimensional integrals over a finite interval, and over a semiinfinite or an infinite interval; and multidimensional integrals. Within each subsection, functions 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 functions 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 function. 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 function; a simple ruleevaluation function may provide much better overall performance. For further guidance, the article by
Lyness (1983) is recommended.
3.1 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 ( nag_1d_quad_vals (d01gac)) 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 function 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 nag_1d_cheb_fit (e02adc) and then use functions nag_1d_cheb_intg (e02ajc) and nag_1d_cheb_eval2 (e02akc) to evaluate its integral (which need not be restricted to the range of the integration points, as is the case for nag_1d_quad_vals (d01gac)). A further alternative is to fit a cubic spline to the data using nag_1d_spline_fit_knots (e02bac) and then to evaluate its integral using nag_1d_spline_intg (e02bdc). 
(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 functions to use will be
nag_1d_quad_gen_1 (d01sjc).
(i) 
Ruleevaluation functions
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 function with a suitable number of abscissae may be used.
If $f\left(x\right)$ is well behaved, apart from a weightfunction of the form
nag_quad_1d_gauss_wgen (d01tcc)
with nag_quad_md_gauss (d01fbc) may be used. 
(ii) 
Automatic adaptive functions
Firstly, several functions 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.
 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
nag_1d_quad_wt_alglog_1 (d01spc);
 if $w\left(x\right)=\frac{1}{xc}$: use
nag_1d_quad_wt_cauchy_1 (d01sqc)
(this integral is called the Hilbert transform of $g$);
 if $w\left(x\right)=\mathrm{cos}\left(\omega x\right)$ or $\mathrm{sin}\left(\omega x\right)$: use nag_1d_quad_wt_trig_1 (d01snc) (this function can also handle certain types of singularities in $g\left(x\right)$).
Secondly, there are some functions for general $f\left(x\right)$. If $f\left(x\right)$ is known to be free of singularities, though it may be oscillatory, nag_1d_quad_osc_1 (d01skc) may be used.
The most powerful of the finite interval integration functions are
nag_1d_quad_gen_1 (d01sjc) and nag_1d_quad_brkpts_1 (d01slc),
which can cope with singularities of several types.
They may be used if none of the more specific situations described above
applies and are
somewhat more reliable, particularly where the integrand has singularities other than at an end point, or has discontinuities or cusps, and are therefore recommended where the integrand is known to be badly
behaved. nag_1d_quad_brkpts_1 (d01slc) is recommended when the locations of singularities are known while nag_1d_quad_gen_1 (d01sjc) is recommended when the nature of the integrand is completely unknown.
If $f\left(x\right)$ has singularities of certain types, discontinuities (including switches in definition) or sharp peaks occurring at known points, the integral should be evaluated separately over each of the subranges or
nag_1d_quad_brkpts_1 (d01slc)
may be used. 

3.2 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, nag_1d_quad_vals (d01gac), should be used. 
(b) 
Integrand defined as a function
(i) 
Rule evaluation functions
If $f\left(x\right)$ behaves approximately like a polynomial in $x$, apart from a weight function of the form:
 ${e}^{\beta x},\beta >0$ (semiinfinite interval, lower limit finite); or
 ${e}^{\beta x},\beta <0$ (semiinfinite interval, upper limit finite); or
 ${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 functions may be used.
nag_1d_quad_gauss_1 (d01tac)
may be used if it is not required to examine the weights and abscissae.

(ii) 
Automatic adaptive functions
nag_1d_quad_inf_1 (d01smc)
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.
nag_1d_quad_inf_wt_trig_1 (d01ssc)
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.
 If the integrand decays rapidly towards an infinite end point, a finite cutoff may be chosen, and the finite range methods applied.
 If the only irregularities occur in the finite part (apart from a singularity at the finite limit, with which
nag_1d_quad_inf_1 (d01smc)
can cope), the range may be divided, with
nag_1d_quad_inf_1 (d01smc) used on the infinite part.
 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.3 Multidimensional Integrals
(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
nag_quad_1d_gauss_wset (d01tbc) or nag_quad_1d_gauss_wgen (d01tcc)
with nag_quad_md_gauss (d01fbc). 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 nag_quad_2d_fin (d01dac) 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.2); 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 functions proves suitable and convenient, the onedimensional functions 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 function 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). Note that, for reasons of efficiency and robustness the integration functions are not defined as recursive, and thus a different library integration function must be used for each dimension. Apart from this restriction, the following combinations are not permitted:
nag_1d_quad_gen_1 (d01sjc) and nag_1d_quad_brkpts_1 (d01slc),
nag_1d_quad_wt_trig_1 (d01snc) and nag_1d_quad_wt_alglog_1 (d01spc),
nag_1d_quad_wt_alglog_1 (d01spc) and nag_1d_quad_wt_cauchy_1 (d01sqc),
nag_1d_quad_wt_trig_1 (d01snc) and nag_1d_quad_wt_cauchy_1 (d01sqc),
nag_1d_quad_wt_trig_1 (d01snc) and nag_1d_quad_inf_wt_trig_1 (d01ssc),
nag_1d_quad_inf_1 (d01smc) and nag_1d_quad_inf_wt_trig_1 (d01ssc),
nag_1d_quad_gen_1 (d01sjc) and nag_1d_quad_osc_1 (d01skc).
Otherwise the full range of onedimensional functions are available, for finite/infinite intervals, constant/variable limits, rule evaluation/automatic strategies etc. 
(b) 
Sag–Szekeres method
nag_quad_md_sphere (d01fdc) 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. 
(c) 
A combinatorial extrapolation method
nag_quad_md_simplex (d01pac) 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. 
(d) 
Automatic functions
(nag_multid_quad_adapt_1 (d01wcc) and nag_multid_quad_monte_carlo_1 (d01xbc))
Both functions are for integrals of the form
nag_multid_quad_monte_carlo_1 (d01xbc)
is an adaptive Monte–Carlo function. This function is usually slow and not recommended for highaccuracy work. It is a robust function that can often be used for lowaccuracy results with highly irregular integrands or when $n$ is large.
nag_multid_quad_adapt_1 (d01wcc)
is an adaptive deterministic function. 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
nag_multid_quad_monte_carlo_1 (d01xbc).
The function will usually work when the integrand is mildly singular and for $n\le 10$ should be used before
nag_multid_quad_monte_carlo_1 (d01xbc).
If it is known in advance that the integrand is highly irregular, it is best to compare results from at least two different functions.
There are many problems for which one or both of the functions 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 allowed by you, and you should set this argument carefully, with reference to the time available and the accuracy desired. 
4 Decision Trees
Tree 1: Onedimensional integrals over a finite interval
Is the functional form of the integrand known? 
_ yes 
Are you concerned with efficiency for simple integrals? 
_ yes 
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}$? 
_ yes 
d01tac, d01tbc or d01tcc and d01fbc or d01gdc 
 

 

no  

 
  

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

 

no  

 
  

Has the integrand discontinuities, sharp peaks or singularities at known points other than the end points? 
_ yes 
Split the range and begin again; or use d01slc 
 

 

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}$?

_ yes 
d01spc 
 

 

no  

 
  

Is the integrand free of singularities, sharp peaks and violent oscillations apart from weight function $\frac{1}{xc}$? 
_ yes 
d01sqc 
 

 

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)$? 
_ yes 
d01snc 
 

 

no  

 
  

Is the integrand free of singularities? 
_ yes 
d01skc 
 

 

no  

 
  

d01sjc 
 

no  

 

d01sjc 
no  

d01gac 
Tree 2: Onedimensional integrals over a semiinfinite or infinite interval
Is the functional form of the integrand known? 
_ yes 
Are you concerned with efficiency for simple integrands? 
_ yes 
Is the integrand smooth (polynomiallike) with no exceptions? 
_ yes 
d01bdc 
 

 

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)? 
_ yes 
d01tac, d01tbc and d01fbc or d01tcc and d01fbc 
 

 

no  

 
  

Has integrand discontinuities, sharp peaks or singularities at known points other than a finite limit? 
_ yes 
Split range; begin again using finite or infinite range trees 
 

 

no  

 
  

Does the integrand oscillate over the entire range? 
_ yes 
Does the integrand decay rapidly towards an infinite limit? 
_ yes 
Use d01smc; or set cutoff and use finite range tree 
 

 

 

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)? 
_ yes 
d01ssc 
 

 

 

no  

 
  
  

Use finiterange integration between the zeros and extrapolate. 
 

 

no  

 
  

d01smc 
 

no  

 

d01smc 
no  

d01gac (integrates over the range of the points supplied) 
Tree 3: Multidimensional integrals
Is dimension $\text{}=2$ and product region? 
_ yes 
d01dac 
no  

Is dimension $\text{}\le 4$ 
_ yes 
Is region an $n$sphere? 
_ yes 
d01fbc with user transformation. 
 

no  

 

Is region a Simplex? 
_ yes 
d01fbc with user transformation or d01pac 
 

no  

 

Is the integrand smooth (polynomiallike) in each dimension apart from weight function? 
_ yes 
d01tbc and d01fbc or d01tcc and d01fbc 
 

no  

 

Is integrand free of extremely bad behaviour? 
_ yes 
d01fdc, d01gdc or d01wcc 
 

no  

 

Is bad behaviour on the boundary? 
_ yes 
d01fdc or d01wcc 
 

no  

 

Compare results from at least two of d01fdc, d01gdc, d01wcc and d01xbc and onedimensional recursive application 
no  

Is region an $n$sphere? 
_ yes 
d01fdc 
no  

Is region a Simplex? 
_ yes 
d01pac 
no  

Is high accuracy required? 
_ yes 
d01fdc with argument tuning 
no  

Is dimension high? 
_ yes 
d01fdc or d01gdc 
no  

d01wcc 
5 Functionality Index
Multidimensional quadrature:   
over a general product region:   
Onedimensional quadrature:   
adaptive integration of a function over a finite interval:   
strategy due to Piessens and de Doncker,   
adaptive integration of a function over an infinite interval or semiinfinite interval:   
integration of a function defined by data values only,   
Weights and abscissae for Gaussian quadrature rules:   
more general choice of rule,   
restricted choice of rule,   
6 Functions Withdrawn or Scheduled for Withdrawal
7 References
Davis P J and Rabinowitz P (1975) Methods of Numerical Integration Academic Press
Gladwell I (1986) Vectorisation of onedimensional quadrature codes Numerical Integration: Recent Developments, and Applications (eds P Keast and G Fairweather) 231–238 D Reidel Publishing Company, Holland
Lyness J N (1983) When not to use an automatic quadrature routine SIAM Rev. 25 63–87
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