INTFCN Function

PV-WAVE Advantage > IMSL Mathematics Reference Guide > Quadrature > INTFCN Function

INTFCN Function

Integrates a user-supplied function. Using different combinations of keywords and parameters, one of several types of integration can be performed including the following:

Default method

Integration of functions based on Gauss-Kronrod rules

Integration of functions with singular points given

Integration of functions with algebraic-logarithmic singularities

Integration of functions over an infinite or semi-infinite interval

Integration of functions containing a sine or cosine factor

Computation of Fourier sine and cosine transforms

Integrals in the Cauchy principle value sense

Integration of smooth functions using a nonadaptive rule

Computation of two-dimensional iterated integrals

Default Method

Usage

result = INTFCN(f, a, b)

Input Parameters

In the default case, the following three parameters are required. If another method of integration is desired, a combination of these parameters, along with additional parameters and keywords must be specified. For a description of the additional parameters and keywords to use for specific methods of integration, see "Optional Methods of Integration".

f—Scalar string specifying the name of a user-supplied function to be integrated. The function f accepts one scalar parameter and returns a single scalar of the same type.

a—Scalar expression specifying the lower limit of integration.

b—Scalar expression specifying the upper limit of integration.

Returned Value

result—An estimate of the desired integral. If no value can be computed, then NaN (Not a Number) is returned.

Global Keywords

The following seven keywords can be used in any combination with each method of integration except the nonadaptive method, which is triggered by the keyword Smooth. Because of this, these global keywords are documented here only and referred to within the Method Keywords subsections of the Optional Method of Integrations section of this routine.

Input Keywords

Double—If present and nonzero, double precision is used.

Err_Abs—Absolute accuracy desired. Default: Err_Abs = SQRT(ε), where ε is the machine precision

Err_Rel—Relative accuracy desired. Default: Err_Rel = SQRT(ε), where ε is the machine precision

Max_Subinter—Number of subintervals allowed. Default: Max_Subinter = 500

Output Keywords

Err_Est—Named variable in which an estimate of the absolute value of the error is stored.

N_Subinter—Named variable into which the number of subintervals generated is stored.

N_Evals—Named variable into which the number of evaluations of f is stored.

Discussion of Default Method

The default method used by INTFCN is a general-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It subdivides the interval [a, b] and uses a 21-point Gauss-Kronrod rule to estimate the integral over each subinterval. The error for each subinterval is estimated by comparison with the 10-point Gauss quadrature rule. The subinterval with the largest estimated error is then bisected, and the same procedure is applied to both halves. The bisection process is continued until either the error criterion is satisfied, the roundoff error is detected, the subintervals become too small, or the maximum number of subintervals allowed is reached. This method uses an extrapolation procedure known as the ε-algorithm. This method is based on the subroutine QAGS by Piessens et al. (1983).

Should the default method fail to produce acceptable results, consider one of the more specialized methods available by using method-specific keywords for this function.

Example

An estimate of:

is computed, then compared to the actual value.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, x^2

- END

; Call INTFCN to compute the integral.

ans = INTFCN('f', 0, 3)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints the following:

; Computed Answer:

;     9.00000

PM, 'Exact - Computed:', 3^2 - ans

; PV-WAVE prints the following:

; Exact - Computed:

;      0.00000

Warning Errors

MATH_ROUNDOFF_CONTAMINATION—Roundoff error, preventing the requested tolerance from being achieved, has been detected.

MATH_PRECISION_DEGRADATION—Degradation in precision has been detected.

MATH_EXTRAPOLATION_ROUNDOFF—Roundoff error in the extrapolation table, preventing requested tolerance from being achieved, has been detected.

MATH_EXTRAPOLATION_PROBLEMS—Extrapolation table, constructed for convergence acceleration of the series formed by the integral contributions of the cycles, does not converge to the requested accuracy.

MATH_BAD_INTEGRAND_BEHAVIOR—Bad integrand behavior occurred in one or more cycles.

Fatal Errors

MATH_DIVERGENT—Integral is probably divergent or slowly convergent.

MATH_MAX_SUBINTERVALS—Maximum number of subintervals allowed has been reached.

MATH_MAX_CYCLES—Maximum number of cycles allowed has been reached.

