Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

## Purpose

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

## Syntax

[vk, res, err, ifail] = d01gc(f, region, npts, vk, nrand, 'ndim', ndim, 'itrans', itrans)
[vk, res, err, ifail] = nag_quad_md_numth(f, region, npts, vk, nrand, 'ndim', ndim, 'itrans', itrans)

## Description

nag_quad_md_numth (d01gc) calculates an approximation to the integral
 d1 dn I = ∫ dx1, … , ∫ dxn f(x1,x2, … ,xn) c1 cn
$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 ci ${c}_{i}$ and di ${d}_{i}$ may be functions of x1 , x2 ,, xi1 ${x}_{1},{x}_{2},\dots ,{x}_{i-1}$, for i = 2 , 3 ,, n $i=2,3,\dots ,n$, with c1 ${c}_{1}$ and d1 ${d}_{1}$ constants. The integral is first of all transformed to an integral over the n$n$-cube [0,1]n ${\left[0,1\right]}^{n}$ by the change of variables
 xi = ci + (di − ci) yi ,   i = 1 , 2 , … , n . $xi = ci + ( di - ci ) yi , i= 1 , 2 ,…, n .$
The method then uses as its basis the number theoretic formula for the n$n$-cube, [0,1]n ${\left[0,1\right]}^{n}$:
 1 1 p ∫ dx1 ⋯ ∫ dxng(x1,x2, … ,xn) = 1/p ∑ g({k(a1)/p}, … ,{k(an)/p}) − E 0 0 k = 1
$∫01 dx1 ⋯ ∫01 dxn g (x1,x2,…,xn) = 1p ∑k=1p g ({ k a1p },…,{ k anp }) - E$
(2)
where {x} $\left\{x\right\}$ denotes the fractional part of x$x$, a1 , a2 ,, an ${a}_{1},{a}_{2},\dots ,{a}_{n}$ are the so-called optimal coefficients, E$E$ is the error, and p$p$ is a prime integer. (It is strictly only necessary that p$p$ be relatively prime to all a1 , a2 ,, an ${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 (x1,x2,,xn) $g\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$ which is assumed to have some degree of periodicity. Depending on the choice of a1 , a2 ,, an ${a}_{1},{a}_{2},\dots ,{a}_{n}$ the contributions from certain groups of Fourier coefficients are eliminated from the error, E$E$. Korobov shows that a1 , a2 ,, an ${a}_{1},{a}_{2},\dots ,{a}_{n}$ can be chosen so that the error satisfies
 E ≤ CK p − α lnαβ p $E≤CK p-α ln αβ ⁡p$ (3)
where α$\alpha$ and C$C$ are real numbers depending on the convergence rate of the Fourier series, β$\beta$ is a constant depending on n$n$, and K$K$ is a constant depending on α$\alpha$ and n$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) $a1 = 1 and ai = ai-1 (mod p)$ (4)
and gives a procedure for calculating the parameter, a$a$, to satisfy the optimal conditions.
In this function the periodisation is achieved by the simple transformation
 xi = yi2 (3 − 2yi) ,   i = 1 , 2 , … , n . $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$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(ai)/p} $\left\{k\frac{{a}_{i}}{p}\right\}$ is replaced by {αi + k(ai)/p} $\left\{{\alpha }_{i}+k\frac{{a}_{i}}{p}\right\}$, for i = 1 , 2 ,, n $i=1,2,\dots ,n$, where each αi ${\alpha }_{i}$, is uniformly distributed over [0,1] $\left[0,1\right]$. 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$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$p$ is a prime number or p$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–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

## Parameters

### Compulsory Input Parameters

1:     f – function handle or string containing name of m-file
f must return the value of the integrand f$f$ at a given point.
[result] = f(ndim, x)

Input Parameters

1:     ndim – int64int32nag_int scalar
n$n$, the number of dimensions of the integral.
2:     x(ndim) – double array
The coordinates of the point at which the integrand f$f$ must be evaluated.

Output Parameters

1:     result – double scalar
The result of the function.
2:     region – function handle or string containing name of m-file
region must evaluate the limits of integration in any dimension.
[c, d] = region(ndim, x, j)

Input Parameters

