Syntax

Description

nag_quad_md_numth (d01gc) calculates an approximation to the integral

d1

dn

I =

∫

dx1, … ,

∫

dxnf(x1,x2, … ,xn)

c1

cn

I=∫c1d1dx1,…,∫cndndxnf(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 cici and didi may be functions of x1,x2, … ,xi − 1x1,x2,…,xi-1, for i = 2,3, … ,ni=2,3,…,n, with c1c1 and d1d1 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 + (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 nn-cube, [0,1]n[0,1]n:

1

1

p

∫

dx1 ⋯

∫

dxng(x1,x2, … ,xn) = 1/p

∑

g({k(a1)/p}, … ,{k(an)/p}) − E

0

0

k = 1

∫01dx1⋯∫01dxng(x1,x2,…,xn)=1p∑k=1pg({ka1p},…,{kanp})-E

(2)

where {x}{x} denotes the fractional part of xx, a1,a2, … ,ana1,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, … ,ana1,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, … ,ana1,a2,…,an the contributions from certain groups of Fourier coefficients are eliminated from the error, EE. Korobov shows that a1,a2, … ,ana1,a2,…,an can be chosen so that the error satisfies

E ≤ CKp − αlnαβp

E≤CKp-α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 = ai − 1(modp)

a1=1 and ai=ai-1(modp)

(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(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 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}{kaip} is replaced by {αi + k(ai)/p}{αi+kaip}, for i = 1,2, … ,ni=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.

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.

Error Indicators and Warnings

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.