MATH_MAX_STEPS—Maximum number of steps allowed have been taken. The integrand is too difficult for this routine.

Optional Methods of Integration

By specifying different sets of parameters and/or keywords, a number of different types of integration can be performed. Internally, the method to be used is determined by examining the combination of parameters and/or keywords used in the call to INTFCN. To specify a specific method of integration, refer to the appropriate discussion.

Integration of Functions Based on Gauss-Kronrod Rules

This method integrates functions using a globally adaptive scheme based on Gauss-Kronrod rules.

Usage

Triggered by the use of keyword Rule.

result = INTFCN(f, a, b, Rule = rule)

Returned Value

result—The value of:

is returned. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in the "Global Keywords", the following keyword is available:

Rule—If specified, the integral is computed using a globally adaptive scheme based on Gauss-Kronrod rules. Refer to Table 5-2: Corresponding Gauss-Kronrod Rules for a list of the Gauss-Kronrod Rules that correspond to the Rule keyword.

Corresponding Gauss-Kronrod Rules
Rule	Gauss-Kronrod Rule
1	7-15 points
2	10-21 points
3	15-31 points
4	20-41 points
5	25-51 points
6	30-61 points

Discussion

This method is a general-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It subdivides the interval [a, b] and uses a (2k+1)-point Gauss-Kronrod rule to estimate the integral over each subinterval. The error for each subinterval is estimated by comparison with the k-point Gauss quadrature rule. The subinterval with the largest estimated error is then bisected, and the same procedure is applied to both halves. The bisection process is continued until either the error criterion is satisfied, roundoff error is detected, the subintervals become too small, or the maximum number of subintervals allowed is reached. This method is based on the subroutine QAG by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed. Since the integrand is oscillatory, Rule = 6 is used. The exact value is 0.50406706. The values of the actual and estimated error are machine dependent.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, SIN(1/x)

- END

; Call INTFCN, with Rule = 6, to compute the integral based on the

; specified Gauss-Kronrod rule.

ans = INTFCN('f', 0, 1, Rule = 6)

; Output the results.

PM, 'Computed Answer:',ans

; PV-WAVE prints the following:

; Computed Answer:

;     0.504051

exact = .50406706

PM, 'EXACT - COMPUTED:', exact - ans

; PV-WAVE prints the following:

; Exact - Computed:

;     1.62125e-05

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Functions with Singular Points Given

This method integrates functions with singularity points given.

Usage

Requires the use of keyword Sing_Pts.

result = INTFCN(f, a, b, Sing_Pts = points)

Returned Value

result—The value of:

is returned. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in the "Global Keywords", the following keyword is available:

Sing_Pts—If present, specifies the abscissas of the singularities. These values should be interior to the interval [a, b].

Discussion

This method is a special-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It subdivides the interval [a, b] into N+1 user-supplied subintervals, where N is the number of singular points, and uses a 21-point Gauss-Kronrod rule to estimate the integral over each subinterval. The error for each subinterval is estimated by comparison with the 10-point Gauss quadrature rule. The subinterval with the largest estimated error is then bisected, and the same procedure is applied to both halves. The bisection process is continued until either the error criterion is satisfied, the roundoff error is detected, the subintervals become too small, or the maximum number of subintervals allowed is reached. This method uses an extrapolation procedure known as the ε-algorithm. This method is based on the subroutine QAGP by Piessens et al. (1983).

Example

The value of:

is computed. The values of the actual and estimated error are machine dependent. Note that this subfunction never evaluates the user-supplied function at the user-supplied breakpoints.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, x^3 * ALOG(ABS((x^2 - 1) * $

   - (x^ 2 - 2)))

- END

; Call INTFCN using keyword Sing_Pts to specify the singular

; points.

ans = INTFCN('f', 0, 3, $

Sing_Pts = [1, SQRT(2)], N_Evals = nevals)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints the following:

; Computed Answer:

;      52.7408

exact = 61 * ALOG(2) + (77/4.) * ALOG(7) - 27

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints the following:

; Exact - Computed:

; -2.67029e-05

PM, 'Number of Function Evaluations:', nevals

; PV-WAVE prints the following:

; Number of Function Evaluations:

