hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_quad_md_numth (d01gc)

Purpose

nag_quad_md_numth (d01gc) calculates an approximation to a definite integral in up to 2020 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,,dxnf(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 ci  and di di  may be functions of x1 , x2 ,, xi1 x1 , x2 ,, xi-1 , for i = 2 , 3 ,, n i= 2 , 3 ,, n , with c1 c1  and d1 d1  constants. The integral is first of all transformed to an integral over the nn-cube [0,1]n [0,1] n  by the change of variables
xi = ci + (dici) 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 nn-cube, [0,1]n [0,1] n :
1 1 p
dx1dxng(x1,x2,,xn) = 1/pg({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} {x}  denotes the fractional part of xx, a1 , a2 ,, an a1 , a2 ,, an  are the so-called optimal coefficients, EE is the error, and pp is a prime integer. (It is strictly only necessary that pp be relatively prime to all a1 , a2 ,, an a1 , a2 ,, an  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 (x1,x2,,xn)  which is assumed to have some degree of periodicity. Depending on the choice of a1 , a2 ,, an a1 , a2 ,, an  the contributions from certain groups of Fourier coefficients are eliminated from the error, EE. Korobov shows that a1 , a2 ,, an a1 , a2 ,, an  can be chosen so that the error satisfies
ECK pα lnαβ p
ECK p-α ln αβ p
(3)
where αα and CC are real numbers depending on the convergence rate of the Fourier series, ββ is a constant depending on nn, and KK is a constant depending on αα and nn. There are a number of procedures for calculating these optimal coefficients. Korobov imposes the constraint that
a1 = 1   and   ai = ai1   (mod p)
a1 = 1   and   ai = ai-1 (mod p)
(4)
and gives a procedure for calculating the parameter, aa, to satisfy the optimal conditions.
In this function the periodisation is achieved by the simple transformation
xi = yi2 (32yi) ,   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 pp 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} { k ai p }  is replaced by {αi + k(ai)/p} { αi+k ai p } , for i = 1 , 2 ,, n i= 1 , 2 ,, n , where each αi αi , is uniformly distributed over [0,1] [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 pp. 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 pp is a prime number or pp 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 ff at a given point.
[result] = f(ndim, x)

Input Parameters

1:     ndim – int64int32nag_int scalar
nn, the number of dimensions of the integral.
2:     x(ndim) – double array
The coordinates of the point at which the integrand ff 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
nn, the number of dimensions of the integral.
2:     x(ndim) – double array
x(1),,x(j1)x1,,xj-1 contain the current values of the first (j1)(j-1) variables, which may be used if necessary in calculating cjcj and djdj.
3:     j – int64int32nag_int scalar
The index jj for which the limits of the range of integration are required.

Output Parameters

1:     c – double scalar
The lower limit cjcj of the range of xjxj.
2:     d – double scalar
The upper limit djdj of the range of xjxj.
3:     npts – int64int32nag_int scalar
The Korobov rule to be used. There are two alternatives depending on the value of npts.
(i) 1npts61npts6.
In this case one of six preset rules is chosen using 21292129, 50035003, 1000710007, 2001120011, 4000940009 or 8002180021 points depending on the respective value of npts being 11, 22, 33, 44, 55 or 66.
(ii) npts > 6npts>6.
npts is the number of actual points to be used with corresponding optimal coefficients supplied in the array vk.
Constraint: npts1npts1.
4:     vk(ndim) – double array
ndim, the dimension of the array, must satisfy the constraint 1ndim201ndim20.
If npts > 6npts>6, vk must contain the nn optimal coefficients (which may be calculated using nag_quad_md_numth_coeff_prime (d01gy) or nag_quad_md_numth_coeff_2prime (d01gz)).
If npts6npts6, 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 33 to 55, is sufficient). The total number of integrand evaluations will be nrand × nptsnrand×npts.
Constraint: nrand1nrand1.

Optional Input Parameters

1:     ndim – int64int32nag_int scalar
Default: The dimension of the array vk.
nn, the number of dimensions of the integral.
Constraint: 1ndim201ndim20.
2:     itrans – int64int32nag_int scalar
Indicates whether the periodising transformation is to be used.
itrans = '0'itrans='0'
The transformation is to be used.
itrans'0'itrans'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: 00

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

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

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,ndim < 1ndim<1,
orndim > 20ndim>20.
  ifail = 2ifail=2
On entry,npts < 1npts<1.
  ifail = 3ifail=3
On entry,nrand < 1nrand<1.

Accuracy

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

Further Comments

The time taken by nag_quad_md_numth (d01gc) will be approximately proportional to nrand × pnrand×p, where pp 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

function nag_quad_md_numth_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



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013