D110: Gaussian Quadrature for Multiple Integrals

Author(s): K.S. Kölbig	Library: MATHLIB
Submitter:	Submitted: 01.12.1988
Language: Fortran	Revised: 15.03.1993

Function subprogram packages RGMLT and DGMLT compute an approximate value of an n-dimensional integral of the form


where .

Each subprogram integrates over only one variable. The integral is computed by means of a set of n nested calls to these subprograms.

On computers other than CDC or Cray, only the double-precision version DGMLT is available. On CDC and Cray computers, only the single-precision version RGMLT is available.

Structure:

FUNCTION subprograms
User Entry Names:
Files Referenced: Unit 6
External References: MTLMTR, ABEND, user-supplied SUBROUTINE subprograms

Usage:

1. Let k be one of the integers . Then, in any arithmetic expression,

   RGMLTk(FSUBk,Ak,Bk,NIk,NGk,X) or

   DGMLTk(FSUBk,Ak,Bk,NIk,NGk,X)

   has the approximate value of the integral with respect to of the function specified below.
   RGMLTk is of type REAL, DGMLTk is of type DOUBLE PRECISION, and the arguments Ak, Bk, and X have the same type as the function name. NIk and NGk are of type INTEGER.

   FSUBk

   Name of a user-supplied SUBROUTINE subprogram, declared EXTERNAL in the calling program.

   Ak,Bk

   End points of the integration interval for variable .

   NIk

   Number of equal subintervals into which the interval (Ak,Bk) is divided.

   NGk

   Number of Gaussian quadrature points to be used in each of the NIk subintervals. (If NGk has any value other than 6 or 8, the value 6 is assumed).

   X

   0ne-dimensional array of dimension .

2. The subroutine FSUBk should be of the form

   SUBROUTINE FSUBk (M,Uk,Fk,X)

   where Uk, Fk and X are of type REAL for RGMLTk and of type DOUBLE PRECISION for DGMLTk, and where M is of type INTEGER.

   M

   An integer , whose value is set by RGMLTk (DGMLTk).

   Uk

   One-dimensional array with declared dimension (*), whose contents is set by RGMLTk (DGMLTk).

   Fk

   One-dimensional array with declared dimension (*), whose contents must be set by FSUBk as described below.

   X

   One-dimensional array which must be the same as the array X appearing as an actual argument in all calls to RGMLT1, RGMLT2, (DGMLT1, DGMLT2, ).

The subprogram RGMLTk (DGMLTk) which calls subroutine FSUBk sets the value of M and places in array Uk a set of M values of the variable . Then, if denotes the function which is to be integrated with respect to , it is the job of subroutine FSUBk to set X(k) to the appropriate value of , to compute for each of these values of (taking the remaining variables from X(k+1), ,X(n) respectively) and place the results in array Fk. (See Examples).

Method:

Integration with respect to each variable is performed by applying the 6- or 8-point Gaussian quadrature formula to each of the equal sub-intervals.

Notes:

1. The time needed to compute an n-dimensional integral by means of these subprograms is approximately

   

   This should be kept in mind when choosing the values of NIk.

2. The accuracy of a particular calculation can be estimated by repeating the integration with different values of NGk (e.g., 8 instead of 6) or by changing NIk, the latter usually being less economical.

Error handling:

Error D110.1: . A message is written on Unit 6, unless subroutine MTLSET (N002) has been called. Execution is halted on a STOP instruction.

Examples:

To calculate (in type REAL) the integral

using 6-point Gaussian quadrature over each of subdivisions of the corresponding interval. In the main program, write for instance

    ...
    EXTERNAL FSUB2
    DIMENSION X(2)
    Q2=RGMLT2(FSUB2,0.,1.,3,6,X)
    ...

For the SUBROUTINE subprograms FSUB2, FSUB1 write for instance

    SUBROUTINE FSUB2(M,U2,F2,X)
    EXTERNAL FSUB1
    DIMENSION U2(*),F2(*),X(2)
    DO 1 L = 1,M
    X(2)=U2(L)
    R=SQRT(X(2))
  1 F2(L)=R*EXP(X(2))*RGMLT1(FSUB1,0.,R,4,6,X)
    RETURN
    END
 
    SUBROUTINE FSUB1(M,U1,F1,X)
    DIMENSION U1(*),F1(*),X(2)
    DO 1 L = 1,M
    X(1)=U1(L)
  1 F1(L)=X(1)*SQRT(X(1)**2+X(2))
    RETURN
    END