NAG FL Interface
d01fdf (md_​sphere)

1 Purpose

d01fdf calculates an approximation to a definite integral in up to 30 dimensions, using the method of Sag and Szekeres (see Sag and Szekeres (1964)). The region of integration is an n-sphere, or by built-in transformation via the unit n-cube, any product region.

2 Specification

Fortran Interface
Subroutine d01fdf ( ndim, f, sigma, region, limit, r0, u, result, ncalls, ifail)
Integer, Intent (In) :: ndim, limit
Integer, Intent (Inout) :: ifail
Integer, Intent (Out) :: ncalls
Real (Kind=nag_wp), External :: f
Real (Kind=nag_wp), Intent (In) :: sigma, r0, u
Real (Kind=nag_wp), Intent (Out) :: result
External :: region
C Header Interface
#include <nag.h>
void  d01fdf_ (const Integer *ndim,
double (NAG_CALL *f)(const Integer *ndim, const double x[]),
const double *sigma,
void (NAG_CALL *region)(const Integer *ndim, const double x[], const Integer *j, double *c, double *d),
const Integer *limit, const double *r0, const double *u, double *result, Integer *ncalls, Integer *ifail)
The routine may be called by the names d01fdf or nagf_quad_md_sphere.

3 Description

d01fdf calculates an approximation to
n-sphere of radius  σ f x1,x2,,xn dx1 dx2 dxn (1)
or, more generally,
c1 d1 dx1 cn dn dxn f x1,,xn (2)
where each ci and di may be functions of xj j<i.
The routine uses the method of Sag and Szekeres (1964), which exploits a property of the shifted p-point trapezoidal rule, namely, that it integrates exactly all polynomials of degree <p (see Krylov (1962)). An attempt is made to induce periodicity in the integrand by making a parameterised transformation to the unit n-sphere. The Jacobian of the transformation and all its direct derivatives vanish rapidly towards the surface of the unit n-sphere, so that, except for functions which have strong singularities on the boundary, the resulting integrand will be pseudo-periodic. In addition, the variation in the integrand can be considerably reduced, causing the trapezoidal rule to perform well.
Integrals of the form (1) are transformed to the unit n-sphere by the change of variables:
xi = yi σr tanh ur 1-r2  
where r2=i=1nyi2 and u is an adjustable parameter.
Integrals of the form (2) are first of all transformed to the n-cube -1,1n by a linear change of variables
xi=di+ci+di-ciyi/2  
and then to the unit sphere by a further change of variables
yi=tanhuzi 1-r  
where r2=i=1nzi2 and u is again an adjustable parameter.
The parameter u in these transformations determines how the transformed integrand is distributed between the origin and the surface of the unit n-sphere. A typical value of u is 1.5. For larger u, the integrand is concentrated toward the centre of the unit n-sphere, while for smaller u it is concentrated toward the perimeter.
In performing the integration over the unit n-sphere by the trapezoidal rule, a displaced equidistant grid of size h is constructed. The points of the mesh lie on concentric layers of radius
ri=h4n+8i-1,  i=1,2,3,.  
The routine requires you to specify an approximate maximum number of points to be used, and then computes the largest number of whole layers to be used, subject to an upper limit of 400 layers.
In practice, the rapidly-decreasing Jacobian makes it unnecessary to include the whole unit n-sphere and the integration region is limited by a user-specified cut-off radius r0<1. The grid-spacing h is determined by r0 and the number of layers to be used. A typical value of r0 is 0.8.
Some experimentation may be required with the choice of r0 (which determines how much of the unit n-sphere is included) and u (which determines how the transformed integrand is distributed between the origin and surface of the unit n-sphere), to obtain best results for particular families of integrals. This matter is discussed further in Section 9.

4 References

Krylov V I (1962) Approximate Calculation of Integrals (trans A H Stroud) Macmillan
Sag T W and Szekeres G (1964) Numerical evaluation of high-dimensional integrals Math. Comput. 18 245–253

5 Arguments

