d01gdc calculates an approximation to the integral
using the Korobov–Conroy number theoretic method (see
Korobov (1957),
Korobov (1963) and
Conroy (1967)). The region of integration defined in
(1) is such that generally
${c}_{i}$ and
${d}_{i}$ may be functions of
${x}_{1},{x}_{2},\dots ,{x}_{i-1}$, for
$i=2,3,\dots ,n$, with
${c}_{1}$ and
${d}_{1}$ constants. The integral is first of all transformed to an integral over the
$n$-cube
${[0,1]}^{n}$ by the change of variables
The method then uses as its basis the number theoretic formula for the
$n$-cube,
${[0,1]}^{n}$:
where
$\left\{x\right\}$ denotes the fractional part of
$x$,
${a}_{1},\dots ,{a}_{n}$ are the so-called optimal coefficients,
$E$ is the error, and
$p$ is a prime integer. (It is strictly only necessary that
$p$ be relatively prime to all
${a}_{1},\dots ,{a}_{n}$ and is in fact chosen to be even for some cases in
Conroy (1967).) The method makes use of properties of the Fourier expansion of
$g({x}_{1},\dots ,{x}_{n})$ which is assumed to have some degree of periodicity. Depending on the choice of
${a}_{1},\dots ,{a}_{n}$ the contributions from certain groups of Fourier coefficients are eliminated from the error,
$E$. Korobov shows that
${a}_{1},\dots ,{a}_{n}$ can be chosen so that the error satisfies
where
$\alpha $ and
$C$ are real numbers depending on the convergence rate of the Fourier series,
$\beta $ is a constant depending on
$n$, and
$K$ is a constant depending on
$\alpha $ and
$n$. There are a number of procedures for calculating these optimal coefficients. Korobov imposes the constraint that
and gives a procedure for calculating the argument,
$a$, to satisfy the optimal conditions.
In this function the periodisation is achieved by the simple transformation
More sophisticated periodisation procedures are available but in practice the degree of periodisation does not appear to be a critical requirement of the method.
An easily calculable error estimate is not available apart from repetition with an increasing sequence of values of
$p$ which can yield erratic results. The difficulties have been studied by
Cranley and Patterson (1976) who have proposed a Monte Carlo error estimate arising from converting
(2) into a stochastic integration rule by the inclusion of a random origin shift which leaves the form of the error
(3) unchanged; i.e., in the formula
(2),
$\left\{k\frac{{a}_{i}}{p}\right\}$ is replaced by
$\{{\alpha}_{i}+k\frac{{a}_{i}}{p}\}$, for
$i=1,2,\dots ,n$, where each
${\alpha}_{i}$, is uniformly distributed over
$[0,1]$. Computing the integral for each of a sequence of random vectors
$\alpha $ allows a ‘standard error’ to be estimated.
This function provides built-in sets of optimal coefficients, corresponding to six different values of
$p$. Alternatively, the optimal coefficients may be supplied by you. Functions
d01gyc and
d01gzc compute the optimal coefficients for the cases where
$p$ is a prime number or
$p$ is a product of two primes, respectively.
Conroy H (1967) Molecular Shroedinger equation VIII. A new method for evaluting multi-dimensional integrals J. Chem. Phys. 47 5307–5318
Cranley R and Patterson T N L (1976) Randomisation of number theoretic methods for mulitple integration SIAM J. Numer. Anal. 13 904–914
Korobov N M (1957) The approximate calculation of multiple integrals using number theoretic methods Dokl. Acad. Nauk SSSR 115 1062–1065
If
${\mathbf{nrand}}>1$, an estimate of the absolute standard error is given by the value, on exit, of
err.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
vecfun and
vecreg must calculate the integrand and limits of integration at a
set of points. For some problems the amount of time spent in these two functions, which must be supplied by you, may account for a significant part of the total computation time.
The time taken will be approximately proportional to
${\mathbf{nrand}}\times p$, where
$p$ is the number of points used, but may depend significantly on the efficiency of the code provided by you in
vecfun and
vecreg.
The exact values of
res and
err on return will depend (within statistical limits) on the sequence of random numbers generated within
d01gdc by calls to
g05sac. Separate runs will produce identical answers.
This example calculates the integral