;         819

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Functions with Algebraic-logarithmic Singularities

This method integrates functions with algebraic-logarithmic singularities.

Method Input Parameters

alpha—The strength of the singularity at a. Must be greater than –1.

beta—Strength of the singularity at b. Must be greater than –1.

Usage

Triggered by the use of the parameters alpha and beta and one of the keywords below in addition to f, a, and b.

result = INTFCN(f, a, b, alpha, beta, /Algebraic)

result = INTFCN(f, a, b, alpha, beta, /Alg_Left_Log)

result = INTFCN(f, a, b, alpha, beta, /Alg_Log)

result = INTFCN(f, a, b, alpha, beta, /Alg_Right_Log)

Returned Value

result—The value of:

is returned, where w (x) is defined by one of the keywords below. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in the "Global Keywords", the following keywords are available. Exactly one of the following keywords must be specified:

Algebraic—If present and nonzero, uses the weight function

. This is the default weight function for this method.

Alg_Left_Log—If present and nonzero, uses the weight function

Alg_Log—If present and nonzero, uses the weight function

Alg_Right_Log—If present and nonzero, uses the weight function

Discussion

This method is a special-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It computes integrals whose integrands have the special form w (x) f (x), where w (x) is a weight function. A combination of modified Clenshaw-Curtis and Gauss-Kronrod formulas is employed. This method is based on the subroutine QAWS, which is fully documented by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, SQRT((1 + x))

- END

ans = $

; Call INTFCN with keyword Alg_Left_Log set and values for the

; method parameters alpha and beta.

INTFCN('f', 0, 1, /Alg_Left_Log, 1.0, .5 )

; Output the results.

PM, 'Computed Answer:', ans

; Computed Answer: -0.213395

exact = (3 * ALOG(2) - 4)/9

PM, 'Exact - Computed:', exact - ans

; Exact - Computed: 1.49012e-08

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Functions Over an Infinite or Semi-infinite Interval

This method integrates functions over an infinite or semi-infinite interval.

Method Input Parameters

bound—The finite limit of integration. If either of the keywords Inf_Bound or Bound_Inf are specified, this parameter is required.

Usage

Triggered by the presence of the function f, a bound bound, and one of the keywords Inf_Inf, Inf_Bound, or Bound_Inf.

result = INTFCN(f, /Inf_Inf)

result = INTFCN(f, bound, /Inf_Bound)

result = INTFCN(f, bound, /Bound_Inf)

Returned Value

result—The value of:

is returned, where a and b are appropriate integration limits. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in the "Global Keywords", the following keywords are available (exactly one of the following keywords must be specified):

Inf_Inf—If present and nonzero, integrates a function over the range ( –infinity, infinity).

Inf_Bound—If present and nonzero, integrates a function over the range ( –infinity,bound ).

Bound_Inf—If present and nonzero, integrates a function over the range ( bound, infinity).

Discussion

This method is a special-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It initially transforms an infinite or semi-infinite interval into the finite interval [0, 1]. It then uses the same strategy that is used when specifying the keyword Sing_Pts. This method is based on the subroutine QAGI by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, ALOG(x)/(1 + (10 * x)^2)

- END

; Call INTFCN with keyword Bound_Inf set. Notice that only lower

; limit of integration is given.

ans = INTFCN('f', 0, /Bound_Inf)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: -0.361689

exact = -!Pi * ALOG(10)/20

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: 5.96046e-08

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Functions Containing a Sine or Cosine Factor

This method integrates functions containing a sine or a cosine factor.

Method Input Parameters

omega—The frequency of the trigonometric weighting function.

Usage

Triggered by the use of parameter omega and one of the keywords Sine or Cosine in addition to f, a, and b.

result = INTFCN(f, a, b, omega, /Sine)

result = INTFCN(f, a, b, omega, /Cosine)

Returned Value

result—The value of:

where the weight function w (ωx) is defined by the keywords below, is returned. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in "Global Keywords", the following keywords are available (exactly one of the following keywords must be specified):

Sine—If present and nonzero, sin (ωx) is used for the integration weight function.

Cosine—If present and nonzero, cos (ωx) is used for the integration weight function.

