), to a user-specified absolute or relative accuracy, by means of a modified Sag–Szekeres method. The routine can handle singularities on the surface or at the centre of the sphere, and returns an error estimate.

2 Specification

Fortran Interface

Subroutine d01jaf (

f, ndim, radius, epsa, epsr, method, icoord, result, esterr, nevals, ifail)

Integer, Intent (In)	::	ndim, method, icoord
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	nevals
Real (Kind=nag_wp), External	::	f
Real (Kind=nag_wp), Intent (In)	::	radius, epsa, epsr
Real (Kind=nag_wp), Intent (Out)	::	result, esterr

C Header Interface

#include <nag.h>

void	d01jaf_ ( double (NAG_CALL f)(const Integer ndim, const double x[]), const Integer ndim, const double radius, const double epsa, const double epsr, const Integer method, const Integer icoord, double result, double esterr, Integer nevals, Integer ifail)

The routine may be called by the names d01jaf or nagf_quad_md_sphere_bad.

3 Description

d01jaf calculates an approximation to the

n

-dimensional integral

I = \int \dots \int_{S} F (x_{1}, \dots, x_{n}) d x_{1} \dots d x_{n}, 2 \leq n \leq 4,

where

S

is the hypersphere

\sqrt{(x_{1}^{2} + \dots + x_{n}^{2})} \leq α < \infty

(the integrand function may also be defined in spherical coordinates). The algorithm is based on the Sag–Szekeres method (see Sag and Szekeres (1964)), applying the product trapezoidal formula after a suitable radial transformation. An improved transformation technique is developed: depending on the behaviour of the function and on the required accuracy, different transformations can be used, some of which are ‘double exponential’, as defined by Takahasi and Mori (1974). The resulting technique allows the routine to deal with integrand singularities on the surface or at the centre of the sphere. When the estimated error of the approximation with mesh size

h

is larger than the tolerated error, the trapezoidal formula with mesh size

h / 2

is calculated. A drawback of this method is the exponential growth of the number of function evaluations in the successive approximations (this number grows with a factor

\approx 2^{n}

). This introduces the restriction

n \leq 4

. Because the convergence rate of the successive approximations is normally better than linear, the error estimate is based on the linear extrapolation of the difference between the successive approximations (see Robinson and de Doncker (1981) and Roose and de Doncker (1981)). For further details of the algorithm, see Roose and de Doncker (1981).

4 References

Robinson I and de Doncker E (1981) Automatic computation of improper integrals over a bounded or unbounded planar region Computing 27 89–284

Roose D and de Doncker E (1981) Automatic integration over a sphere J. Comput. Appl. Math. 7 203–224

Sag T W and Szekeres G (1964) Numerical evaluation of high-dimensional integrals Math. Comput. 18 245–253

Takahasi H and Mori M (1974) Double Exponential Formulas for Numerical Integration 9 Publ. RIMS, Kyoto University 721–741

5 Arguments

1: $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: $x (ndim)$ – Real (Kind=nag_wp) array Input: On entry: the coordinates of the point at which the integrand $f$ must be evaluated. These coordinates are given in Cartesian or spherical polar form according to the value of icoord.

f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01jaf 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 d01jaf. If your code inadvertently does return any NaNs or infinities, d01jaf is likely to produce unexpected results.

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases d01jaf may return useful information.

$ifail = 1$: The required accuracy could not be achieved.

$ifail = 2$: The required accuracy could not be achieved.

$ifail = 3$: The required accuracy could not be achieved.
If $method = 0$ , $1$ or $2$ , setting $method = 3$ or $4$ may help.

$ifail = 4$: On entry, $icoord = ⟨ value ⟩$ .
Constraint: $icoord \leq 2$ .

On entry, $icoord = ⟨ value ⟩$ .
Constraint: $icoord \geq 0$ .

On entry, $icoord = 2$ and $method = ⟨ value ⟩$ .
Constraint: when $icoord = 2$ , $method = 3$ or $4$ .

On entry, $method = ⟨ value ⟩$ .
Constraint: $method \leq 6$ .

On entry, $method = ⟨ value ⟩$ .
Constraint: $method \geq 0$ .

On entry, $method = 3$ and $icoord = ⟨ value ⟩$ .
Constraint: when $method = 3$ , $icoord = 2$ .

