3 Description

nag_quad_md_numth_vec (d01gdc) calculates an approximation to the integral

I=∫c1d1⋯∫cndnfx1,…,xndxn…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,1n by the change of variables

xi=ci+di-ciyi, i=1,2,…,n.

The method then uses as its basis the number theoretic formula for the n-cube, 0,1n:

∫01⋯∫01gx1,…,xndxn⋯dx1=1p∑k=1pgka1p,…,kanp-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 gx1,…,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

E≤CKp-α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-1modp

(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=yi23-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), kaip is replaced by αi+kaip, 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.

On exit: fv[i-1] must contain the value of the integrand of the ith point, i.e., fv[i-1]=fXi,1,Xi,2,…,Xi,ndim, for i=1,2,…,m.

4:
m – IntegerInput

On entry:
the number of points m at which the integrand is to be evaluated.

5:
comm – Nag_Comm *

Pointer to structure of type Nag_Comm; the following members are relevant to vecfun.

user – double *

iuser – Integer *

p – Pointer

The type Pointer will be void *. Before calling nag_quad_md_numth_vec (d01gdc) you may allocate memory and initialize these pointers with various quantities for use by vecfun when called from nag_quad_md_numth_vec (d01gdc) (see Section 3.2.1 in the Essential Introduction).

3:
vecreg – function, supplied by the userExternal Function

vecreg must evaluate the limits of integration in any dimension for a set of points.

On exit:
d[i-1] must be set to the upper limit of the range for Xi,j, for i=1,2,…,m.

6:
m – IntegerInput

On entry:
the number of points m at which the limits of integration must be specified.

7:
comm – Nag_Comm *

Pointer to structure of type Nag_Comm; the following members are relevant to vecreg.

user – double *

iuser – Integer *

p – Pointer

The type Pointer will be void *. Before calling nag_quad_md_numth_vec (d01gdc) you may allocate memory and initialize these pointers with various quantities for use by vecreg when called from nag_quad_md_numth_vec (d01gdc) (see Section 3.2.1 in the Essential Introduction).

4:
npts – IntegerInput

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

(i)

1≤npts≤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.

If npts≤6, vk contains the n optimal coefficients used by the preset rule.

6:
nrand – IntegerInput

On entry: 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 1≤npts≤6, 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:
nrand≥1.

7:
transform – Nag_BooleanInput

On entry: indicates whether the periodising transformation is to be used.

transform=Nag_TRUE

The transformation is to be used.

transform=Nag_FALSE

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).

Suggested value:
transform=Nag_TRUE.

8:
res – double *Output

On exit: the approximation to the integral I.

9:
err – double *Output

On exit: the standard error as computed from nrand sample values. If nrand=1, then err contains zero.

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.

7 Accuracy

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

8 Further Comments

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.

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 (d01gdc) by calls to nag_rand_basic (g05sac). Separate runs will produce identical answers.