Max_Moments—A scalar expression specifying an upper bound on the number of Chebyshev moments that can be stored. Increasing (decreasing) this number may increase (decrease) execution speed and space used. Default: Max_Moments = 21

Discussion

This method is a special-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It computes integrals whose integrands have the special form w(x) f(x), where w(x) is either cos(ωx) or sin(ωx). Depending on the subinterval’s length in relation to the size of ω, either a modified Clenshaw-Curtis procedure or a Gauss-Kronrod 7/15 rule is employed to approximate the integral on a subinterval. This method is based on the subroutine QAWO by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed. The following is the exact answer:

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, x^2

- END

; Call INTFCN with Sine set and value for method parameter omega.

ans = INTFCN('f', 0, 1, 3 * !Pi, /Sine)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: 0.101325

exact = ((3 * !Pi)^2 - 2)/((3 * !pi)^3) - 2/(3 * !Pi)^3

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: 0.00000

Error Handling

See "Warning Errors" and "Fatal Errors".

Computation of Fourier Sine or Cosine Transforms

This method computes Fourier sine or cosine transforms.

Method Input Parameters

omega—The frequency of the trigonometric weighting function.

Usage

Triggered by the use of parameter omega and one of the keywords Sine or Cosine in addition to f and a.

result = INTFCN(f, a, omega, /Sine)

result = INTFCN(f, a, omega, /Cosine)

Returned Value

result—The value of:

where the weight function w (ωx) is defined by the keywords below, is returned. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in "Global Keywords", the following keywords are available (exactly one of the keywords Sine or Cosine must be specified):

Sine—If present and nonzero, sin (ωx) is used for the integration weight function.

Cosine—If present and nonzero, cos (ωx) is used for the integration weight function.

Max_Moments—Number of subintervals allowed in the partition of each cycle. Default: Max_Moments = 21

Method Output Keywords

N_Cycles—Named variable into which the number of cycles is stored.

Discussion

This method is a special-purpose integrator that uses a globally adaptive scheme to reduce the absolute error. It computes integrals whose integrands have the special form w (x) f (x), where w (x) is either cos (ωx) or sin (ωx). The integration interval is always semi-infinite of the form [a, infinity]. This method is based on the subroutine QAWF by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed. Notice that the function is coded to protect for the singularity at zero.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - IF (x EQ 0) THEN RETURN, x $

      - ELSE RETURN, 1/SQRT(x)

- END

; Call INTFCN with keyword Cosine set and a value for the method

; specific parameter omega.

ans = INTFCN('f', 0, !Pi/2, /Cosine)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: 0.101325

exact = 1.0

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: 0.00000

Error Handling

See "Warning Errors" and "Fatal Errors".

Integrals in the Cauchy Principle Value Sense

This method computes integrals of the form:

in the Cauchy principal value sense.

Method Input Parameters

c—The singular point must not equal a or b.

Usage

Triggered by the use of parameter c and keyword Cauchy in addition to f, a, and b.

result = INTFCN(f, a, b, c, /Cauchy)

Returned Value

result—The value of:

is returned. If no value can be computed, NaN is returned.

Method Input Keywords

In addition to the keywords listed in "Global Keywords", the following keyword is available.

Cauchy—If present and nonzero, computes integrals of the form:

in the Cauchy principal value sense.

Discussion

This method uses a globally adaptive scheme in an attempt to reduce the absolute error. It computes integrals whose integrands have the special form w (x) f (x), where w (x) = 1/(x – c). If c lies in the interval of integration, then the integral is interpreted as a Cauchy principal value. A combination of modified Clenshaw-Curtis and Gauss-Kronrod formulas is employed. The method is an implementation of the subroutine QAWC by Piessens et al. (1983).

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The Cauchy principal value of:

is computed.

.RUN

; Define the function to be integrated.

- FUNCTION f, x

   - RETURN, 1/(5 * x^3 + 6)

- END

; Call INTFCN with keyword Cauchy set.

ans = INTFCN('f', -1, 5, 0, /Cauchy)

; Output the results.

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: -0.0899440

exact = ALOG(125/631.)/18

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: 1.49012e-08

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Smooth Functions Using Nonadaptive Rule

This method integrates smooth functions using a nonadaptive rule.

Usage

