intFcn¶

Integrates a function using a globally adaptive scheme based on Gauss-Kronrod rules.

Synopsis¶

intFcn (fcn, a, b)

Required Arguments¶

float fcn (float x) (Input): User-supplied function to be integrated.
float a (Input): Lower limit of integration.
float b (Input): Upper limit of integration.

Return Value¶

The value of

\[\int_a^b\textit{fcn}(x)dx\]

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

Optional Arguments¶

rule, int (Input): Choice of quadrature rule.

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

Default: rule = 1

errAbs, float (Input)

Absolute accuracy desired.

Default: \(\mathit{errAbs} = \sqrt{\varepsilon}\) where ɛ is the machine precision

errRel, float (Input)

Relative accuracy desired.

Default: \(\mathit{errRel} = \sqrt{\varepsilon}\) where ɛ is the machine precision

errEst (Output)

An estimate of the absolute value of the error.

maxSubinter, int (Input)

Number of subintervals allowed.

Default: maxSubinter = 500

nSubinter (Output)

The number of subintervals generated.

nEvals, (Output)

The number of evaluations of fcn.

Description¶

The function 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 (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. The function intFcn is based on the subroutine QAG by Piessens et al. (1983).

Should intFcn fail to produce acceptable results, consider one of the more specialized functions documented in this chapter.

Examples¶

Example 1¶

The value of

\[\int_0^2 xe^xdx = e^2 + 1\]

is computed. Since the integrand is not oscillatory, all of the default values are used. The values of the actual and estimated error are machine dependent.

from __future__ import print_function
from numpy import *
from pyimsl.math.intFcn import intFcn

exact = 8.389


def fcn(x):
    res = x * (exp(x))
    return res


# Evaluate the integral
q = intFcn(fcn, 0.0, 2.0)

# Print the result and the exact answer
print("integral = %10.3f\nexact = %10.3f" % (q, exact))

Output¶

integral =      8.389
exact =      8.389

Example 2¶

The value of

\[\int_0^1 \sin(1/x)dx\]

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.

from __future__ import print_function
from numpy import *
from pyimsl.math.intFcn import intFcn

err_est = []
err_abs = 0.0001
exact = 0.50406706


def fcn(x):
    # Compute sin(1/x), avoiding division by zero
    if ((x) > 1.0e-5):
        res = sin(1.0 / (x))
    else:
        res = 0.0
    return res


# Integrate fcn(x) from 0 to 1
q = intFcn(fcn, 0.0, 1.0, errAbs=err_abs,
           rule=6, errEst=err_est)
error = q - exact

# Print the result and the exact answer
print(" integral = %10.3f\n    exact = %10.3f\n    error = %10.3f" %
      (q, exact, error))
print("   err_est = %g" % err_est[0])

Output¶

 integral =      0.504
    exact =      0.504
    error =      0.000
   err_est = 9.86967e-05

Warning Errors¶

`IMSL_ROUNDOFF_CONTAMINATION`	Roundoff error has been detected. The requested tolerances, “`errAbs`” = # and “`errRel`” cannot be reached.
`IMSL_PRECISION_DEGRADATION`	Precision is degraded due to too fine a subdivision relative to the requested tolerance. This may be due to bad integrand behavior in the interval (#,#). Higher precision may alleviate this problem.

Fatal Errors¶

IMSL_MAX_SUBINTERVALS The maximum number of subintervals allowed “maxSubinter” has been reached.

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.