D01 Chapter Contents
D01 Chapter Introduction
NAG Library Manual

# NAG Library Routine DocumentD01GCF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

## 1  Purpose

D01GCF calculates an approximation to a definite integral in up to $20$ dimensions, using the Korobov–Conroy number theoretic method.

## 2  Specification

 SUBROUTINE D01GCF ( NDIM, F, REGION, NPTS, VK, NRAND, ITRANS, RES, ERR, IFAIL)
 INTEGER NDIM, NPTS, NRAND, ITRANS, IFAIL REAL (KIND=nag_wp) F, VK(NDIM), RES, ERR EXTERNAL F, REGION

## 3  Description

D01GCF calculates an approximation to the integral
 $I= ∫ c1 d1 dx1 ,…, ∫ cn dn dxn f x1,x2,…,xn$ (1)
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
 $xi = ci + di - ci yi , i= 1 , 2 ,…, n .$
The method then uses as its basis the number theoretic formula for the $n$-cube, ${\left[0,1\right]}^{n}$:
 $∫01 dx1 ⋯ ∫01 dxn g x1,x2,…,xn = 1p ∑k=1p g k a1p ,…, k anp - E$ (2)
where $\left\{x\right\}$ denotes the fractional part of $x$, ${a}_{1},{a}_{2},\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},{a}_{2},\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},{x}_{2},\dots ,{x}_{n}\right)$ which is assumed to have some degree of periodicity. Depending on the choice of ${a}_{1},{a}_{2},\dots ,{a}_{n}$ the contributions from certain groups of Fourier coefficients are eliminated from the error, $E$. Korobov shows that ${a}_{1},{a}_{2},\dots ,{a}_{n}$ can be chosen so that the error satisfies
 $E≤CK p-α ln αβ ⁡p$ (3)
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
 $a1 = 1 and ai = ai-1 mod p$ (4)
and gives a procedure for calculating the parameter, $a$, to satisfy the optimal conditions.
In this routine the periodisation is achieved by the simple transformation
 $xi = yi2 3-2yi , i= 1 , 2 ,…, n .$
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 $\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 routine provides built-in sets of optimal coefficients, corresponding to six different values of $p$. Alternatively, the optimal coefficients may be supplied by you. Routines D01GYF and D01GZF compute the optimal coefficients for the cases where $p$ is a prime number or $p$ is a product of two primes, respectively.

## 4  References

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
Korobov N M (1963) Number Theoretic Methods in Approximate Analysis Fizmatgiz, Moscow

## 5  Parameters

1:     NDIM – INTEGERInput
On entry: $n$, the number of dimensions of the integral.
Constraint: $1\le {\mathbf{NDIM}}\le 20$.
2:     F – REAL (KIND=nag_wp) FUNCTION, supplied by the user.External Procedure
F must return the value of the integrand $f$ at a given point.
The specification of F is:
 FUNCTION F ( NDIM, X)
 REAL (KIND=nag_wp) F
 INTEGER NDIM REAL (KIND=nag_wp) X(NDIM)
1:     NDIM – INTEGERInput
On entry: $n$, the number of dimensions of the integral.
2:     X(NDIM) – REAL (KIND=nag_wp) arrayInput
On entry: the coordinates of the point at which the integrand $f$ must be evaluated.
F must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D01GCF is called. Parameters denoted as Input must not be changed by this procedure.
3:     REGION – SUBROUTINE, supplied by the user.External Procedure
REGION must evaluate the limits of integration in any dimension.
The specification of REGION is:
 SUBROUTINE REGION ( NDIM, X, J, C, D)
 INTEGER NDIM, J REAL (KIND=nag_wp) X(NDIM), C, D
1:     NDIM – INTEGERInput
On entry: $n$, the number of dimensions of the integral.
2:     X(NDIM) – REAL (KIND=nag_wp) arrayInput
On entry: ${\mathbf{X}}\left(1\right),\dots ,{\mathbf{X}}\left(j-1\right)$ contain the current values of the first $\left(j-1\right)$ variables, which may be used if necessary in calculating ${c}_{j}$ and ${d}_{j}$.
3:     J – INTEGERInput
On entry: the index $j$ for which the limits of the range of integration are required.
4:     C – REAL (KIND=nag_wp)Output
On exit: the lower limit ${c}_{j}$ of the range of ${x}_{j}$.
5:     D – REAL (KIND=nag_wp)Output
On exit: the upper limit ${d}_{j}$ of the range of ${x}_{j}$.
REGION must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D01GCF is called. Parameters denoted as Input must not be changed by this procedure.
4:     NPTS – INTEGERInput
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$. NPTS is the number of actual points to be used with corresponding optimal coefficients supplied in the array VK.
Constraint: ${\mathbf{NPTS}}\ge 1$.
5:     VK(NDIM) – REAL (KIND=nag_wp) arrayInput/Output
On entry: if ${\mathbf{NPTS}}>6$, VK must contain the $n$ optimal coefficients (which may be calculated using D01GYF or D01GZF).
If ${\mathbf{NPTS}}\le 6$, VK need not be set.
On exit: if ${\mathbf{NPTS}}>6$, VK is unchanged.
If ${\mathbf{NPTS}}\le 6$, VK contains the $n$ optimal coefficients used by the preset rule.
6:     NRAND – INTEGERInput
On entry: the number of random samples to be generated in the error estimation (generally a small value, say $3$ to $5$, is sufficient). The total number of integrand evaluations will be ${\mathbf{NRAND}}×{\mathbf{NPTS}}$.
Constraint: ${\mathbf{NRAND}}\ge 1$.
7:     ITRANS – INTEGERInput
On entry: indicates whether the periodising transformation is to be used.
${\mathbf{ITRANS}}=\text{'0'}$
The transformation is to be used.
${\mathbf{ITRANS}}\ne \text{'0'}$
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 F).
Suggested value: ${\mathbf{ITRANS}}=\text{'0'}$.
8:     RES – REAL (KIND=nag_wp)Output
On exit: the approximation to the integral $I$.
9:     ERR – REAL (KIND=nag_wp)Output
On exit: the standard error as computed from NRAND sample values. If ${\mathbf{NRAND}}=1$, then ERR contains zero.
10:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
On exit: ${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6  Error Indicators and Warnings

If on entry ${\mathbf{IFAIL}}={\mathbf{0}}$ or $-{\mathbf{1}}$, explanatory error messages are output on the current error message unit (as defined by X04AAF).
Errors or warnings detected by the routine:
${\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$.

## 7  Accuracy

An estimate of the absolute standard error is given by the value, on exit, of ERR.

The time taken by D01GCF will be approximately proportional to ${\mathbf{NRAND}}×p$, where $p$ is the number of points used.
The exact values of RES and ERR on return will depend (within statistical limits) on the sequence of random numbers generated within D01GCF by calls to G05SAF. Separate runs will produce identical answers.

## 9  Example

This example calculates the integral
 $∫01 ∫01 ∫01 ∫01 cos 0.5+2 x1 + x2 + x3 + x4 - 4 dx1 dx2 dx3 dx4 .$

### 9.1  Program Text

Program Text (d01gcfe.f90)

None.

### 9.3  Program Results

Program Results (d01gcfe.r)