Triggered by the use of keyword Smooth in addition to f, a, and b.

result = INTFCN(f, a, b, /Smooth)

Returned Value

result—The value of:

is returned. If no value can be computed, NaN is returned.

Method Input Keywords

Because this method is nonadaptive, there are fewer options with the algorithm. For this method, all keywords described in "Global Keywords" do not apply. A complete list of the available keywords is given below. This method requires the use of keyword Smooth.

Smooth—If present and nonzero, uses a nonadaptive rule to compute the integral.

Double—If present and nonzero, uses double precision.

Err_Abs—Absolute accuracy desired. Default: Err_Abs = SQRT(ε), where ε is the machine precision

Err_Rel—Relative accuracy desired. Default: Err_Rel = SQRT(ε), where ε is the machine precision

Method Output Keywords

Err_Est—Named variable into which an estimate of the absolute value of the error is stored.

Discussion

This method is designed to integrate smooth functions. It implements a nonadaptive quadrature procedure based on nested Paterson rules of order 10, 21, 43, and 87. These rules are positive quadrature rules with degree of accuracy 19, 31, 64, and 130, respectively. This method applies these rules successively, estimating the error until either the error estimate satisfies the user-supplied constraints or the last rule is applied.

This method is not very robust, but for certain smooth functions, it can be efficient. This method is based on the subroutine QNG by Piessens et al. (1983). If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

The value of:

is computed.

.RUN

; Define the function to integrate.

- FUNCTION f, x

   - RETURN, x * EXP(x)

- END

; Call INTFCN with keyword Smooth set.

ans = INTFCN('f', 0, 2, /Smooth)

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: 8.38906

exact = EXP(2) + 1

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: 9.53674e-07

Error Handling

See "Warning Errors" and "Fatal Errors".

Integration of Two-dimensional Iterated Integrals

This method integrates two-dimensional iterated integrals.

Method Input Parameters

f—Scalar string specifying the name of a user-supplied PV‑WAVE function to be integrated. Function f accepts two scalar parameters and returns a single scalar of the same type.

a—Scalar expression specifying the lower limit of the outer integral.

b—Scalar expression specifying the upper limit of the outer integral.

h—Scalar string specifying the name of a user-supplied PV‑WAVE function used to evaluate the lower limit of the inner integral. Function h accepts one scalar parameter and returns a single scalar of the same type.

g—Scalar string specifying the name of a user-supplied PV‑WAVE function used to evaluate the upper limit of the inner integral. Function g accepts one scalar parameter and returns a single scalar of the same type.

Usage

Triggered by the use of the parameters g and h and keyword Two_Dimensional in addition to f, a, and b.

result = INTFCN(f, a, b, g, h, /Two_Dimensional)

Returned Value

result—The value of:

is returned. If no value can be computed, NaN is returned.

Method Keywords

In addition to the keywords listed in "Global Keywords", the following keyword is available and must be specified for this method:

Two_Dimensional—If present and nonzero, integrates a two-dimensional iterated integral.

Discussion

This method approximates the following two-dimensional iterated integral:

The lower-numbered rules are used for less smooth integrands, while the higher-order rules are more efficient for smooth (oscillatory) integrands.

If this method is used, the function should be coded to protect endpoint singularities if they exist.

Example

In this example, the value of the integral:

is computed.

.RUN

; Define the function to be integrated.

- FUNCTION f, x, y

   - RETURN, SIN(x + y)

- END

.RUN

; Define the function for the lower limit of the inner integral.

- FUNCTION g, x

   - RETURN, x

- END

.RUN

; Define the function for the upper limit of the inner integral.

- FUNCTION h, x

   - RETURN, 2 * x

- END

; Call INTFCN with keyword Two_Dimensional set and the names

; of the functions defining the limits of the inner integral.

ans = INTFCN('f',0,1,'g','h',/Two_Dimensional)

PM, 'Computed Answer:', ans

; PV-WAVE prints: Computed Answer: 0.407609

exact = -SIN(3)/3 + SIN(2)/2

PM, 'Exact - Computed:', exact - ans

; PV-WAVE prints: Exact - Computed: -5.96046e-08

Error Handling

See "Warning Errors" and "Fatal Errors".