1: ndim Integer Input
On entry: n, the number of dimensions of the integral.
Constraint: 1ndim30.
2: f real (Kind=nag_wp) Function, supplied by the user. External Procedure
f must return the value of the integrand f at a given point.
The specification of f is:
Fortran Interface
Function f ( ndim, x)
Real (Kind=nag_wp) :: f
Integer, Intent (In) :: ndim
Real (Kind=nag_wp), Intent (In) :: x(ndim)
C Header Interface
double  f_ (const Integer *ndim, const double x[])
1: ndim Integer Input
On entry: n, the number of dimensions of the integral.
2: xndim Real (Kind=nag_wp) array Input
On entry: the coordinates of the point at which the integrand f must be evaluated.
f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01fdf is called. Arguments denoted as Input must not be changed by this procedure.
Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01fdf. If your code inadvertently does return any NaNs or infinities, d01fdf is likely to produce unexpected results.
3: sigma Real (Kind=nag_wp) Input
On entry: indicates the region of integration.
sigma0.0
The integration is carried out over the n-sphere of radius sigma, centred at the origin.
sigma<0.0
The integration is carried out over the product region described by region.
4: region Subroutine, supplied by the NAG Library or the user. External Procedure
If sigma<0.0, region must evaluate the limits of integration in any dimension.
The specification of region is:
Fortran Interface
Subroutine region ( ndim, x, j, c, d)
Integer, Intent (In) :: ndim, j
Real (Kind=nag_wp), Intent (In) :: x(ndim)
Real (Kind=nag_wp), Intent (Out) :: c, d
C Header Interface
void  region_ (const Integer *ndim, const double x[], const Integer *j, double *c, double *d)
1: ndim Integer Input
On entry: n, the number of dimensions of the integral.
2: xndim Real (Kind=nag_wp) array Input
On entry: x1,,xj-1 contain the current values of the first j-1 variables, which may be used if necessary in calculating cj and dj.
3: j Integer Input
On entry: the index j for which the limits of the range of integration are required.
4: c Real (Kind=nag_wp) Output
On exit: the lower limit cj of the range of xj.
5: d Real (Kind=nag_wp) Output
On exit: the upper limit dj of the range of xj.
region must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01fdf is called. Arguments denoted as Input must not be changed by this procedure.
Note: region should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01fdf. If your code inadvertently does return any NaNs or infinities, d01fdf is likely to produce unexpected results.
If sigma0.0, region is not called by d01fdf, but a dummy routine must be supplied (d01fdv may be used).
5: limit Integer Input
On entry: the approximate maximum number of integrand evaluations to be used.
Constraint: limit100.
6: r0 Real (Kind=nag_wp) Input
On entry: the cut-off radius on the unit n-sphere, which may be regarded as an adjustable parameter of the method.
Suggested value: a typical value is r0=0.8. (See also Section 9.)
Constraint: 0.0<r0<1.0.
7: u Real (Kind=nag_wp) Input
On entry: must specify an adjustable parameter of the transformation to the unit n-sphere.
Suggested value: a typical value is u=1.5. (See also Section 9.)
Constraint: u>0.0.
8: result Real (Kind=nag_wp) Output
On exit: the approximation to the integral I.
9: ncalls Integer Output
On exit: the actual number of integrand evaluations used. (See also Section 9.)
10: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1 or 1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this argument, the recommended value is 0. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, ndim=value.
Constraint: ndim30.
On entry, ndim=value.
Constraint: ndim1.
ifail=2
On entry, limit=value.
Constraint: limit100.
ifail=3
On entry, r0=value.
Constraint: r0<1.0.
On entry, r0=value.
Constraint: r0>0.0.
ifail=4
On entry, u=value.
Constraint: u>0.0.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

No error estimate is returned, but results may be verified by repeating with an increased value of limit (provided that this causes an increase in the returned value of ncalls).

8 Parallelism and Performance

d01fdf is not threaded in any implementation.

9 Further Comments

The time taken by d01fdf will be approximately proportional to the returned value of ncalls, which, except in the circumstances outlined in (b) below, will be close to the given value of limit.
  1. (a)Choice of r0 and u
    If the chosen combination of r0 and u is too large in relation to the machine accuracy it is possible that some of the points generated in the original region of integration may transform into points in the unit n-sphere which lie too close to the boundary surface to be distinguished from it to machine accuracy (despite the fact that r0<1). To be specific, the combination of r0 and u is too large if
    ur0 1-r02 >0.3465t-1,   if ​sigma0.0,  
    or
    ur0 1-r0 > 0.3465t- 1,   if ​ sigma< 0.0,  
    where t is the number of bits in the mantissa of a real number.
    The contribution of such points to the integral is neglected. This may be justified by appeal to the fact that the Jacobian of the transformation rapidly approaches zero towards the surface. Neglect of these points avoids the occurrence of overflow with integrands which are infinite on the boundary.
  2. (b)Values of limit and ncalls
    limit is an approximate upper limit to the number of integrand evaluations, and may not be chosen less than 100. There are two circumstances when the returned value of ncalls (the actual number of evaluations used) may be significantly less than limit.
    Firstly, as explained in (a), an unsuitably large combination of r0 and u may result in some of the points being unusable. Such points are not included in the returned value of ncalls.
    Secondly, no more than 400 layers will ever be used, no matter how high limit is set. This places an effective upper limit on ncalls as follows:
    n=1: 56 n=2: 1252 n=3: 23690 n=4: 394528 n=5: 5956906  

10 Example

This example calculates the integral
s dx1dx2dx3 σ2-r2 =22.2066  
where s is the 3-sphere of radius σ, r2=x12+x22+x32 and σ=1.5. Both sphere-to-sphere and general product region transformations are used. For the former, we use r0=0.9 and u=1.5; for the latter, r0=0.8 and u=1.5.

10.1 Program Text

Program Text (d01fdfe.f90)

10.2 Program Data

None.

10.3 Program Results

Program Results (d01fdfe.r)