D103 Adaptive Gaussian Quadrature

Function subprograms GAUSS, DGAUSS and QGAUSS compute, to an attempted specified accuracy, the value of the integral

$I=\int$_A^Bf(x)dx.

The quadruple-precision version QGAUSS is available only on computers which support a REAL*16 Fortran data type.

Structure:

FUNCTION subprograms
User Entry Names: GAUSS, DGAUSS, QGAUSS
Files Referenced: Unit 6
External References: MTLMTR (N002), ABEND (Z035), user-supplied FUNCTION subprogram

Usage:

In any arithmetic expression, GAUSS(F,A,B,EPS), DGAUSS(F,A,B,EPS) or QGAUSS(F,A,B,EPS)

has the approximate value of the integral I.

F

Name of a user-supplied FUNCTION subprogram, declared EXTERNAL in the calling program. This subprogram must set $F(X)=\; f(X)$ .

A,B

End-points of integration interval. Note that B may be less than A.

EPS

Accuracy parameter (see Accuracy).

Method:

For any interval $[a,b]$ we define $g$₈(a,b) and $g$₁₆(a,b) to be the 8-point and 16-point Gaussian quadrature approximations to

$\int$_a^bf(x)dx

and define

$r(a,b)\; =\{|g$₁₆(a,b) - g₈(a,b)|1+|g₁₆(a,b)|}.

Then, with G = GAUSS or DGAUSS,

$G\; =\sum$_i=1^kg₁₆(x_i-1,x_i),

where, starting with $x$₀=A and finishing with $x$_k=B , the subdivision points $x$_i (i=1,2,...) are given by

$x$_i= x_i-1+ λ(B-x_i-1),

with $\lambda$ equal to the first member of the sequence $1,\{12\},\{14\},...$ for which $r(x$_i-1,x_i) < EPS . If, at any stage in the process of subdivision, the ratio

$q=|\{x$_i-x_i-1B-A}|

is so small that $1+0.005q$ is indistinguishable from 1 to machine accuracy, an error exit occurs with the function value set equal to zero.

Accuracy:

Unless there is severe cancellation of positive and negative values of $f(x)$ over the interval $[A,B]$ , the argument EPS may be considered as specifying a bound on the relative error of I in the case |I|>1, and a bound on the absolute error in the case |I|<1. More precisely, if k is the number of sub-intervals contributing to the approximation (see Method), and if

$I$_abs= ∫_A^B|f(x)|dx,

then the relation

$\{|G\; -\; I|I$_abs+k}< EPS

will nearly always be true, provided the routine terminates without printing an error message. For functions f having no singularities in the closed interval $[A,B]$

the accuracy will usually be much higher than this.

Error handling:

Error D103.1: The requested accuracy cannot be obtained (see Method). The function value is set equal to zero, and a message is written on Unit 6 unless subroutine MTLSET (N002) has been called.

Notes:

Values of the function $f(x)$ at the interval end-points A and B are not required. The subprogram may therefore be used when these values are undefined.
$•$

D104