NAG FL Interface
d01uaf (dim1_gauss_vec)

d01uaf computes an estimate of the definite integral of a function of known analytical form, using a Gaussian quadrature formula with a specified number of abscissae. Formulae are provided for a finite interval (Gauss–Legendre), a semi-infinite interval (Gauss–Laguerre, rational Gauss), and an infinite interval (Gauss–Hermite).

2 Specification

Fortran Interface

Subroutine d01uaf (

key, a, b, n, f, dinest, iuser, ruser, ifail)

Integer, Intent (In)	::	key, n
Integer, Intent (Inout)	::	iuser(*), ifail
Real (Kind=nag_wp), Intent (In)	::	a, b
Real (Kind=nag_wp), Intent (Inout)	::	ruser(*)
Real (Kind=nag_wp), Intent (Out)	::	dinest
External	::	f

C Header Interface

#include <nag.h>

void	d01uaf_ (const Integer key, const double a, const double b, const Integer n, void (NAG_CALL f)(const double x[], const Integer nx, double fv[], Integer iflag, Integer iuser[], double ruser[]), double dinest, Integer iuser[], double ruser[], Integer *ifail)

The routine may be called by the names d01uaf or nagf_quad_dim1_gauss_vec.

3 Description

3.1 General

d01uaf evaluates an estimate of the definite integral of a function

f (x)

, over a finite or infinite range, by

n

-point Gaussian quadrature (see Davis and Rabinowitz (1975), Fröberg (1970), Ralston (1965) or Stroud and Secrest (1966)). The integral is approximated by a summation of the product of a set of weights and a set of function evaluations at a corresponding set of abscissae

x_{i}

. For adjusted weights, the function values correspond to the values of the integrand

f

, and hence the sum will be

\sum_{i = 1}^{n} w_{i} f (x_{i})

where the

w_{i}

are called the weights, and the

x_{i}

the abscissae. A selection of values of

n

is available. (See Section 5.)

Where applicable, normal weights may instead be used, in which case the corresponding weight function

ω

is factored out of the integrand as

f (x) = ω (x) g (x)

and hence the sum will be

\sum_{i = - 1}^{n} {\bar{w}}_{i} g (x),

where the normal weights

{\bar{w}}_{i} = w_{i} ω (x_{i})

are computed internally.

d01uaf uses a vectorized f to evaluate the integrand or normalized integrand at a set of abscissae,

x_{i}

, for

i = 1, 2, \dots, n_{x}

. If adjusted weights are used, the integrand

f (x_{i})

must be evaluated otherwise the normalized integrand

g (x_{i})

must be evaluated.

3.2 Both Limits Finite

\int_{a}^{b} f (x) d x .

The Gauss–Legendre weights and abscissae are used, and the formula is exact for any function of the form:

f (x) = \sum_{i = 0}^{2 n - 1} c_{i} x^{i} .

The formula is appropriate for functions which can be well approximated by such a polynomial over

[a, b]

. It is inappropriate for functions with algebraic singularities at one or both ends of the interval, such as

{(1 + x)}^{- 1 / 2}

[−1, 1]

3.3 One Limit Infinite

\int_{a}^{\infty} f (x) d x or \int_{- \infty}^{a} f (x) d x .

Two quadrature formulae are available for these integrals.

(a)The Gauss–Laguerre formula is exact for any function of the form:

$f (x) = e^{- b x} \sum_{i = 0}^{2 n - 1} c_{i} x^{i} .$

This formula is appropriate for functions decaying exponentially at infinity; the argument $b$ should be chosen if possible to match the decay rate of the function.

If the adjusted weights are selected, the complete integrand $f (x)$ should be provided through f.

If the normal form is selected, the contribution of $e^{- b x}$ is accounted for internally, and f should only return $g (x)$ , where $f (x) = e^{- b x} g (x)$ .

If $b < 0$ is supplied, the interval of integration will be $[a, \infty)$ . Otherwise if $b > 0$ is supplied, the interval of integration will be $(- \infty, a]$ .
(b)The rational Gauss formula is exact for any function of the form:

$f (x) = \sum_{i = 2}^{2 n + 1} \frac{c_{i}}{{(x + b)}^{i}} = \frac{\sum_{i = 0}^{2 n - 1} c_{2 n + 1 - i} {(x + b)}^{i}}{{(x + b)}^{2 n + 1}} .$

This formula is likely to be more accurate for functions having only an inverse power rate of decay for large $x$ . Here the choice of a suitable value of $b$ may be more difficult; unfortunately a poor choice of $b$ can make a large difference to the accuracy of the computed integral.

Only the adjusted form of the rational Gauss formula is available, and as such, the complete integrand $f (x)$ must be supplied in f.

If $a + b < 0$ , the interval of integration will be $[a, \infty)$ . Otherwise if $a + b > 0$ , the interval of integration will be $(- \infty, a]$ .

3.4 Both Limits Infinite

