Note: before using this routine, please read the Usersâ€™ Note for your implementation to check the interpretation of bold italicised terms andother implementation-dependent details.
D01GCF calculates an approximation to a definite integral in up to 20 dimensions, using the Korobovâ€“Conroy number theoretic method.
SUBROUTINE D01GCF(NDIM, FUNCTN, REGION, NPTS, VK, NRAND, ITRANS, RES,
This routine calculates an approximation to the integral
using the Korobovâ€“Conroy number The region of integration defined such that generally ci and di may be functions ofx1; x2; . . . ; xiÃ€1, for i Â¼ 2; 3; . . . ; n, with c1 and d1 constants. The integral is first of all transformed to anintegral over the n-cube Â½0; 1ÂŠn by the change of variables
The method then uses as its basis the number theoretic formula for the n-cube, Â½0; 1ÂŠn:
where fxg denotes the fractional part of x, a1; a2; . . . ; an are the so-called optimal coefficients, E is theerror, and p is a prime integer. (It is strictly only necessary that p be relatively prime to all a1; a2; . . . ; anand is in fact chosen to be even use of properties ofthe Fourier expansion of gÃ°x1; x2; . . . ; xnÃž which is assumed to have some degree of periodicity. Depending on the choice of a1; a2; . . . ; an the contributions from certain groups of Fourier coefficients areeliminated from the error, E. Korobov shows that a1; a2; . . . ; an can be chosen so that the error satisfies
where and C are real numbers depending on the convergence rate of the Fourier series, is a constantdepending on n, and K is a constant depending on and n.
calculating these optimal coefficients. Korobov imposes the constraint that
and gives a procedure for calculating the parameter, a, to satisfy the optimal conditions.
In this routine the periodisation is achieved by the simple transformation
More sophisticated periodisation procedures are available but in practice the degree of periodisation doesnot appear to be a critical requirement of the method.
An easily calculable error estimate is not available apart from repetition with an increasing sequence ofvalues of p which can yield erratic results. The difficulties have been studied by Cranley and Pattersonproposed a Monte Carlo error estimate arising from a stochasticintegration rule by the inclusion of a random origin shift which leaves the form of the
Computing the integral for each of a sequence of random vectors allows a
â€˜standard errorâ€™ to be estimated.
This routine provides built-in sets of optimal coefficients, corresponding to six different values of p. Alternatively the optimal coefficients may be supplied by the user.
compute the optimal coefficients for the cases where p is a prime number or p is a product of 2 primes,respectively.
Korobov N M (1957) The approximate calculation of multiple integrals using number theoretic methodsDokl. Acad. Nauk SSSR 115 1062â€“1065
Korobov N M (1963) Number Theoretic Methods in Approximate Analysis Fizmatgiz, Moscow
Conroy H (1967) Molecular Shroedinger equation VIII. A new method for evaluting multi-dimensionalintegrals J. Chem. Phys. 47 5307â€“5318
Cranley R and Patterson T N L (1976) Randomisation of number theoretic methods for mulitple integrationSIAM J. Numer. Anal. 13 904â€“914
On entry: the number of dimensions of the integral, n.
FUNCTN â€“ real FUNCTION, supplied by the user.
FUNCTN must return the value of the integrand f at a given point.
On entry: the number of dimensions of the integral, n.
On entry: the co-ordinates of the point at which the integrand must be evaluated.
FUNCTN must be declared as EXTERNAL in the (sub)program from which D01GCF is called. Parameters denoted as Input must not be changed by this procedure.
REGION â€“ SUBROUTINE, supplied by the user.
REGION must evaluate the limits of integration in any dimension.
On entry: the number of dimensions of the integral, n.
On entry: XÃ°1Ãž; . . . ; XÃ°j Ã€ 1Ãž contain the current values of the first Ã°j Ã€ 1Ãž variables,which may be used if necessary in calculating cj and dj.
On entry: the index j for which the limits of the range of integration are required.
On exit: the lower limit cj of the range of xj.
On exit: the upper limit dj of the range of xj.
REGION must be declared as EXTERNAL in the (sub)program from which D01GCF is called. Parameters denoted as Input must not be changed by this procedure.
On entry: the Korobov rule to be used. There are two alternatives depending on the value of NPTS.
In this case one of six preset rules is chosen using 2129, 5003, 10007, 20011, 40009 or 80021points depending on the respective value of NPTS being 1, 2, 3, 4, 5 or 6.
NPTS is the number of actual points to be used with corresponding optimal coefficientssupplied in the
On entry: if 6, VK must contain the n optimal coefficients (which may be calculated using6, VK need not be set.
On exit: if 6, VK is unchanged; if 6, VK contains the n optimal coefficients usedby the preset rule.
On entry: the number of random samples to be generated in the error estimation (generally a smallvalue, say 3 to 5 is sufficient). The total number of integrand evaluations will be NRAND Ã‚
On entry: indicates whether the periodising transformation is to be used:
if ITRANS Â¼ 0, the transformation is to be used;
if ITRANS 6Â¼ 0, the transformation is to be suppresssed (to cover cases where the integrandmay already be periodic or where the user desires to specify a particular transformation in the
On exit: an estimate of the value of the integral.
On exit: the standard error as If 1, then ERRcontains zero.
On entry: IFAIL must be set to 0, Ã€1 or 1. Users who are unfamiliar with this parameter should
On exit: IFAIL Â¼ 0 unless the routine detects an
For environments where it might be inappropriate to halt program execution when an error isdetected, the value Ã€1 or 1 is recommended. If the output of error messages is undesirable, then thevalue 1 is recommended. Otherwise, for users not familiar with this parameter the recommendedvalue is 0. When the value Ã€1 or 1 is used it is essential to test the value of IFAIL on exit.
If on entry IFAIL Â¼ 0 or Ã€1, explanatory error messages are output on the current error message unit (as
Errors or warnings detected by the routine:
An estimate of the absolute standard error is given by the value, on
The time taken by the routine will be approximately proportional to p, where p is the numberof points used.
The exact by D01GCF will depend (within statistical limits) on thesequence of random numbers generated within the routine that the results
returned by D01GCF in separate runs are identical, users before callingD01GCF; to ensure that they are
cosÃ°0:5 Ã¾ 2Ã°x1 Ã¾ x2 Ã¾ x3 Ã¾ x4Ãž Ã€ 4Ãž dx1 dx2 dx3 dx4:
Note: the listing of the example program presented below uses bold italicised terms to denote precision-dependent details. Please read theUsersâ€™ Note for your implementation to check the interpretation of these the results produced may not be identical for all implementations.
. Executable Statements . WRITE (NOUT,*) â€™D01GCF Example Program Resultsâ€™NPTS = 2ITRANS = 0NRAND = 4IFAIL = 0
CALL D01GCF(NDIM,FUNCT,REGION,NPTS,VK,NRAND,ITRANS,RES,ERR,IFAIL)
WRITE (NOUT,*)WRITE (NOUT,99999) â€™Result =â€™, RES, â€™
. Executable Statements . A = 0.0e0B = 1.0e0RETURNEND
. Executable Statements . SUM = 0.0e0DO 20 J = 1, NDIM
FUNCT = COS(0.5e0+2.0e0*SUM-real(NDIM))RETURNEND

Advice on Replacement Calls for Withdrawn/Superseded Routines

Implementation-specific Information

C05 - Roots of One or More Transcendental Equations

D02 - Ordinary Differential Equations

D02 - Ordinary Differential Equations

E04 - Minimizing or Maximizing a Function

F06 - Linear Algebra Support Routines

F08 - Least-squares and Eigenvalue Problems (LAPACK)