NAG FL Interface
d01fdf (md_sphere)
1
Purpose
d01fdf calculates an approximation to a definite integral in up to
dimensions, using the method of Sag and Szekeres (see
Sag and Szekeres (1964)). The region of integration is an
-sphere, or by built-in transformation via the unit
-cube, any product region.
2
Specification
Fortran Interface
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) |
|
C++ Header Interface
#include <nag.h> extern "C" {
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
or, more generally,
where each
and
may be functions of
.
The routine uses the method of
Sag and Szekeres (1964), which exploits a property of the shifted
-point trapezoidal rule, namely, that it integrates exactly all polynomials of degree
(see
Krylov (1962)). An attempt is made to induce periodicity in the integrand by making a parameterised transformation to the unit
-sphere. The Jacobian of the transformation and all its direct derivatives vanish rapidly towards the surface of the unit
-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
-sphere by the change of variables:
where
and
is an adjustable parameter.
Integrals of the form
(2) are first of all transformed to the
-cube
by a linear change of variables
and then to the unit sphere by a further change of variables
where
and
is again an adjustable parameter.
The parameter in these transformations determines how the transformed integrand is distributed between the origin and the surface of the unit -sphere. A typical value of is . For larger , the integrand is concentrated toward the centre of the unit -sphere, while for smaller it is concentrated toward the perimeter.
In performing the integration over the unit
-sphere by the trapezoidal rule, a displaced equidistant grid of size
is constructed. The points of the mesh lie on concentric layers of radius
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
layers.
In practice, the rapidly-decreasing Jacobian makes it unnecessary to include the whole unit -sphere and the integration region is limited by a user-specified cut-off radius . The grid-spacing is determined by and the number of layers to be used. A typical value of is .
Some experimentation may be required with the choice of
(which determines how much of the unit
-sphere is included) and
(which determines how the transformed integrand is distributed between the origin and surface of the unit
-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:
– Integer
Input
-
On entry: , the number of dimensions of the integral.
Constraint:
.
-
2:
– real (Kind=nag_wp) Function, supplied by the user.
External Procedure
-
f must return the value of the integrand
at a given point.
The specification of
f is:
Fortran Interface
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[]) |
|
C++ Header Interface
#include <nag.h> extern "C" {
double |
f_ (const Integer &ndim, const double x[]) |
}
|
-
1:
– Integer
Input
-
On entry: , the number of dimensions of the integral.
-
2:
– Real (Kind=nag_wp) array
Input
-
On entry: the coordinates of the point at which the integrand 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:
– Real (Kind=nag_wp)
Input
-
On entry: indicates the region of integration.
- The integration is carried out over the -sphere of radius sigma, centred at the origin.
- The integration is carried out over the product region described by region.
-
4:
– Subroutine, supplied by the NAG Library or the user.
External Procedure
-
If
,
region must evaluate the limits of integration in any dimension.
The specification of
region is:
Fortran Interface
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) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
region_ (const Integer &ndim, const double x[], const Integer &j, double &c, double &d) |
}
|
-
1:
– Integer
Input
-
On entry: , the number of dimensions of the integral.
-
2:
– Real (Kind=nag_wp) array
Input
-
On entry: contain the current values of the first variables, which may be used if necessary in calculating and .
-
3:
– Integer
Input
-
On entry: the index for which the limits of the range of integration are required.
-
4:
– Real (Kind=nag_wp)
Output
-
On exit: the lower limit of the range of .
-
5:
– Real (Kind=nag_wp)
Output
-
On exit: the upper limit of the range of .
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
,
region is not called by
d01fdf,
but a dummy routine must be supplied (
d01fdv may be used).
-
5:
– Integer
Input
-
On entry: the approximate maximum number of integrand evaluations to be used.
Constraint:
.
-
6:
– Real (Kind=nag_wp)
Input
-
On entry: the cut-off radius on the unit -sphere, which may be regarded as an adjustable parameter of the method.
Suggested value:
a typical value is
. (See also
Section 9.)
Constraint:
.
-
7:
– Real (Kind=nag_wp)
Input
-
On entry: must specify an adjustable parameter of the transformation to the unit -sphere.
Suggested value:
a typical value is
. (See also
Section 9.)
Constraint:
.
-
8:
– Real (Kind=nag_wp)
Output
-
On exit: the approximation to the integral .
-
9:
– Integer
Output
-
On exit: the actual number of integrand evaluations used. (See also
Section 9.)
-
10:
– Integer
Input/Output
-
On entry:
ifail must be set to
,
or
to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of means that an error message is printed while a value of means that it is not.
If halting is not appropriate, the value
or
is recommended. If message printing is undesirable, then the value
is recommended. Otherwise, the value
is recommended.
When the value or is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
On entry, .
Constraint: .
On entry, .
Constraint: .
-
On entry, .
Constraint: .
-
On entry, .
Constraint: .
On entry, .
Constraint: .
-
On entry, .
Constraint: .
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.
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.
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.
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.
-
(a)Choice of and
If the chosen combination of
and
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
-sphere which lie too close to the boundary surface to be distinguished from it to machine accuracy (despite the fact that
). To be specific, the combination of
and
is too large if
or
where
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.
-
(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
. 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
and
may result in some of the points being unusable. Such points are not included in the returned value of
ncalls.
Secondly, no more than
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
is the
-sphere of radius
,
and
. Both sphere-to-sphere and general product region transformations are used. For the former, we use
and
; for the latter,
and
.
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results