NAG Library Function Document
nag_quad_md_sphere (d01fdc)
1
Purpose
nag_quad_md_sphere (d01fdc) 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 builtin transformation via the unit
$n$cube, any product region.
2
Specification
#include <nag.h> 
#include <nagd01.h> 
void 
nag_quad_md_sphere (Integer ndim,
double 
(*f)(Integer ndim,
const double x[],
Nag_Comm *comm),


double sigma,
void 
(*region)(Integer ndim,
const double x[],
Integer j,
double *c,
double *d,
Nag_Comm *comm),


Integer limit,
double r0,
double u,
double *result,
Integer *ncalls,
Nag_Comm *comm,
NagError *fail) 

3
Description
nag_quad_md_sphere (d01fdc) calculates an approximation to
or, more generally,
where each
${c}_{i}$ and
${d}_{i}$ may be functions of
${x}_{j}$ $\left(j<i\right)$.
The function 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
$\text{}<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 pseudoperiodic. 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:
where
${r}^{2}={\displaystyle \sum _{i=1}^{n}}{y}_{i}^{2}$ and
$u$ is an adjustable parameter.
Integrals of the form
(2) are first of all transformed to the
$n$cube
${\left[1,1\right]}^{n}$ by a linear change of variables
and then to the unit sphere by a further change of variables
where
${r}^{2}={\displaystyle \sum _{i=1}^{n}}{z}_{i}^{2}$ 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
The function 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 rapidlydecreasing Jacobian makes it unnecessary to include the whole unit $n$sphere and the integration region is limited by a userspecified cutoff radius ${r}_{0}<1$. The gridspacing $h$ is determined by ${r}_{0}$ and the number of layers to be used. A typical value of ${r}_{0}$ is $0.8$.
Some experimentation may be required with the choice of
${r}_{0}$ (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 highdimensional integrals Math. Comput. 18 245–253
5
Arguments
 1:
$\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
Constraint:
$1\le {\mathbf{ndim}}\le 30$.
 2:
$\mathbf{f}$ – function, supplied by the userExternal Function

f must return the value of the integrand
$f$ at a given point.
The specification of
f is:
double 
f (Integer ndim,
const double x[],
Nag_Comm *comm)


 1:
$\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
 2:
$\mathbf{x}\left[{\mathbf{ndim}}\right]$ – const doubleInput

On entry: the coordinates of the point at which the integrand $f$ must be evaluated.
 3:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
f.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
nag_quad_md_sphere (d01fdc) you may allocate memory and initialize these pointers with various quantities for use by
f when called from
nag_quad_md_sphere (d01fdc) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: f should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_quad_md_sphere (d01fdc). If your code inadvertently
does return any NaNs or infinities,
nag_quad_md_sphere (d01fdc) is likely to produce unexpected results.
 3:
$\mathbf{sigma}$ – doubleInput

On entry: indicates the region of integration.
 ${\mathbf{sigma}}\ge 0.0$
 The integration is carried out over the $n$sphere of radius sigma, centred at the origin.
 ${\mathbf{sigma}}<0.0$
 The integration is carried out over the product region described by region.
 4:
$\mathbf{region}$ – function, supplied by the userExternal Function

If
${\mathbf{sigma}}<0.0$,
region must evaluate the limits of integration in any dimension.
The specification of
region is:
void 
region (Integer ndim,
const double x[],
Integer j,
double *c,
double *d,
Nag_Comm *comm)


 1:
$\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
 2:
$\mathbf{x}\left[{\mathbf{ndim}}\right]$ – const doubleInput

On entry: ${\mathbf{x}}\left[0\right],\dots ,{\mathbf{x}}\left[j2\right]$ contain the current values of the first $\left(j1\right)$ variables, which may be used if necessary in calculating ${c}_{j}$ and ${d}_{j}$.
 3:
$\mathbf{j}$ – IntegerInput

On entry: the index $j$ for which the limits of the range of integration are required.
 4:
$\mathbf{c}$ – double *Output

On exit: the lower limit ${c}_{j}$ of the range of ${x}_{j}$.
 5:
$\mathbf{d}$ – double *Output

On exit: the upper limit ${d}_{j}$ of the range of ${x}_{j}$.
 6:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
region.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
nag_quad_md_sphere (d01fdc) you may allocate memory and initialize these pointers with various quantities for use by
region when called from
nag_quad_md_sphere (d01fdc) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: region should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_quad_md_sphere (d01fdc). If your code inadvertently
does return any NaNs or infinities,
nag_quad_md_sphere (d01fdc) is likely to produce unexpected results.
If
${\mathbf{sigma}}\ge 0.0$,
region is not called by
nag_quad_md_sphere (d01fdc),
but the
NAG defined null function pointer NULLFN must be supplied.
 5:
$\mathbf{limit}$ – IntegerInput

On entry: the approximate maximum number of integrand evaluations to be used.
Constraint:
${\mathbf{limit}}\ge 100$.
 6:
$\mathbf{r0}$ – doubleInput

On entry: the cutoff radius on the unit $n$sphere, which may be regarded as an adjustable parameter of the method.
Suggested value:
a typical value is
${\mathbf{r0}}=0.8$. (See also
Section 9.)
Constraint:
$0.0<{\mathbf{r0}}<1.0$.
 7:
$\mathbf{u}$ – doubleInput

On entry: must specify an adjustable parameter of the transformation to the unit $n$sphere.
Suggested value:
a typical value is
${\mathbf{u}}=1.5$. (See also
Section 9.)
Constraint:
${\mathbf{u}}>0.0$.
 8:
$\mathbf{result}$ – double *Output

On exit: the approximation to the integral $I$.
 9:
$\mathbf{ncalls}$ – Integer *Output

On exit: the actual number of integrand evaluations used. (See also
Section 9.)
 10:
$\mathbf{comm}$ – Nag_Comm *

The NAG communication argument (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
 11:
$\mathbf{fail}$ – NagError *Input/Output

The NAG error argument (see
Section 3.7 in How to Use the NAG Library and its Documentation).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_INT

On entry, ${\mathbf{limit}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{limit}}\ge 100$.
On entry, ${\mathbf{ndim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ndim}}\le 30$.
On entry, ${\mathbf{ndim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ndim}}\ge 1$.
 NE_INTERNAL_ERROR

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.
See
Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
 NE_REAL

On entry, ${\mathbf{r0}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{r0}}<1.0$.
On entry, ${\mathbf{r0}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{r0}}>0.0$.
On entry, ${\mathbf{u}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{u}}>0.0$.
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
nag_quad_md_sphere (d01fdc) is not threaded in any implementation.
The time taken by
nag_quad_md_sphere (d01fdc) 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.
(a) 
Choice of ${r}_{0}$ and $u$ If the chosen combination of ${r}_{0}$ 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 ${r}_{0}<1$). To be specific, the combination of ${r}_{0}$ and $u$ is too large if
or
where $t$ is the number of bits in the mantissa of a double 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. 
(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 ${r}_{0}$ 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:

10
Example
This example calculates the integral
where
$s$ is the
$3$sphere of radius
$\sigma $,
${r}^{2}={x}_{1}^{2}+{x}_{2}^{2}+{x}_{3}^{2}$ and
$\sigma =1.5$. Both spheretosphere and general product region transformations are used. For the former, we use
${r}_{0}=0.9$ and
$u=1.5$; for the latter,
${r}_{0}=0.8$ and
$u=1.5$.
10.1
Program Text
Program Text (d01fdce.c)
10.2
Program Data
None.
10.3
Program Results
Program Results (d01fdce.r)