IDL Analyst Reference Guide: Quadrature

IMSL_INTFCNHYPER

The IMSL_INTFCNHYPER function integrates a function on a hyper-rectangle as follows:

Syntax

Result = IMSL_INTFCNHYPER(f, a, b [, ERR_ABS=value] [, ERR_EST=variable] [, ERR_REL=value] [, MAX_EVALS=value])

Return Value

The value of the hyper-rectangle function is returned. If no value can be computed, the floating-point value NaN (Not a Number) is returned.

Arguments

f

A scalar string specifying the name of a user-supplied function to be integrated. The function f accepts an array of data points at which the function is to be evaluated and returns the scalar value of the function.

a

A vector specifying the lower limit of integration.

b

A vector specifying the upper limit of integration.

Keywords

ERR_ABS

Set this keyword to a value specifying the accuracy desired. Default: ERR_ABS=SQRT(e), where e is the machine precision

ERR_EST

Set this keyword equal to a named variable that will contain an estimate of the absolute value of the error.

ERR_REL

Set this keyword to a value specifying the relative accuracy desired. Default: ERR_REL=SQRT(e), where e is the machine precision

MAX_EVALS

Set this keyword to a scalar value specifying the number of evaluations allowed. Default: MAX_EVALS = 1,000,000 for n £ 2 and MAX_EVALS = 256ⁿ for n > 2, where n is the number of independent variables of f.

Discussion

The IMSL_INTFCNHYPER function approximates the following n-dimensional iterated integral:

An estimate of the error is returned in the optional keyword ERR_EST. The approximation is achieved by iterated applications of product Gauss formulas. The integral is first estimated by a two-point, tensor-product formula in each direction. Then, for ( i = 0, ..., n – 1 ), the function calculates a new estimate by doubling the number of points in the i-th direction, but halving the number immediately afterwards if the new estimate does not change appreciably. This process is repeated until either one complete sweep results in no increase in the number of sample points in any dimension, the number of Gauss points in one direction exceeds 256, or the number of function evaluations needed to complete a sweep exceeds MAX_EVALS.

Example

This example computes the integral of:

on an expanding cube. The values of the error estimates are machine dependent. The exact integral over R is p^3/2.

.RUN  
; Define the function to be integrated.  
FUNCTION f, x  
   RETURN, EXP(-TOTAL(x^2))  
END  
  
limit = !Pi^1.5  
; Compute the exact value of the integral.  
PM, '    Limit:', limit  
       Limit:      5.56833  
FOR i = 1, 6 DO  BEGIN $  
   a = [-i/2., -i/2., -i/2.] &$  
   b = [i/2., i/2., i/2.] &$  
   ans = IMSL_INTFCNHYPER('f', a, b) &$  
   PRINT, 'integral = ', ans, ' limit = ', limit  
   ; Compute values of the integral over expanding cubes and  
   ; output the results after each call to IMSL_INTFCNHYPER.  
    integral = 0.785213 limit = 5.56833  
    integral =  3.33231 limit = 5.56833  
    integral =  5.02107 limit = 5.56833  
    integral =  5.49055 limit = 5.56833  
    integral =  5.56135 limit = 5.56833  
    integral =  5.56771 limit = 5.56833