d01gc:: Quadrature (NAG Toolbox)

nag_quad_md_numth (d01gc) calculates an approximation to a definite integral in up to

20

dimensions, using the Korobov–Conroy number theoretic method.

nag_quad_md_numth (d01gc) calculates an approximation to the integral

I = \int_{c_{1}}^{d_{1}} d x_{1}, \dots, \int_{c_{n}}^{d_{n}} d x_{n} f (x_{1}, x_{2}, \dots, x_{n})

(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

{[0, 1]}^{n}

by the change of variables

x_{i} = c_{i} + (d_{i} - c_{i}) y_{i}, i = 1, 2, \dots, n .

The method then uses as its basis the number theoretic formula for the

n

-cube,

{[0, 1]}^{n}

\int_{0}^{1} d x_{1} \dots \int_{0}^{1} d x_{n} g (x_{1}, x_{2}, \dots, x_{n}) = \frac{1}{p} \sum_{k = 1}^{p} g (\{k \frac{a_{1}}{p}\}, \dots, \{k \frac{a_{n}}{p}\}) - E

(2)

where

\{x\}

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 (x_{1}, x_{2}, \dots, x_{n})

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 \leq C K p^{- α} \ln^{α β} p

(3)

where

α

and

C

are real numbers depending on the convergence rate of the Fourier series,

β

is a constant depending on

n

, and

K

is a constant depending on

α

and

n

. There are a number of procedures for calculating these optimal coefficients. Korobov imposes the constraint that

a_{1} = 1 and a_{i} = a^{i - 1} (\mod p)

(4)

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

x_{i} = y_{i}^{2} (3 - 2 y_{i}), i = 1, 2, \dots, 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),

\{k \frac{a_{i}}{p}\}

is replaced by

\{α_{i} + k \frac{a_{i}}{p}\}

, for

i = 1, 2, \dots, n

, where each

α_{i}

, is uniformly distributed over

[0, 1]

. Computing the integral for each of a sequence of random vectors

α

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 nag_quad_md_numth_coeff_prime (d01gy) and nag_quad_md_numth_coeff_2prime (d01gz) 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

Korobov N M (1963) Number Theoretic Methods in Approximate Analysis Fizmatgiz, Moscow

f must return the value of the integrand

f

at a given point.

n

, the number of dimensions of the integral.

The coordinates of the point at which the integrand

f

must be evaluated.

The value of

f (x)

evaluated at x.

region must evaluate the limits of integration in any dimension.

n

, the number of dimensions of the integral.

x (1), \dots, x (j - 1)

contain the current values of the first

(j - 1)

variables, which may be used if necessary in calculating

c_{j}

and

d_{j}

The index

j

for which the limits of the range of integration are required.

The lower limit

c_{j}

of the range of

x_{j}

The upper limit

d_{j}

of the range of

x_{j}

The Korobov rule to be used. There are two alternatives depending on the value of npts.

(i)	$1 \leq npts \leq 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)	$npts > 6$ . npts is the number of actual points to be used with corresponding optimal coefficients supplied in the array vk.

Constraint:

npts \geq 1

npts > 6

, vk must contain the

n

optimal coefficients (which may be calculated using nag_quad_md_numth_coeff_prime (d01gy) or nag_quad_md_numth_coeff_2prime (d01gz)).

npts \leq 6

, vk need not be set.

The number of random samples to be generated in the error estimation (generally a small value, say

3

5

, is sufficient). The total number of integrand evaluations will be

nrand \times npts

Constraint:

nrand \geq 1

Default: the dimension of the array vk.

n

, the number of dimensions of the integral.

Constraint:

1 \leq ndim \leq 20

Indicates whether the periodising transformation is to be used.

$itrans ='0'$: The transformation is to be used.
$itrans \neq'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).

npts > 6

, vk is unchanged.

npts \leq 6

, vk contains the

n

optimal coefficients used by the preset rule.

The approximation to the integral

I

The standard error as computed from nrand sample values. If

nrand = 1

, then err contains zero.

ifail = 0

unless the function detects an error (see Error Indicators and Warnings).

Errors or warnings detected by the function:

An unexpected error has been triggered by this routine. Please contact NAG.

Your licence key may have expired or may not have been installed correctly.

Dynamic memory allocation failed.

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

The time taken by nag_quad_md_numth (d01gc) will be approximately proportional to

nrand \times 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 nag_quad_md_numth (d01gc) by calls to nag_rand_dist_uniform01 (g05sa). Separate runs will produce identical answers.

This example calculates the integral

\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \cos (0.5 + 2 (x_{1} + x_{2} + x_{3} + x_{4}) - 4) d x_{1} d x_{2} d x_{3} d x_{4} .

On entry,	$ndim < 1$ ,
or	$ndim > 20$ .

NAG Toolbox: nag_quad_md_numth (d01gc)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example