d01gd calculates an approximation to a definite integral in up to $20$ dimensions, using the Korobov–Conroy number theoretic method. This method is designed to be particularly efficient on vector processors.

# Syntax

C# |
---|

public static void d01gd( int ndim, D01..::..D01GD_VECFUN vecfun, D01..::..D01GD_VECREG vecreg, int npts, double[] vk, int nrand, int itrans, out double res, out double err, out int ifail ) |

Visual Basic |
---|

Public Shared Sub d01gd ( _ ndim As Integer, _ vecfun As D01..::..D01GD_VECFUN, _ vecreg As D01..::..D01GD_VECREG, _ npts As Integer, _ vk As Double(), _ nrand As Integer, _ itrans As Integer, _ <OutAttribute> ByRef res As Double, _ <OutAttribute> ByRef err As Double, _ <OutAttribute> ByRef ifail As Integer _ ) |

Visual C++ |
---|

public: static void d01gd( int ndim, D01..::..D01GD_VECFUN^ vecfun, D01..::..D01GD_VECREG^ vecreg, int npts, array<double>^ vk, int nrand, int itrans, [OutAttribute] double% res, [OutAttribute] double% err, [OutAttribute] int% ifail ) |

F# |
---|

static member d01gd : ndim : int * vecfun : D01..::..D01GD_VECFUN * vecreg : D01..::..D01GD_VECREG * npts : int * vk : float[] * nrand : int * itrans : int * res : float byref * err : float byref * ifail : int byref -> unit |

#### Parameters

- ndim
- Type: System..::..Int32
*On entry*: $n$, the number of dimensions of the integral.*Constraint*: $1\le {\mathbf{ndim}}\le 20$.

- vecfun
- Type: NagLibrary..::..D01..::..D01GD_VECFUNvecfun must evaluate the integrand at a specified set of points.
A delegate of type D01GD_VECFUN.

- vecreg
- Type: NagLibrary..::..D01..::..D01GD_VECREGvecreg must evaluate the limits of integration in any dimension for a set of points.
A delegate of type D01GD_VECREG.

- npts
- Type: System..::..Int32
*On entry*: the Korobov rule to be used. There are two alternatives depending on the value of npts.(i) $1\le {\mathbf{npts}}\le 6$. In this case one of six preset rules is chosen using $2129$, $5003$, $10007$, $20011$, $40009$ or $80021$ points depending on the respective value of npts being $1$, $2$, $3$, $4$, $5$ or $6$.(ii) ${\mathbf{npts}}>6$. *Constraint*: ${\mathbf{npts}}\ge 1$.

- vk
- Type: array<System..::..Double>[]()[][]An array of size [ndim]

- nrand
- Type: System..::..Int32
*On entry*: the number of random samples to be generated (generally a small value, say $3$ to $5$, is sufficient). The estimate, res, of the value of the integral returned by the method is then the average of nrand calculations with different random origin shifts. If ${\mathbf{npts}}>6$, the total number of integrand evaluations will be ${\mathbf{nrand}}\times {\mathbf{npts}}$. If $1\le {\mathbf{npts}}\le 6$, then the number of integrand evaluations will be ${\mathbf{nrand}}\times p$, where $p$ is the number of points corresponding to the six preset rules. For reasons of efficiency, these values are calculated a number at a time in vecfun.*Constraint*: ${\mathbf{nrand}}\ge 1$.

- itrans
- Type: System..::..Int32
*On entry*: indicates whether the periodising transformation is to be used.- ${\mathbf{itrans}}=0$
- The transformation is to be used.
- ${\mathbf{itrans}}=1$
- The transformation is to be suppressed (to cover cases where the integrand may already be periodic or where you want to specify a particular transformation in the definition of vecfun).

*Suggested value*: ${\mathbf{itrans}}=0$.

- res
- Type: System..::..Double%
*On exit*: the approximation to the integral $I$.

- err
- Type: System..::..Double%

- ifail
- Type: System..::..Int32%
*On exit*: ${\mathbf{ifail}}={0}$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

# Description

d01gd 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 ${\left[0,1\right]}^{n}$ by the change of variables