\int_{- \infty}^{+ \infty} f (x) d x .

The Gauss–Hermite weights and abscissae are used, and the formula is exact for any function of the form:

f (x) = e^{- b {(x - a)}^{2}} \sum_{i = 0}^{2 n - 1} c_{i} x^{i},

where

b > 0

. Again, for general functions not of this exact form, the argument

b

should be chosen to match if possible the decay rate at

\pm \infty

If the adjusted weights are selected, the complete integrand

f (x)

should be provided through f.

If the normal form is selected, the contribution of

e^{- b {(x - a)}^{2}}

is accounted for internally, and f should only return

g (x)

, where

f (x) = e^{- b {(x - a)}^{2}} g (x)

4 References

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

Fröberg C E (1970) Introduction to Numerical Analysis Addison–Wesley

Ralston A (1965) A First Course in Numerical Analysis pp. 87–90 McGraw–Hill

Stroud A H and Secrest D (1966) Gaussian Quadrature Formulas Prentice–Hall

5 Arguments

1: $key$ – Integer Input

On entry: indicates the quadrature formula.

$key = 0$: Gauss–Legendre quadrature on a finite interval, using normal weights.
$key = 3$: Gauss–Laguerre quadrature on a semi-infinite interval, using normal weights.
$key = −3$: Gauss–Laguerre quadrature on a semi-infinite interval, using adjusted weights.
$key = 4$: Gauss–Hermite quadrature on an infinite interval, using normal weights.
$key = −4$: Gauss–Hermite quadrature on an infinite interval, using adjusted weights.
$key = −5$: Rational Gauss quadrature on a semi-infinite interval, using adjusted weights.

Constraint:

key = 0

3

−3

4

−4

−5

2: $a$ – Real (Kind=nag_wp) Input

3: $b$ – Real (Kind=nag_wp) Input

On entry: the quantities

a

and

b

as described in the appropriate subsection of Section 3.

Constraints:

Rational Gauss: $a + b \neq 0.0$ ;
Gauss–Laguerre: $b \neq 0.0$ ;
Gauss–Hermite: $b > 0$ .

4: $n$ – Integer Input

On entry:

n

, the number of abscissae to be used.

Constraint:

n = 1

2

3

4

5

6

8

10

12

14

16

20

24

32

48

64

If the soft fail option is used, the answer is evaluated for the largest valid value of n less than the requested value.

5: $f$ – Subroutine, supplied by the user. External Procedure

f must return the value of the integrand

f

, or the normalized integrand

g

, at a set of points.

The specification of f is:

Fortran Interface

Subroutine f (

x, nx, fv, iflag, iuser, ruser)

Integer, Intent (In)	::	nx
Integer, Intent (Inout)	::	iflag, iuser(*)
Real (Kind=nag_wp), Intent (In)	::	x(nx)
Real (Kind=nag_wp), Intent (Inout)	::	ruser(*)
Real (Kind=nag_wp), Intent (Out)	::	fv(nx)

C Header Interface

void	f (const double x[], const Integer nx, double fv[], Integer iflag, Integer iuser[], double ruser[])

1: $x (nx)$ – Real (Kind=nag_wp) array Input: On entry: the abscissae, $x_{i}$ , for $i = 1, 2, \dots, n_{x}$ at which function values are required.
2: $nx$ – Integer Input: On entry: $n_{x}$ , the number of abscissae.
3: $fv (nx)$ – Real (Kind=nag_wp) array Output: On exit: if adjusted weights are used, the values of the integrand $f$ . $fv (i) = f (x_{i})$ , for $i = 1, 2, \dots, n_{x}$ .
Otherwise the values of the normalized integrand $g$ . $fv (i) = g (x_{i})$ , for $i = 1, 2, \dots, n_{x}$ .
4: $iflag$ – Integer Input/Output: On entry: $iflag = 0$ .

On exit: set $iflag < 0$ if you wish to force an immediate exit from d01uaf with $ifail = - 1$ .
5: $iuser (*)$ – Integer array User Workspace
6: $ruser (*)$ – Real (Kind=nag_wp) array User Workspace: f is called with the arguments iuser and ruser as supplied to d01uaf. You should use the arrays iuser and ruser to supply information to f.

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

Some points to bear in mind when coding f are mentioned in Section 7.

6: $dinest$ – Real (Kind=nag_wp) Output

On exit: the estimate of the definite integral.

7: $iuser (*)$ – Integer array User Workspace

8: $ruser (*)$ – Real (Kind=nag_wp) array User Workspace

iuser and ruser are not used by d01uaf, but are passed directly to f and may be used to pass information to this routine.

9: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

−1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. 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:

Note: in some cases d01uaf may return useful information.

$ifail = 1$: The $n$ -point rule is not among those stored.
On entry: $n = ⟨ value ⟩$ .
$n$ -point rule used: $n = ⟨ value ⟩$ .

