intFcnSingPts¶

Integrates a function with singularity points given.

Synopsis¶

intFcnSingPts (fcn, a, b, points)

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.
float points[] (Input): The abscissas of the singularities. These values should be interior to the interval [a, b].

Return Value¶

The value of

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

is returned. 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¶

The function intFcnSingPts is a special-purpose integrator that uses a globally adaptive scheme in order to reduce the absolute error. It subdivides the interval [a, b] into npoints + 1 user-supplied subintervals 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 intFcnSingPts is based on the subroutine QAGP by Piessens et al. (1983).

Examples¶

Example 1¶

The value of

\[\int_0^3 x^3 \ln \left| \left(x^2-1\right)\left(x^2-2\right)\right|dx = 61 \ln 2 + \frac{77}{4} \ln 7 - 27\]

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

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


def fcn(x):
    return x * x * x * (log(fabs((x * x - 1.) * (x * x - 2.))))


# Set singular points
points = array([1.0, sqrt(2.0)])

# Evaluate the integral
err_est = []
n_evals = []
q = intFcnSingPts(fcn, 0.0, 3.0, points)

# Print the result and the exact answer
exact = 61. * log(2.) + (77. / 4) * log(7.) - 27.
print("integral  = %10.3f\nexact     = %10.3f" % (q, exact))

Output¶

integral  =     52.741
exact     =     52.741

Example 2¶

The value of

\[\int_0^3 x^3 \ln \left| \left(x^2-1\right)\left(x^2-2\right)\right|dx = 61 \ln 2 + \frac{77}{4} \ln 7 - 27\]

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

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


def fcn(x):
    return x * x * x * (log(fabs((x * x - 1.) * (x * x - 2.))))


# Evaluate the integral and get the error estimate and the
# number of evaluations
err_est = []
points = array([1.0, sqrt(2.0)])
n_evals = []
q = intFcnSingPts(fcn, 0.0, 3.0, points,
                  errEst=err_est,
                  nEvals=n_evals)

# Print the result and the exact answer
exact = 61. * log(2.) + (77. / 4) * log(7.) - 27.
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))
print("The number of function evaluations  =  %d" % n_evals[0])

Output¶

integral  =     52.741
exact     =     52.741
error estimate   = 2.624763e-07
exact error      = 3.552714e-14
The number of function evaluations  =  1029

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 has been detected in the extrapolation table. The tolerances, “`errAbs`” = # and “`errRel`” = # cannot be reached.

Fatal Errors¶

IMSL_DIVERGENT Integral is probably divergent or slowly convergent.

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.