Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

SGLGQ and DGLGQ--Numerical Quadrature Performed on a Function Using Gauss-Laguerre Quadrature

These functions approximate the integral of a real valued function over a semi-infinite interval, using the Gauss-Laguerre Quadrature method of specified order.

Table 164. Data Types

Syntax

On Entry

subf

is the user-supplied subroutine that evaluates the integrand function. The subroutine should be defined with three arguments: t, y, and n. For details, see Programming Considerations for the SUBF Subroutine.

Specified as: subf must be declared as an external subroutine in your application program. It can be whatever name you choose.

a

has the following meaning, where:

If b > 0, it is the lower limit of integration.

If b < 0, it is the upper limit of integration.

Specified as: a number of the data type indicated in Table 164.

b

is the scaling constant b for the exponential. Specified as: a number of the data type indicated in Table 164; b > 0 or b < 0.

n

is the order of the quadrature method to be used. Specified as: a fullword integer; n = 1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 20, 24, 32, 40, 48, or 64.

On Return

Function value

is the approximation of the integral. Returned as: a number of the data type indicated in Table 164.

Notes

1. Declare the DGLGQ function in your program as returning a long-precision real number. Declare the SGLGQ function, if necessary, as returning a short-precision real number.
2. The subroutine specified for subf must be declared as external in your program. Also, data types used by subf must agree with the data types specified by this ESSL subroutine. The variable x, described under Function, and the argument n correspond to the subf arguments t and n, respectively. For details on how to set up the subroutine, see Programming Considerations for the SUBF Subroutine.

Function

The integral is approximated for a real valued function over a semi-infinite interval, using the Gauss-Laguerre Quadrature method of specified order. The region of integration is:

(a, infinity)    if b > 0

(-infinity, a)    if b < 0

The method of order n is theoretically exact for integrals of the following form, where f is a polynomial of degree less than 2n:

The method of order n is a good approximation when your integrand is closely approximated by a function of the form f(x)e^-bx, where f is a polynomial of degree less than 2n. See references [26] and [92]. The result is returned as the function value.

Error Conditions

Computational Errors

None

Input-Argument Errors

1. b = 0
2. n is not an allowable value, as listed in the syntax for this argument.

Example 1

This example shows how to compute the integral of the function f given by:

f(x) = sin (3.0x)e^-1.5x

over the interval (-2.0, infinity), using the Gauss-Laguerre method with 20 points:

The user-supplied subroutine FUN1, which evaluates the integrand function, is coded in Fortran as follows:

   SUBROUTINE FUN1 (T,Y,N)
   INTEGER*4 N
   REAL*4 T(*),Y(*)
   DO 1 I=1,N
1     Y(I)=SIN(3.0*T(I))*EXP(-1.5*T(I))
   RETURN
   END

Program Statements and Input

EXTERNAL FUN1
        .
        .
        .
              SUBF     A     B    N
               |       |     |    |
XINT = SGLGQ( FUN1 , -2.0 , 1.5 , 20 )
        .
        .
        .

Output

XINT     =  5.891

Example 2

This example shows how to compute the integral of the function f given by:

f(x) = sin (3.0x)e^1.5x

over the interval (-infinity, -2.0), using the Gauss-Laguerre method with 20 points:

The user-supplied subroutine FUN2, which evaluates the integrand function, is coded in Fortran as follows:

   SUBROUTINE FUN2 (T,Y,N)
   INTEGER*4 N
   REAL*4 T(*),Y(*),TEMP
   DO 1 I=1,N
1     Y(I)=SIN(3.0*T(I))*EXP(1.5*T(I))
   RETURN
   END

Program Statements and Input

EXTERNAL FUN2
        .
        .
        .
              SUBF     A      B    N
               |       |      |    |
XINT = SGLGQ( FUN2 , -2.0 , -1.5 , 20 )
        .
        .
        .
FUN2     =  (see above)

Output

XINT     =  -0.011