next up previous index
Next: D104 Cauchy Principal Up: CERNLIB Previous: D102 Adaptive Gaussian

D103 Adaptive Gaussian Quadrature

Routine ID: D103
Author(s): K.S. KölbigLibrary: MATHLIB
Submitter: Submitted: 02.05.1966
Language: FortranRevised: 15.03.1993

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


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


FUNCTION subprograms
Files Referenced: Unit 6
External References: MTLMTR (N002), ABEND (Z035), user-supplied FUNCTION subprogram


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.

Name of a user-supplied FUNCTION subprogram, declared EXTERNAL in the calling program. This subprogram must set F(X)= f(X) .
End-points of integration interval. Note that B may be less than A.
Accuracy parameter (see Accuracy).
GAUSS is of type REAL, DGAUSS is of type DOUBLE PRECISION, QGAUSS is of type REAL*16, and the arguments F, A, B, EPS and X (in F) have the same type as the function name.


For any interval [a,b] we define g8(a,b) and g16(a,b) to be the 8-point and 16-point Gaussian quadrature approximations to


and define

r(a,b) ={|g16(a,b) - g8(a,b)|1+|g16(a,b)|}.

Then, with G = GAUSS or DGAUSS,

G =∑i=1kg16(xi-1,xi),

where, starting with x0=A and finishing with xk=B , the subdivision points xi (i=1,2,...) are given by

xi= xi-1+ λ(B-xi-1),

with λ equal to the first member of the sequence 1,{12},{14},... for which r(xi-1,xi) < EPS . If, at any stage in the process of subdivision, the ratio


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.


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

Iabs= ∫AB|f(x)|dx,

then the relation

{|G - I|Iabs+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.


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.


next up previous index
Next: D104 Cauchy Principal Up: CERNLIB Previous: D102 Adaptive Gaussian

Janne Saarela
Mon Apr 3 15:06:23 METDST 1995