hide long namesshow long names
hide short namesshow short names
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)

Purpose

nag_quad_md_sphere (d01fd) calculates an approximation to a definite integral in up to 3030 dimensions, using the method of Sag and Szekeres (see Sag and Szekeres (1964)). The region of integration is an nn-sphere, or by built-in transformation via the unit nn-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
 f(x1,x2,,xn)dx1dx2dxn
n-sphere of radius  σ
n-sphere of radius  σ f (x1,x2,,xn) dx1 dx2 dxn
(1)
or, more generally,
d1 dn
dx1dxnf(x1,,xn)
c1 cn
c1 d1 dx1 cn dn dxn f (x1,,xn)
(2)
where each cici and didi may be functions of xjxj (j < i)(j<i).
The function uses the method of Sag and Szekeres (1964), which exploits a property of the shifted pp-point trapezoidal rule, namely, that it integrates exactly all polynomials of degree < p<p (see Krylov (1962)). An attempt is made to induce periodicity in the integrand by making a parameterised transformation to the unit nn-sphere. The Jacobian of the transformation and all its direct derivatives vanish rapidly towards the surface of the unit nn-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 nn-sphere by the change of variables:
xi = yi σ/r tanh((ur)/(1r2))
xi = yi σr tanh( ur 1-r2 )
where r2 = i = 1nyi2r2=i=1nyi2 and uu is an adjustable parameter.
Integrals of the form (2) are first of all transformed to the nn-cube [1,1]n[-1,1]n by a linear change of variables
xi = ((di + ci) + (dici)yi) / 2
xi=((di+ci)+(di-ci)yi)/2
and then to the unit sphere by a further change of variables
yi = tanh((uzi)/(1r))
yi=tanh(uzi 1-r )
where r2 = i = 1nzi2r2=i=1nzi2 and uu is again an adjustable parameter.
The parameter uu in these transformations determines how the transformed integrand is distributed between the origin and the surface of the unit nn-sphere. A typical value of uu is 1.51.5. For larger uu, the integrand is concentrated toward the centre of the unit nn-sphere, while for smaller uu it is concentrated toward the perimeter.
In performing the integration over the unit nn-sphere by the trapezoidal rule, a displaced equidistant grid of size hh is constructed. The points of the mesh lie on concentric layers of radius
ri = h/4sqrt(n + 8(i1)),  i = 1,2,3,.
ri=h4n+8(i-1),  i=1,2,3,.
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 400400 layers.
In practice, the rapidly-decreasing Jacobian makes it unnecessary to include the whole unit nn-sphere and the integration region is limited by a user-specified cut-off radius r0 < 1r0<1. The grid-spacing hh is determined by r0r0 and the number of layers to be used. A typical value of r0r0 is 0.80.8.
Some experimentation may be required with the choice of r0r0 (which determines how much of the unit nn-sphere is included) and uu (which determines how the transformed integrand is distributed between the origin and surface of the unit nn-sphere), to obtain best results for particular families of integrals. This matter is discussed further in Section [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
nn, the number of dimensions of the integral.
Constraint: 1ndim301ndim30.
2:     f – function handle or string containing name of m-file
f must return the value of the integrand ff at a given point.
[result] = f(ndim, x)

Input Parameters

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

Output Parameters

1:     result – double scalar
The result of the function.
3:     sigma – double scalar
Indicates the region of integration.
sigma0.0sigma0.0
The integration is carried out over the nn-sphere of radius sigma, centred at the origin.
sigma < 0.0sigma<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
If sigma < 0.0sigma<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
nn, the number of dimensions of the integral.
2:     x(ndim) – double array
x(1),,x(j1)x1,,xj-1 contain the current values of the first (j1)(j-1) variables, which may be used if necessary in calculating cjcj and djdj.
3:     j – int64int32nag_int scalar
The index jj for which the limits of the range of integration are required.

Output Parameters

1:     c – double scalar
The lower limit cjcj of the range of xjxj.
2:     d – double scalar
The upper limit djdj of the range of xjxj.
If sigma0.0sigma0.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: limit100limit100.

Optional Input Parameters

1:     r0 – double scalar
The cut-off radius on the unit nn-sphere, which may be regarded as an adjustable parameter of the method.
Default: 0.80.8 
Constraint: 0.0 < r0 < 1.00.0<r0<1.0.
2:     u – double scalar
Must specify an adjustable parameter of the transformation to the unit nn-sphere.
Default: 1.51.5 
Constraint: u > 0.0u>0.0.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

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

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,ndim < 1ndim<1,
orndim > 30ndim>30.
  ifail = 2ifail=2
On entry,limit < 100limit<100.
  ifail = 3ifail=3
On entry,r00.0r00.0,
orr01.0r01.0.
  ifail = 4ifail=4
On entry,u0.0u0.0.

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 r0r0 and uu 
If the chosen combination of r0r0 and uu 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 nn-sphere which lie too close to the boundary surface to be distinguished from it to machine accuracy (despite the fact that r0 < 1r0<1). To be specific, the combination of r0r0 and uu is too large if
(ur0)/(1r02) > 0.3465(t1),   if ​sigma0.0,
ur0 1-r02 >0.3465(t-1),   if ​sigma0.0,
or
(ur0)/(1r0) > 0.3465(t1),   if ​ sigma < 0.0,
ur0 1-r0 > 0.3465(t- 1),   if ​ sigma< 0.0,
where tt 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 100100. 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 r0r0 and uu may result in some of the points being unusable. Such points are not included in the returned value of ncalls.
Secondly, no more than 400400 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
n=1: 56 n=2: 1252 n=3: 23690 n=4: 394528 n=5: 5956906

Example

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

function result = f(ndim,x)
  result = 2.25;
  for i = 1:double(ndim)
    result = result - x(i)^2;
  end
  result=1/sqrt(abs(result));
function [c,d] = region(ndim, x, j)
      c = -1.5;
      d = 1.5;
      if (j > 1)
         sm = 2.25;
         for i = 1:double(j-1)
            sm = sm - x(i)*x(i);
         end
         d = sqrt(abs(sm));
         c = -d;
      end
 

result =

   22.1679


ncalls =

                 8026


ifail =

                    0


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

function result = f(ndim,x)
  result = 2.25;
  for i = 1:double(ndim)
    result = result - x(i)^2;
  end
  result=1/sqrt(abs(result));
function [c,d] = region(ndim, x, j)
      c = -1.5;
      d = 1.5;
      if (j > 1)
         sm = 2.25;
         for i = 1:double(j-1)
            sm = sm - x(i)*x(i);
         end
         d = sqrt(abs(sm));
         c = -d;
      end
 

result =

   22.1679


ncalls =

                 8026


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013