$ifail = 2$: Underflow occurred in calculation of normal weights.
Reduce n or use adjusted weights: $n = ⟨ value ⟩$ .

$ifail = 3$: No nonzero weights were generated for the provided parameters.

$ifail = 11$: On entry, $key = ⟨ value ⟩$ .
Constraint: $key = 0$ , $3$ , $−3$ , $4$ , $−4$ or $−5$ .

$ifail = 12$

The value of a and/or b is invalid for the chosen key. Either:

The value of a and/or b is invalid.
On entry, $key = ⟨ value ⟩$ .
On entry, $a = ⟨ value ⟩$ and $b = ⟨ value ⟩$ .
Constraint: $| a + b | > 0.0$ .
The value of a and/or b is invalid.
On entry, $key = ⟨ value ⟩$ .
On entry, $a = ⟨ value ⟩$ and $b = ⟨ value ⟩$ .
Constraint: $| b | > 0.0$ .
The value of a and/or b is invalid.
On entry, $key = ⟨ value ⟩$ .
On entry, $a = ⟨ value ⟩$ and $b = ⟨ value ⟩$ .
Constraint: $b > 0.0$ .

$ifail = 14$: On entry, $n = ⟨ value ⟩$ .
Constraint: $n > 0$ .

$ifail = - 1$: Exit requested from f with $iflag = ⟨ value ⟩$ .

$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

The accuracy depends on the behaviour of the integrand, and on the number of abscissae used. No tests are carried out in d01uaf to estimate the accuracy of the result. If such an estimate is required, the routine may be called more than once, with a different number of abscissae each time, and the answers compared. It is to be expected that for sufficiently smooth functions a larger number of abscissae will give improved accuracy.

Alternatively, the range of integration may be subdivided, the integral estimated separately for each sub-interval, and the sum of these estimates compared with the estimate over the whole range.

The coding of f may also have a bearing on the accuracy. For example, if a high-order Gauss–Laguerre formula is used, and the integrand is of the form

f (x) = e^{- b x} g (x)

it is possible that the exponential term may underflow for some large abscissae. Depending on the machine, this may produce an error, or simply be assumed to be zero. In any case, it would be better to evaluate the expression with

f (x) = sgn (g (x)) \times \exp (- b x + \ln | g (x) |)

Another situation requiring care is exemplified by

\int_{- \infty}^{+ \infty} e^{- x^{2}} x^{m} d x = 0, m ​ odd.

The integrand here assumes very large values; for example, when

m = 63

, the peak value exceeds

3 \times 10^{33}

. Now, if the machine holds floating-point numbers to an accuracy of

k

significant decimal digits, we could not expect such terms to cancel in the summation leaving an answer of much less than

10^{33 - k}

(the weights being of order unity); that is, instead of zero we obtain a rather large answer through rounding error. Such situations are characterised by great variability in the answers returned by formulae with different values of

n

In general, you should be aware of the order of magnitude of the integrand, and should judge the answer in that light.

8 Parallelism and Performance

d01uaf is currently neither directly nor indirectly threaded. In particular, the user-supplied argument f is not called from within a parallel region initialized inside d01uaf.

The user-supplied argument f uses a vectorized interface, allowing for the required vector of function values to be evaluated in parallel; for example by placing appropriate OpenMP directives in the code for the user-supplied argument f.

9 Further Comments

The time taken by d01uaf depends on the complexity of the expression for the integrand and on the number of abscissae required.

10 Example

This example evaluates the integrals

\int_{0}^{1} \frac{4}{1 + x^{2}} d x = π

by Gauss–Legendre quadrature;

\int_{2}^{\infty} \frac{1}{x^{2} \ln x} d x = 0.378671

by rational Gauss quadrature with

b = 0

;

\int_{2}^{\infty} \frac{e^{- x}}{x} d x = 0.048901

by Gauss–Laguerre quadrature with

b = 1

; and

\int_{- \infty}^{+ \infty} e^{- 3 x^{2} - 4 x - 1} d x = \int_{- \infty}^{+ \infty} e^{- 3 {(x + 1)}^{2}} e^{2 x + 2} d x = 1.428167

by Gauss–Hermite quadrature with

a = −1

and

b = 3

The formulae with

n = 2, 4, 8, 16, 32

and

64

are used in each case. Both adjusted and normal weights are used for Gauss–Laguerre and Gauss–Hermite quadrature.

10.1 Program Text

Program Text (d01uafe.f90)

10.2 Program Data

None.

10.3 Program Results

Program Results (d01uafe.r)

NAG Library Manual, Mark 28.4

Interfaces: FL CL CPP AD PY MB

NAG FL Interface Introduction

D01 (Quad) Chapter Contents

D01 (Quad) Chapter Introduction

d01ua: FL CL CPP AD PY MB

NAG FL Interfaced01uaf (dim1_​gauss_​vec)

▸▿ Contents

1 Purpose