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_bad (d01ja)

Purpose

nag_quad_md_sphere_bad (d01ja) attempts to evaluate an integral over an nn-dimensional sphere (n = 2n=2, 33, or 44), to a user-specified absolute or relative accuracy, by means of a modified Sag–Szekeres method. The function can handle singularities on the surface or at the centre of the sphere, and returns an error estimate.

Syntax

[result, esterr, nevals, ifail] = d01ja(f, ndim, radius, epsa, epsr, icoord, 'method', method)
[result, esterr, nevals, ifail] = nag_quad_md_sphere_bad(f, ndim, radius, epsa, epsr, icoord, 'method', method)

Description

nag_quad_md_sphere_bad (d01ja) calculates an approximation to the nn-dimensional integral
I = SF(x1,,xn)dx1dxn,  2n4,
I=SF(x1,,xn)dx1dxn,  2n4,
where SS is the hypersphere
sqrt((x12 + + xn2))α <
(x12++xn2)α<
(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 function 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 hh is larger than the tolerated error, the trapezoidal formula with mesh size h / 2h/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 2n2n). This introduces the restriction n4n4. 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).

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

Parameters

Compulsory Input Parameters

1:     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. These coordinates are given in Cartesian or spherical polar form according to the value of icoord.

Output Parameters

1:     result – double scalar
The result of the function.
2:     ndim – int64int32nag_int scalar
nn, the dimension of the sphere.
Constraint: 2ndim42ndim4.
3:     radius – double scalar
αα, the radius of the sphere.
Constraint: radius0.0radius0.0.
4:     epsa – double scalar
The requested absolute tolerance. If epsa < 0.0epsa<0.0, its absolute value is used. See Section [Accuracy].
5:     epsr – double scalar
The requested relative tolerance.
epsr < 0.0epsr<0.0
Its absolute value is used.
epsr < 10 × (machine precision)epsr<10×(machine precision)
The latter value is used as epsr by the function. See Section [Accuracy].
6:     icoord – int64int32nag_int scalar
Must specify which kind of coordinates are used in f.
icoord = 0icoord=0
Cartesian coordinates xixi, for i = 1,2,,ni=1,2,,n.
icoord = 1icoord=1
Spherical coordinates (see Section [Spherical Polar Coordinates]): x(1) = ρx1=ρ; x(i) = θi1xi=θi-1, for i = 2,3,,ni=2,3,,n.
icoord = 2icoord=2,
Special spherical polar coordinates (see Section [Machine Dependencies]), with the additional transformation ρ = αλρ=α-λ: x(1) = λ = αρx1=λ=α-ρ; x(i) = θi1xi=θi-1, for i = 2,3,,ni=2,3,,n.
Constraint: icoord = 0icoord=0, 11 or 22.
If method = 3method=3 or 44, icoord = 2icoord=2

Optional Input Parameters

1:     method – int64int32nag_int scalar
Must specify the transformation to be used by the function. The choice depends on the behaviour of the integrand and on the required accuracy.
For well-behaved functions and functions with mild singularities on the surface of the sphere only:
method = 1method=1
Low accuracy required.
method = 2method=2
High accuracy required.
For functions with severe singularities on the surface of the sphere only:
method = 3method=3
Low accuracy required.
method = 4method=4
High accuracy required.
(in this case icoord must be set to icoord = 2icoord=2, and the function defined in special spherical coordinates).
For functions with a singularity at the centre of the sphere (and possibly with singularities on the surface as well):
method = 5method=5
Low accuracy required.
method = 6method=6
High accuracy required.
method = 0method=0 can be used as a default value and is equivalent to:
  • method = 1method=1 if epsr > 106epsr>10-6, and
  • method = 2method=2 if epsr106epsr10-6.
