intFcnSing¶

Integrates a function, which may have endpoint singularities, using a globally adaptive scheme based on Gauss-Kronrod rules.

Synopsis¶

intFcnSing (fcn, a, b)

Required Arguments¶

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

Return Value¶

An estimate of

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

If no value can be computed, NaN is returned.

Optional Arguments¶

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¶

This function is designed to handle functions with endpoint singularities. However, the performance on functions that are well-behaved at the endpoints is also quite good.

The function intFcnSing is a general-purpose integrator that uses a globally adaptive scheme in order 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, roundoff error is detected, the subintervals become too small, or the maximum number of subintervals allowed is reached. This function uses an extrapolation procedure known as the ɛ-algorithm.

The function intFcnSing is based on the subroutine QAGS by Piessens et al. (1983).

Examples¶

Example 1¶

The value of

\[\int_0^1 \ln (x) x^{-1/2} dx = -4\]

is estimated.

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


def fcn(x):
    ret = log(x) / sqrt(x)
    return ret


err_est = []

# Evaluate the integral
q = intFcnSing(fcn, 0.0, 1.0)

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

Output¶

integral  =     -4.000
exact     =     -4.000

Example 2¶

The value of

\[\int_0^1 \ln (x) x^{-1/2} dx = -4\]

is again estimated. The values of the actual and estimated errors are printed as well. Note that these numbers are machine dependent. Furthermore, usually the error estimate is pessimistic. That is, the actual error is usually smaller than the error estimate as is in this example.

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


def fcn(x):
    ret = log(x) / sqrt(x)
    return ret


err_est = []

# Evaluate the integral
q = intFcnSing(fcn, 0.0, 1.0, errEst=err_est)

# Print the result and the exact answer
exact = -4.0
exact_err = fabs(exact - q)
print("integral  = %10.3f\nexact     = %10.3f" % (q, exact))
print("error estimate = %e\nexact error    = %e"
      % (err_est[0], exact_err))

Output¶

integral  =     -4.000
exact     =     -4.000
error estimate = 2.273737e-13
exact error    = 2.620126e-14

Warning Errors¶

`IMSL_ROUNDOFF_CONTAMINATION`	Roundoff error, preventing the requested tolerance from being achieved, has been detected.
`IMSL_PRECISION_DEGRADATION`	A degradation in precision has been detected.
`IMSL_EXTRAPOLATION_ROUNDOFF`	Roundoff error in the extrapolation table, preventing the requested tolerance from being achieved, has been detected.

Fatal Errors¶

`IMSL_DIVERGENT`	Integral is probably divergent or slowly convergent.
`IMSL_PRECISION_DEGRADATION`	Integral is probably divergent or slowly convergent.
`IMSL_MAX_SUBINTERVALS`	The maximum number of subintervals allowed has been reached.
`IMSL_STOP_USER_FCN`	Request from user supplied function to stop algorithm. User flag = “#”.