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_vec (d01gd)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_quad_md_numth_vec (d01gd) calculates an approximation to a definite integral in up to 20 dimensions, using the Korobov–Conroy number theoretic method. This function is designed to be particularly efficient on vector processors.

Syntax

[vk, res, err, ifail] = d01gd(vecfun, vecreg, npts, vk, nrand, 'ndim', ndim, 'itrans', itrans)
[vk, res, err, ifail] = nag_quad_md_numth_vec(vecfun, vecreg, npts, vk, nrand, 'ndim', ndim, 'itrans', itrans)

Description

nag_quad_md_numth_vec (d01gd) calculates an approximation to the integral
I= c1 d1 cn dn f x1,,xn dxn dx1 (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  and di  may be functions of x1 , x2 ,, xi-1 , for i= 2 , 3 ,, n , with c1  and d1  constants. The integral is first of all transformed to an integral over the n-cube 0,1 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, 0,1 n :
01 01 g x1,,xn dxn dx1 = 1p k=1p g k a1p ,, k anp - E (2)
where x  denotes the fractional part of x, a1 ,, an  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 a1 ,, 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,,xn  which is assumed to have some degree of periodicity. Depending on the choice of a1 ,, an  the contributions from certain groups of Fourier coefficients are eliminated from the error, E. Korobov shows that a1 ,, an  can be chosen so that the error satisfies
ECK 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
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), k ai p  is replaced by αi+k ai p , for i= 1 , 2 ,, 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.
This function is designed to be particularly efficient on vector processors, although it is very important that you also code vecfun and vecreg efficiently.

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:     vecfun – function handle or string containing name of m-file
vecfun must evaluate the integrand at a specified set of points.
[fv] = vecfun(ndim, x, m)

Input Parameters

1:     ndim int64int32nag_int scalar
n, the number of dimensions of the integral.
2:     xmndim – double array
The coordinates of the m points at which the integrand must be evaluated. xij contains the jth coordinate of the ith point.
3:     m int64int32nag_int scalar
The number of points m at which the integrand is to be evaluated.

Output Parameters

1:     fvm – double array
fvi must contain the value of the integrand of the ith point, i.e., fvi=f xi1,xi2,,xindim, for i=1,2,,m.
2:     vecreg – function handle or string containing name of m-file
vecreg must evaluate the limits of integration in any dimension for a set of points.
[c, d] = vecreg(ndim, x, j, m)

Input Parameters

1:     ndim int64int32nag_int scalar
n, the number of dimensions of the integral.
2:     xmndim – double array
For i=1,2,,m, xi1, xi2,,xij-1 contain the current values of the first j-1 coordinates of the ith point, which may be used if necessary in calculating the m values of cj and dj.
3:     j int64int32nag_int scalar
The index j for which the limits of the range of integration are required.
4:     m int64int32nag_int scalar
The number of points m at which the limits of integration must be specified.

Output Parameters

1:     cm – double array
ci must be set to the lower limit of the range for xij, for i=1,2,,m.
2:     dm – double array
di must be set to the upper limit of the range for xij, for i=1,2,,m.
3:     npts int64int32nag_int scalar
The Korobov rule to be used. There are two alternatives depending on the value of npts.
(i) 1npts6.
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: npts1.
4:     vkndim – double array
If 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 npts6, vk need not be set.
5:     nrand int64int32nag_int scalar
The number of random samples to be generated (generally a small value, say 3 to 5, is sufficient). The estimate, res, of the value of the integral returned by the function is then the average of nrand calculations with different random origin shifts. If npts>6, the total number of integrand evaluations will be nrand×npts. If 1npts6, then the number of integrand evaluations will be nrand×p, where p is the number of points corresponding to the six preset rules. For reasons of efficiency, these values are calculated a number at a time in vecfun.
Constraint: nrand1.

Optional Input Parameters

1:     ndim int64int32nag_int scalar
Default: the dimension of the array vk.
n, the number of dimensions of the integral.
Constraint: 1ndim20.
2:     itrans int64int32nag_int scalar
Default: 0
Indicates whether the periodising transformation is to be used.
itrans=0
The transformation is to be used.
itrans0
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 vecfun).

Output Parameters

1:     vkndim – double array
If npts>6, vk is unchanged.
If npts6, vk contains the n optimal coefficients used by the preset rule.
2:     res – double scalar
The approximation to the integral I.
3:     err – double scalar
The standard error as computed from nrand sample values. If nrand=1, then err contains zero.
4:     ifail int64int32nag_int scalar
ifail=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
On entry,ndim<1,
orndim>20.
   ifail=2
On entry,npts<1.
   ifail=3
On entry,nrand<1.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

Accuracy

If nrand>1, an estimate of the absolute standard error is given by the value, on exit, of err.

Further Comments

nag_quad_md_numth_vec (d01gd) performs the same computation as nag_quad_md_numth (d01gc). However, the interface has been modified so that it can perform more efficiently on machines with vector processing capabilities. In particular, vecfun and vecreg must calculate the integrand and limits of integration at a set of points. For some problems the amount of time spent in these two functions, which must be supplied by you, may account for a significant part of the total computation time. For this reason it is vital that you consider the possibilities for vectorization in the code supplied for these two functions.
The time taken will be approximately proportional to nrand×p, where p is the number of points used, but may depend significantly on the efficiency of the code provided by you in vecfun and vecreg.
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_vec (d01gd) 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 d01gd_example


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

npts = int64(2);
vk = zeros(4,1);
nrand = int64(4);

[vk, res, err, ifail] = d01gd( ...
			       @vecfun, @vecreg, npts, vk, nrand);

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



function fv = vecfun(ndim, x, m)
  fv = zeros(m,1);

  for i=1:m
    fv(i) = cos(0.5 + 2*sum(x(i,1:ndim)) - double(ndim));
  end

function [c,d] = vecreg(ndim, x, j, m)
  c = zeros(m,1);
  d = ones(m,1);
d01gd example results

Result =       0.43999,  Standard error =   1.89e-06

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–2015