1:     ndim – int64int32nag_int scalar
n$n$, the number of dimensions of the integral.
2:     x(ndim) – double array
x(1),,x(j1)${\mathbf{x}}\left(1\right),\dots ,{\mathbf{x}}\left(j-1\right)$ contain the current values of the first (j1)$\left(j-1\right)$ variables, which may be used if necessary in calculating cj${c}_{j}$ and dj${d}_{j}$.
3:     j – int64int32nag_int scalar
The index j$j$ for which the limits of the range of integration are required.

Output Parameters

1:     c – double scalar
The lower limit cj${c}_{j}$ of the range of xj${x}_{j}$.
2:     d – double scalar
The upper limit dj${d}_{j}$ of the range of xj${x}_{j}$.
3:     npts – int64int32nag_int scalar
The Korobov rule to be used. There are two alternatives depending on the value of npts.
 (i) 1 ≤ npts ≤ 6$1\le {\mathbf{npts}}\le 6$. In this case one of six preset rules is chosen using 2129$2129$, 5003$5003$, 10007$10007$, 20011$20011$, 40009$40009$ or 80021$80021$ points depending on the respective value of npts being 1$1$, 2$2$, 3$3$, 4$4$, 5$5$ or 6$6$. (ii) npts > 6${\mathbf{npts}}>6$. npts is the number of actual points to be used with corresponding optimal coefficients supplied in the array vk.
Constraint: npts1${\mathbf{npts}}\ge 1$.
4:     vk(ndim) – double array
ndim, the dimension of the array, must satisfy the constraint 1ndim20$1\le {\mathbf{ndim}}\le 20$.
If npts > 6${\mathbf{npts}}>6$, vk must contain the n$n$ optimal coefficients (which may be calculated using nag_quad_md_numth_coeff_prime (d01gy) or nag_quad_md_numth_coeff_2prime (d01gz)).
If npts6${\mathbf{npts}}\le 6$, vk need not be set.
5:     nrand – int64int32nag_int scalar
The number of random samples to be generated in the error estimation (generally a small value, say 3$3$ to 5$5$, is sufficient). The total number of integrand evaluations will be ${\mathbf{nrand}}×{\mathbf{npts}}$.
Constraint: nrand1${\mathbf{nrand}}\ge 1$.

### Optional Input Parameters

1:     ndim – int64int32nag_int scalar
Default: The dimension of the array vk.
n$n$, the number of dimensions of the integral.
Constraint: 1ndim20$1\le {\mathbf{ndim}}\le 20$.
2:     itrans – int64int32nag_int scalar
Indicates whether the periodising transformation is to be used.
itrans = '0'${\mathbf{itrans}}=\text{'0'}$
The transformation is to be used.
itrans'0'${\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).
Default: 0$0$

None.

### Output Parameters

1:     vk(ndim) – double array
If npts > 6${\mathbf{npts}}>6$, vk is unchanged.
If npts6${\mathbf{npts}}\le 6$, vk contains the n$n$ optimal coefficients used by the preset rule.
2:     res – double scalar
The approximation to the integral I$I$.
3:     err – double scalar
The standard error as computed from nrand sample values. If nrand = 1${\mathbf{nrand}}=1$, then err contains zero.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, ndim < 1${\mathbf{ndim}}<1$, or ndim > 20${\mathbf{ndim}}>20$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, npts < 1${\mathbf{npts}}<1$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, nrand < 1${\mathbf{nrand}}<1$.

## Accuracy

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 × p${\mathbf{nrand}}×p$, where p$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.

## Example

npts = int64(2);
vk = zeros(4,1);
nrand = int64(4);
[vkOut, res, err, ifail] = nag_quad_md_numth(@f, @region, npts, vk, nrand)

function result = f(ndim,x)
result = cos(0.5+2*sum(x)-double(ndim));
function [c,d] = region(ndim, x, j)
c=0;
d=1;

vkOut =

1
792
1889
191

res =

0.4400

err =

1.8894e-06

ifail =

0

function d01gc_example
npts = int64(2);
vk = zeros(4,1);
nrand = int64(4);
[vkOut, res, err, ifail] = d01gc(@f, @region, npts, vk, nrand)

function result = f(ndim,x)
result = cos(0.5+2*sum(x)-double(ndim));
function [c,d] = region(ndim, x, j)
c=0;
d=1;

vkOut =

1
792
1889
191

res =

0.4400

err =

1.8894e-06

ifail =

0