Integer type: int32 int64 nag_int show int32 show int32 show int64 show int64 show nag_int show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox

NAG Toolbox: nag_quad_md_sphere (d01fd)

▸▿ Contents

1 Purpose

2 Syntax

3 Description

4 References

▸▿ 5 Parameters

5.1 Compulsory Input Parameters

5.2 Optional Input Parameters

5.3 Output Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

Purpose

nag_quad_md_sphere (d01fd) 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.

Syntax

[result, ncalls, ifail] = d01fd(ndim, f, sigma, region, limit, 'r0', r0, 'u', u)

[result, ncalls, ifail] = nag_quad_md_sphere(ndim, f, sigma, region, limit, 'r0', r0, 'u', u)

Description

nag_quad_md_sphere (d01fd) calculates an approximation to

\int_{n -sphere of radius σ} f (x_{1}, x_{2}, \dots, x_{n}) d x_{1} d x_{2} \dots d x_{n}

(1)

or, more generally,

\int_{c_{1}}^{d_{1}} d x_{1} \dots \int_{c_{n}}^{d_{n}} d x_{n} f (x_{1}, \dots, x_{n})

(2)

where each

c_{i}

and

d_{i}

may be functions of

x_{j}

(j < i)

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

< 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:

x_{i} = y_{i} \frac{σ}{r} \tanh (\frac{u r}{1 - r^{2}})

where

r^{2} = \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

{[- 1, 1]}^{n}

by a linear change of variables

x_{i} = ((d_{i} + c_{i}) + (d_{i} - c_{i}) y_{i}) / 2

and then to the unit sphere by a further change of variables

y_{i} = \tanh (\frac{u z_{i}}{1 - r})

where

r^{2} = \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

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

r_{i} = \frac{h}{4} \sqrt{n + 8 (i - 1)}, i = 1, 2, 3, \dots .

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 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

r_{0} < 1

. The grid-spacing

h

is determined by

r_{0}

and the number of layers to be used. A typical value of

r_{0}

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 Further Comments.

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

Parameters

Compulsory Input Parameters

1: $ndim$ – int64int32nag_int scalar

n

, the number of dimensions of the integral.

Constraint:

1 \leq ndim \leq 30

2: $f$ – function handle or string containing name of m-file

f must return the value of the integrand

f

at a given point.

[result] = f(ndim, x)

Input Parameters

1: $ndim$ – int64int32nag_int scalar: $n$ , the number of dimensions of the integral.
2: $x (ndim)$ – double array: The coordinates of the point at which the integrand $f$ must be evaluated.

Output Parameters

1: $result$ – double scalar: The value of $f (x)$ evaluated at x.

3: $sigma$ – double scalar

Indicates the region of integration.

$sigma \geq 0.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$ – function handle or string containing name of m-file

sigma < 0.0

, region must evaluate the limits of integration in any dimension.

[c, d] = region(ndim, x, j)

Input Parameters

1: $ndim$ – int64int32nag_int scalar: $n$ , the number of dimensions of the integral.
2: $x (ndim)$ – double array: $x (1), \dots, x (j - 1)$ contain the current values of the first $(j - 1)$ variables, which may be used if necessary in calculating $c_{j}$ and $d_{j}$ .
3: $j$ – int64int32nag_int scalar: The index $j$ for which the limits of the range of integration are required.

Output Parameters

1: $c$ – double scalar: The lower limit $c_{j}$ of the range of $x_{j}$ .
2: $d$ – double scalar: The upper limit $d_{j}$ of the range of $x_{j}$ .

sigma \geq 0.0

, region is not called by nag_quad_md_sphere (d01fd), string 'd01fdv'

5: $limit$ – int64int32nag_int scalar

The approximate maximum number of integrand evaluations to be used.

Constraint:

limit \geq 100

Optional Input Parameters

1: $r0$ – double scalar: Suggested value: a typical value is $r0 = 0.8$ . (See also Further Comments.)
Default: $0.8$
The cut-off radius on the unit $n$ -sphere, which may be regarded as an adjustable parameter of the method.

Constraint: $0.0 < r0 < 1.0$ .
2: $u$ – double scalar: Suggested value: a typical value is $u = 1.5$ . (See also Further Comments.)
Default: $1.5$
Must specify an adjustable parameter of the transformation to the unit $n$ -sphere.

Constraint: $u > 0.0$ .

Output Parameters

1: $result$ – double scalar: The approximation to the integral $I$ .
2: $ncalls$ – int64int32nag_int scalar: The actual number of integrand evaluations used. (See also Further Comments.)
3: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$

On entry,	$ndim < 1$ ,
or	$ndim > 30$ .

$ifail = 2$

On entry,

limit < 100

$ifail = 3$

On entry,	$r0 \leq 0.0$ ,
or	$r0 \geq 1.0$ .

$ifail = 4$

On entry,

u \leq 0.0

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

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).

Further Comments

The time taken by nag_quad_md_sphere (d01fd) 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

\frac{u r_{0}}{1 - r_{0}^{2}} > 0.3465 (t - 1),   if ​ sigma \geq 0.0,

\frac{u r_{0}}{1 - r_{0}} > 0.3465 (t - 1),   if ​ sigma < 0.0,

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:

\begin{array}{lr} n = 1 : & 56 \\ n = 2 : & 1252 \\ n = 3 : & 23690 \\ n = 4 : & 394528 \\ n = 5 : & 5956906 \end{array}

Example

This example calculates the integral

\underset{s}{\int \int \int} \frac{d x_{1} d x_{2} d x_{3}}{\sqrt{σ^{2} - r^{2}}} = 22.2066

where

s

is the

3

-sphere of radius

σ

r^{2} = x_{1}^{2} + x_{2}^{2} + x_{3}^{2}

and

σ = 1.5

. Both sphere-to-sphere 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

Open in the MATLAB editor: d01fd_example

function d01fd_example


fprintf('d01fd example results\n\n');

ndim = int64(3);
sigma = 1.5;
limit = int64(8000);
r0 = 0.9;
[result, ncalls, ifail] = ...
  d01fd(...
        ndim, @f, sigma, @region, limit, 'r0', r0);

fprintf('Sphere-to-sphere transformation\n\n');
fprintf('Estimated value of the integral = %9.3f\n', result);
fprintf('Number of integrand evaluations = %5d\n\n', ncalls);

sigma = -1;
r0 = 0.8;
[result, ncalls, ifail] = ...
  d01fd(...
        ndim, @f, sigma, @region, limit, 'r0', r0);
fprintf('Product region transformation\n\n');
fprintf('Estimated value of the integral = %9.3f\n', result);
fprintf('Number of integrand evaluations = %5d\n', ncalls);




function result = f(ndim,x)
  result=1/sqrt(abs(9/4 - dot(x,x)));

function [c,d] = region(ndim, x, j)
  d = sqrt(abs(9/4-dot(x(1:j-1),x(1:j-1))));
  c = -d;

d01fd example results

Sphere-to-sphere transformation

Estimated value of the integral =    22.168
Number of integrand evaluations =  8026

Product region transformation

Estimated value of the integral =    22.137
Number of integrand evaluations =  8026

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox