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$ 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
 $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 argument, $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 .$
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 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.

## 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:     $\mathrm{f}$ – function handle or string containing name of m-file
f must return the value of the integrand $f$ at a given point.
[result] = f(ndim, x)

Input Parameters

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

Output Parameters

1:     $\mathrm{result}$ – double scalar
The value of $f\left(x\right)$ evaluated at x.
2:     $\mathrm{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:     $\mathrm{ndim}$int64int32nag_int scalar
$n$, the number of dimensions of the integral.
2:     $\mathrm{x}\left({\mathbf{ndim}}\right)$ – double array
${\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:     $\mathrm{j}$int64int32nag_int scalar
The index $j$ for which the limits of the range of integration are required.

Output Parameters

1:     $\mathrm{c}$ – double scalar
The lower limit ${c}_{j}$ of the range of ${x}_{j}$.
2:     $\mathrm{d}$ – double scalar
The upper limit ${d}_{j}$ of the range of ${x}_{j}$.
3:     $\mathrm{npts}$int64int32nag_int scalar
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$.
4:     $\mathrm{vk}\left({\mathbf{ndim}}\right)$ – double array
If ${\mathbf{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)).
If ${\mathbf{npts}}\le 6$, vk need not be set.
5:     $\mathrm{nrand}$int64int32nag_int scalar
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$.

### Optional Input Parameters

1:     $\mathrm{ndim}$int64int32nag_int scalar
Default: the dimension of the array vk.
$n$, the number of dimensions of the integral.
Constraint: $1\le {\mathbf{ndim}}\le 20$.
2:     $\mathrm{itrans}$int64int32nag_int scalar
Default: $0$
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).

### Output Parameters

1:     $\mathrm{vk}\left({\mathbf{ndim}}\right)$ – double array
If ${\mathbf{npts}}>6$, vk is unchanged.
If ${\mathbf{npts}}\le 6$, vk contains the $n$ optimal coefficients used by the preset rule.
2:     $\mathrm{res}$ – double scalar
The approximation to the integral $I$.
3:     $\mathrm{err}$ – double scalar
The standard error as computed from nrand sample values. If ${\mathbf{nrand}}=1$, then err contains zero.
4:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\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}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

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

## Example

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

fprintf('d01gc example results\n\n');

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

fprintf('Result = %13.5f,  Standard error = %10.2e\n', res, err);

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;
```
```d01gc example results

Result =       0.43999,  Standard error =   1.89e-06
```