The distinction between low and high required accuracies, as mentioned above, depends also on the behaviour of the function. Roughly one may assume the critical value of epsa and epsr to be 10610-6, but the critical value will be smaller for a well-behaved integrand and larger for an integrand with severe singularities.
Default: 00 
Constraint: method = 0method=0, 11, 22, 33, 44, 55 or 66.
If icoord = 2icoord=2, method = 3method=3 or 44

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     result – double scalar
The approximation to the integral II.
2:     esterr – double scalar
An estimate of the modulus of the absolute error.
3:     nevals – int64int32nag_int scalar
The number of function evaluations used.
4:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_quad_md_sphere_bad (d01ja) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W ifail = 1ifail=1
The required accuracy cannot be achieved within a limiting number of function evaluations (which is set by the function).
W ifail = 2ifail=2
The required accuracy cannot be achieved because of round-off error.
W ifail = 3ifail=3
The required accuracy cannot be achieved because the maximum accuracy with respect to the machine constants nag_machine_precision (x02aj) and nag_machine_real_safe (x02am) has been attained. If this maximum accuracy is rather low (compared with nag_machine_precision (x02aj)), the cause of the problem is a severe singularity on the boundary or at the centre of the sphere. If method = 0method=0, 11 or 22, then setting method = 3method=3 or 44 may help.
  ifail = 4ifail=4
On entry,ndim < 2ndim<2 or ndim > 4ndim>4,
orradius < 0.0radius<0.0,
ormethod0method0, 11, 22, 33, 44, 55 or 66,
oricoord0icoord0, 11 or 22,
oricoord = 2icoord=2 and method3method3 or 44,
ormethod = 3method=3 or 44 and icoord2icoord2.
No calculations have been performed. result and esterr are set to 0.00.0.

Accuracy

You can specify an absolute and/or a relative tolerance, setting epsa and epsr. The function attempts to calculate an approximation result such that
|Iresult|max {epsa,epsr × |I|}.
|I-result|max{epsa,epsr×|I|}.
If 0ifail30ifail3, esterr returns an estimate of, but not necessarily a bound for, |Iresult||I-result|.

Further Comments

Timing

Timing depends on the integrand and the accuracy required.

Spherical Polar Coordinates

Cartesian coordinates are related to the spherical polar coordinates by:
x1 = ρ . sinθ1sinθn2 . sinθn1
x2 = ρ . sinθ1sinθn2 . cosθn1
x3 = ρ . sinθ1cosθn2
xn = ρ . cosθ1
x1 = ρ.sinθ1sinθn-2.sinθn-1 x2 = ρ.sinθ1sinθn-2.cosθn-1 x3 = ρ.sinθ1cosθn-2 xn = ρ.cosθ1
where 0 < θi < π0<θi<π, for i = 1,2,,n2i=1,2,,n-2 and 0 < θn1 < 2π0<θn-1<2π.

Machine Dependencies

As a consequence of the transformation technique, the severity of the singularities which can be handled by nag_quad_md_sphere_bad (d01ja) depends on the precision and range of real numbers on the machine. method = 3method=3 or 44 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 = 3method=3 or 44 if nag_quad_md_sphere_bad (d01ja) terminates with ifail = 3ifail=3 when called with method = 0method=0, 11 or 22.
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 nag_quad_md_sphere_bad (d01ja) may be very close to the surface or to the centre of the sphere.
Example:
Note that nag_quad_md_sphere_bad (d01ja) ensures that λ = x(1) > x02amλ=x1>x02am, but underflow could occur in the computation of λ2λ2.

Example

function nag_quad_md_sphere_bad_example
ndim = int64(2);
radius = 1;
epsa = 0;
epsr = 5e-05;
icoord = int64(1);
[result, esterr, nevals, ifail] = ...
    nag_quad_md_sphere_bad(@f, ndim, radius, epsa, epsr, icoord)

function result = f(ndim, x)
  a = (1-x(1))*(1+x(1));
  if (a == 0)
    result = 0;
  else
    result = 1/sqrt(a);
  end
 

result =

    6.2832


esterr =

   1.9784e-04


nevals =

                  193


ifail =

                    0


function d01ja_example
ndim = int64(2);
radius = 1;
epsa = 0;
epsr = 5e-05;
icoord = int64(1);
[result, esterr, nevals, ifail] = ...
    d01ja(@f, ndim, radius, epsa, epsr, icoord)

function result = f(ndim, x)
  a = (1-x(1))*(1+x(1));
  if (a == 0)
    result = 0;
  else
    result = 1/sqrt(a);
  end
 

result =

    6.2832


esterr =

   1.9784e-04


nevals =

                  193


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