next up previous index
Next: D113 Adaptive Complex Up: CERNLIB Previous: D108 Trapezoidal Rule

D110 Gaussian Quadrature for Multiple Integrals

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

Function subprogram packages RGMLT and DGMLT compute an approximate value of an n-dimensional integral of the form
In = ∫anbndxnφn(xn)∫an-1(xn)bn-1(xn)dxn-1φn-1(xn-1,xn) ...
φ2(x2,...,xn)∫a1(x2,...,xn)b1(x2,...,xn)dx1φ1(x1,...,xn),
where 1 ≤n ≤6 .

Each subprogram integrates over only one variable. The integral In

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: RGMLT1, RGMLT2, RGMLT3, RGMLT4, RGMLT5, RGMLT6,
DGMLT1, DGMLT2, DGMLT3, DGMLT4, DGMLT5, DGMLT6
Files Referenced: Unit 6
External References: MTLMTR (N002), ABEND (Z035), User supplied SUBROUTINE subprograms

Usage:

  1. Let k be one of the integers 1,2,...,6 . 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 xk 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 xk .
    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 ≥n .

  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 (≤64) , 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 xk . Then, if fk(xk,...,xn) denotes the function which is to be integrated with respect to xk , it is the job of subroutine FSUBk to set X(k) to the appropriate value of xk , to compute fk for each of these values of xk

(taking the remaining variables xk+1,...,xn 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

    (NG1*NG2*NGn)*(NI1*NI2*NIn).

    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: NIk ≤0 . 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

Q2 = ∫01dx2x2ex20x2dx1x1x12+x2 = {13}(22-1)(e-2)

using 6-point Gaussian quadrature over each of n2= 3, n1= 4

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

D113



next up previous index
Next: D113 Adaptive Complex Up: CERNLIB Previous: D108 Trapezoidal Rule


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