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
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
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
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
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
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 [ainfinity]. 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
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
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
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