The method then uses as its basis the number theoretic formula for the $n$-cube, ${\left[0,1\right]}^{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\left({x}_{1},\dots ,{x}_{n}\right)$ 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 parameter, $a$, to satisfy the optimal conditions.

$$I=\underset{{c}_{1}}{\overset{{d}_{1}}{\int}}\cdots \underset{{c}_{n}}{\overset{{d}_{n}}{\int}}f\left({x}_{1},\dots ,{x}_{n}\right)d{x}_{n}\dots d{x}_{1}$$ | (1) |

$${x}_{i}={c}_{i}+\left({d}_{i}-{c}_{i}\right){y}_{i}\text{, \hspace{1em}}i=1,2,\dots ,n\text{.}$$ |

$$\underset{0}{\overset{1}{\int}}\cdots \underset{0}{\overset{1}{\int}}g\left({x}_{1},\dots ,{x}_{n}\right)d{x}_{n}\cdots d{x}_{1}=\frac{1}{p}\sum _{k=1}^{p}g\left(\left\{k\frac{{a}_{1}}{p}\right\},\dots ,\left\{k\frac{{a}_{n}}{p}\right\}\right)-E$$ | (2) |

$$E\le CK{p}^{-\alpha}{\mathrm{ln}}^{\alpha \beta}\u200ap$$ | (3) |

$${a}_{1}=1\text{\hspace{1em} and \hspace{1em}}{a}_{i}={a}^{i-1}\left(\mathrm{mod}\text{}p\right)$$ | (4) |

In this method 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.

$${x}_{i}={y}_{i}^{2}\left(3-2{y}_{i}\right)\text{, \hspace{1em}}i=1,2,\dots ,n\text{.}$$ |

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 $\left\{{\alpha}_{i}+k\frac{{a}_{i}}{p}\right\}$, for $i=1,2,\dots ,n$, where each ${\alpha}_{i}$, is uniformly distributed over $\left[0,1\right]$. Computing the integral for each of a sequence of random vectors $\alpha $ allows a ‘standard error’ to be estimated.

This method provides built-in sets of optimal coefficients, corresponding to six different values of $p$. Alternatively, the optimal coefficients may be supplied by you. Methods d01gy and d01gz compute the optimal coefficients for the cases where $p$ is a prime number or $p$ is a product of two primes, respectively.

# References

Conroy H (1967) Molecular Shroedinger equation VIII. A new method for evaluting multi-dimensional integrals

*J. Chem. Phys.***47**5307–5318Cranley R and Patterson T N L (1976) Randomisation of number theoretic methods for mulitple integration

*SIAM J. Numer. Anal.***13**904–914Korobov N M (1957) The approximate calculation of multiple integrals using number theoretic methods

*Dokl. Acad. Nauk SSSR***115**1062–1065Korobov N M (1963)

*Number Theoretic Methods in Approximate Analysis*Fizmatgiz, Moscow# Error Indicators and Warnings

Errors or warnings detected by the method:

- ${\mathbf{ifail}}=1$
On entry, ${\mathbf{ndim}}<1$, or ${\mathbf{ndim}}>20$.

- ${\mathbf{ifail}}=2$
On entry, ${\mathbf{npts}}<1$.

- ${\mathbf{ifail}}=3$
On entry, ${\mathbf{nrand}}<1$.

- ${\mathbf{ifail}}=-9000$
- An error occured, see message report.
- ${\mathbf{ifail}}=-6000$
- Invalid Parameters $\u2329\mathit{\text{value}}\u232a$
- ${\mathbf{ifail}}=-8000$
- Negative dimension for array $\u2329\mathit{\text{value}}\u232a$
- ${\mathbf{ifail}}=-6000$
- Invalid Parameters $\u2329\mathit{\text{value}}\u232a$

# Accuracy

If ${\mathbf{nrand}}>1$, an estimate of the absolute standard error is given by the value, on exit, of err.

# Parallelism and Performance

None.

# Further Comments

d01gd performs the same computation as (D01GCF not in this release). However, the interface has been modified so that it can perform more efficiently on machines with vector processing capabilities. In particular,
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 methods, which must be supplied by you, may account for a significant part of the total computation time.
For this reason it is vital that you consider the possibilities for vectorization in the code supplied for these two methods.

# Example

This example calculates the integral

$$\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1}{\int}}\mathrm{cos}\left(0.5+2\left({x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}\right)-4\right)d{x}_{1}d{x}_{2}d{x}_{3}d{x}_{4}\text{.}$$ |