Real (Kind=nag_wp)	::	d01fbf
Integer, Intent (In)	::	ndim, nptvec(ndim), lwa
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), External	::	fun
Real (Kind=nag_wp), Intent (In)	::	weight(lwa), abscis(lwa)

C Header Interface

#include nagmk26.h

double

d01fbf_ ( const Integer *ndim, const Integer nptvec[], const Integer *lwa, const double weight[], const double abscis[],
double (NAG_CALL *fun)( const Integer *ndim, const double x[]),
Integer *ifail)

3

Description

d01fbf approximates a multidimensional integral by evaluating the summation

\sum_{i_{1} = 1}^{l_{1}} w_{1, i_{1}} \sum_{i_{2} = 1}^{l_{2}} w_{2, i_{2}} \dots \sum_{i_{n} = 1}^{l_{n}} w_{n, i_{n}} f (x_{1, i_{1}}, x_{2, i_{2}}, \dots, x_{n, i_{n}})

given the weights

w_{j, i_{j}}

and abscissae

x_{j, i_{j}}

for a multidimensional product integration rule (see Davis and Rabinowitz (1975)). The number of dimensions may be anything from

1

20

The weights and abscissae for each dimension must have been placed in successive segments of the arrays weight and abscis; for example, by calling d01bcf or d01tbf once for each dimension using a quadrature formula and number of abscissae appropriate to the range of each

x_{j}

and to the functional dependence of

f

x_{j}

If normal weights are used, the summation will approximate the integral

\int w_{1} (x_{1}) \int w_{2} (x_{2}) \dots \int w_{n} (x_{n}) f (x_{1}, x_{2}, \dots, x_{n}) d x_{n} \dots d x_{2} d x_{1}

where

w_{j} (x)

is the weight function associated with the quadrature formula chosen for the

j

th dimension; while if adjusted weights are used, the summation will approximate the integral

\int \int \dots \int f (x_{1}, x_{2}, \dots, x_{n}) d x_{n} \dots d x_{2} d x_{1} .

You must supply a subroutine to evaluate

f (x_{1}, x_{2}, \dots, x_{n})

at any values of

x_{1}, x_{2}, \dots, x_{n}

within the range of integration.

4

References

Davis P J and Rabinowitz P (1975) Methods of Numerical Integration Academic Press

5

Arguments

1: $ndim$ – IntegerInput

On entry:

n

, the number of dimensions of the integral.

Constraint:

1 \leq ndim \leq 20

2: $nptvec (ndim)$ – Integer arrayInput

On entry:

nptvec (j)

must specify the number of points in the

j

th dimension of the summation, for

j = 1, 2, \dots, n

3: $lwa$ – IntegerInput

On entry: the dimension of the arrays weight and abscis as declared in the (sub)program from which d01fbf is called.

Constraint:

lwa \geq nptvec (1) + nptvec (2) + \dots + nptvec (ndim)

4: $weight (lwa)$ – Real (Kind=nag_wp) arrayInput

On entry: must contain in succession the weights for the various dimensions, i.e.,

weight (k)

contains the

i

th weight for the

j

th dimension, with

k = nptvec (1) + nptvec (2) + \dots + nptvec (j - 1) + i .

5: $abscis (lwa)$ – Real (Kind=nag_wp) arrayInput

On entry: must contain in succession the abscissae for the various dimensions, i.e.,

abscis (k)

contains the

i

th abscissa for the

j

th dimension, with

k = nptvec (1) + nptvec (2) + \dots + nptvec (j - 1) + i .

6: $fun$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

fun must return the value of the integrand

f

at a specified point.

The specification of fun is:

Fortran Interface

Function fun (

ndim, x)

Real (Kind=nag_wp)	::	fun
Integer, Intent (In)	::	ndim
Real (Kind=nag_wp), Intent (In)	::	x(ndim)

C Header Interface

#include nagmk26.h

double

fun ( const Integer *ndim, const double x[])

1: $ndim$ – IntegerInput: On entry: $n$ , the number of dimensions of the integral.
2: $x (ndim)$ – Real (Kind=nag_wp) arrayInput: On entry: the coordinates of the point at which the integrand $f$ must be evaluated.

fun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01fbf is called. Arguments denoted as Input must not be changed by this procedure.

Note: fun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01fbf. If your code inadvertently does return any NaNs or infinities, d01fbf is likely to produce unexpected results.

7: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 ​ or ​ 1

. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 ​ or ​ 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, if you are not familiar with this argument, the recommended value is

0

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

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:

$ifail = 1$

On entry,	$ndim < 1$ ,
or	$ndim > 20$ ,
or	$lwa < nptvec (1) + nptvec (2) + \dots + nptvec (ndim)$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

The accuracy of the computed multidimensional sum depends on the weights and the integrand values at the abscissae. If these numbers vary significantly in size and sign then considerable accuracy could be lost. If these numbers are all positive, then little accuracy will be lost in computing the sum.

8

Parallelism and Performance

d01fbf is not threaded in any implementation.

9

Further Comments

The total time taken by d01fbf will be proportional to

T \times nptvec (1) \times nptvec (2) \times \dots \times nptvec (ndim),

where

T

is the time taken for one evaluation of fun.

10

Example

This example evaluates the integral

\int_{1}^{2} \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{1}^{\infty} \frac{{(x_{1} x_{2} x_{3})}^{6}}{{(x_{4} + 2)}^{8}} e^{- 2 x_{2}} e^{- 0.5 x_{3}^{2}} d x_{4} d x_{3} d x_{2} d x_{1}

using adjusted weights. The quadrature formulae chosen are:

$x_{1}$ : Gauss–Legendre, $a = 1.0$ , $b = 2.0$ ,
$x_{2}$ : Gauss–Laguerre, $a = 0.0$ , $b = 2.0$ ,
$x_{3}$ : Gauss–Hermite, $a = 0.0$ , $b = 0.5$ ,
$x_{4}$ : rational Gauss, $a = 1.0$ , $b = 2.0$ .

Four points are sufficient in each dimension, as this integral is in fact a product of four one-dimensional integrals, for each of which the chosen four-point formula is exact.