On entry, $method = 4$ and $icoord = ⟨ value ⟩$ .
Constraint: when $method = 4$ , $icoord = 2$ .

On entry, $ndim = ⟨ value ⟩$ .
Constraint: $ndim \leq 4$ .

On entry, $ndim = ⟨ value ⟩$ .
Constraint: $ndim \geq 2$ .

On entry, $radius = ⟨ value ⟩$ .
Constraint: $radius \geq 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

You can specify an absolute and/or a relative tolerance, setting epsa and epsr. The routine attempts to calculate an approximation result such that

| I - result | \leq \max {epsa, epsr \times | I |} .

0 \leq ifail \leq 3

, esterr returns an estimate of, but not necessarily a bound for,

| I - result |

8 Parallelism and Performance

d01jaf is not threaded in any implementation.

9 Further Comments

9.1 Timing

Timing depends on the integrand and the accuracy required.

9.2 Spherical Polar Coordinates

Cartesian coordinates are related to the spherical polar coordinates by:

\begin{array}{cl} x_{1} & = & ρ . \sin θ_{1} \dots \sin θ_{n - 2} . \sin θ_{n - 1} \\ x_{2} & = & ρ . \sin θ_{1} \dots \sin θ_{n - 2} . \cos θ_{n - 1} \\ x_{3} & = & ρ . \sin θ_{1} \dots \cos θ_{n - 2} \\ ⋮ \\ x_{n} & = & ρ . \cos θ_{1} \end{array}

where

0 < θ_{i} < π

, for

i = 1, 2, \dots, n - 2

and

0 < θ_{n - 1} < 2 π

9.3 Machine Dependencies

As a consequence of the transformation technique, the severity of the singularities which can be handled by d01jaf depends on the precision and range of real numbers on the machine.

method = 3

4

must be used when the singularity on the surface is ‘severe’ in view of the requested accuracy and machine precision. In practice one has to set

method = 3

4

if d01jaf terminates with

ifail = 3

when called with

method = 0

1

2

When integrating a function with a severe singular behaviour on the surface of the sphere, the additional transformation

ρ = α - λ

helps to avoid the loss of significant figures due to round-off error in the calculation of the integration nodes which are very close to the surface. For these points, the value of

λ

can be computed more accurately than the value of

ρ

. Naturally, care must be taken that the function subprogram does not contain expressions of the form

α - λ

, which could cause a large round-off error in the calculation of the integrand at the boundary of the sphere.

Care should be taken to avoid underflow and/or overflow problems in the function subprogram, because some of the integration nodes used by d01jaf may be very close to the surface or to the centre of the sphere.

Example:

suppose the function

$f (ρ) = {(1 - ρ^{2})}^{- 0.7}$

is to be integrated over the unit sphere, with $method = 3$ or $4$ . Then icoord should be set to 2; the transformation $ρ = 1 - λ$ gives f(ρ)=(2λ-λ2) -0.7; and f could be coded thus:
```
f = 1.0
a = x(1)
If (a.gt.0.0) f = 1.0/(a*(2.0-a))**0.7
Return
```

Note that d01jaf ensures that

λ = x (1) > x02amf

, but underflow could occur in the computation of

λ^{2}

10 Example

This example evaluates the integrals

\int \dots \int_{S} \frac{1}{\sqrt{1 - ρ^{2}}} d x_{1} \dots d x_{n}

where

ρ = \sqrt{\sum_{i = 1}^{n} x_{i}^{2}}

, and

S

is the unit sphere of dimension

n = 2

4

The exact values (to

12

decimal places) are

6.28318530718

and

13.1594725348

10.1 Program Text

Program Text (d01jafe.f90)

10.2 Program Data

None.

10.3 Program Results

Program Results (d01jafe.r)

NAG Library Manual, Mark 27.3

Interfaces: FL CL CPP AD

NAG FL Interface Introduction

D01 (Quad) Chapter Contents

D01 (Quad) Chapter Introduction

d01ja: FL CL CPP AD

NAG FL Interfaced01jaf (md_​sphere_​bad)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 Spherical Polar Coordinates

9.3 Machine Dependencies

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
d01jaf (md_sphere_bad)