API Reference¶

Public API¶

integrationdifferentiation — SciPy-backed numerical integration and differentiation.

All public symbols are re-exported from integrationdifferentiation.core so callers can simply write:

from integrationdifferentiation import integrate_adaptive, gauss_quadrature_rule

This package is a self-contained set of integration and differentiation utilities inspired by the IMSL Fortran Numerical Library Integration and Differentiation chapter, backed by scipy.integrate and numpy.

class integrationdifferentiation.IntegrationResult(value, error, success, message, n_fev=0)¶

Bases: object

Result of a numerical integration procedure.

Parameters:

value (float)
error (float)
success (bool)
message (str)
n_fev (int)

value¶

Integral estimate.

Type:: float

error¶

Estimated absolute error bound.

Type:: float

success¶

Whether the integration succeeded.

Type:: bool

message¶

Description of the termination condition.

Type:: str

n_fev¶

Number of function evaluations performed.

Type:: int

class integrationdifferentiation.GaussQuadratureRule(points, weights, n)¶

Bases: object

A Gauss quadrature rule (abscissas and weights).

Parameters:

points (ndarray)
weights (ndarray)
n (int)

points¶

Quadrature abscissas sorted in ascending order.

Type:: np.ndarray

weights¶

Corresponding quadrature weights.

Type:: np.ndarray

n¶

Number of quadrature points.

Type:: int

integrationdifferentiation.integrate_adaptive(f, a, b, abs_tol=1.49e-08, rel_tol=1.49e-08, max_fev=500)¶

General-purpose adaptive quadrature (IMSL QDAG).

Uses global adaptive integration (QUADPACK DQAGSE) via scipy.integrate.quad(). Suitable for smooth integrands with finite or infinite limits.

Parameters:

f (Callable) – Integrand: f(x) -> float.
a (float) – Lower integration limit (may be -np.inf).
b (float) – Upper integration limit (may be np.inf).
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
rel_tol (float) – Relative error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol or rel_tol is not positive.

integrationdifferentiation.integrate_adaptive_endpoint_singular(f, a, b, abs_tol=1.49e-08, max_fev=500)¶

Adaptive quadrature with endpoint singularities (IMSL QDAGS).

Uses scipy.integrate.quad() which automatically handles integrable algebraic and logarithmic endpoint singularities (QUADPACK DQAGSE).

Parameters:

f (Callable) – Integrand: f(x) -> float. May diverge at a or b provided the singularity is integrable.
a (float) – Lower integration limit.
b (float) – Upper integration limit.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_adaptive_singular_points(f, a, b, points, abs_tol=1.49e-08, max_fev=500)¶

Adaptive quadrature with specified interior singularity points (IMSL QDAGP).

Splits the integration interval at each specified singularity point and applies adaptive integration to each sub-interval using QUADPACK DQAGPE.

Parameters:

f (Callable) – Integrand: f(x) -> float.
a (float) – Lower integration limit.
b (float) – Upper integration limit.
points (list[float]) – Interior points where f may have a singularity or discontinuity. Must satisfy a < points[i] < b.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive or any point is outside (a, b).

integrationdifferentiation.integrate_adaptive_1d(f, a, b, abs_tol=1.49e-08, max_fev=500)¶

Adaptive 1D quadrature with possible internal singularity (IMSL QDAG1D).

Uses scipy.integrate.quad() with a finer local control on singularity handling. Equivalent to QDAG but exposes a separate entry point consistent with the IMSL function family.

Parameters:

f (Callable) – Integrand: f(x) -> float.
a (float) – Lower integration limit.
b (float) – Upper integration limit.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_adaptive_infinite(f, a=-inf, b=inf, abs_tol=1.49e-08, max_fev=500)¶

Adaptive quadrature on an infinite or semi-infinite interval (IMSL QDAGI).

Uses scipy.integrate.quad() with ±np.inf limits. The integration is transformed to a finite interval internally by QUADPACK DQAGIE.

Parameters:

f (Callable) – Integrand: f(x) -> float. Must decay sufficiently fast for the integral to converge.
a (float) – Lower limit. Use -np.inf for a fully infinite lower bound. Defaults to -np.inf.
b (float) – Upper limit. Use np.inf for a fully infinite upper bound. Defaults to np.inf.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_weighted_oscillatory(f, a, b, omega, weight='sin', abs_tol=1.49e-08, max_fev=500)¶

Adaptive weighted oscillatory quadrature (IMSL QDAWO).

Computes ∫_a^b f(x) * w(x) dx where w(x) is either sin(omega*x) or cos(omega*x). Uses QUADPACK DQAWO via scipy.integrate.quad() with the appropriate weight and wvar arguments.

Parameters:

f (Callable) – Non-oscillatory part of the integrand: f(x) -> float.
a (float) – Lower integration limit (must be finite).
b (float) – Upper integration limit (must be finite).
omega (float) – Oscillation frequency omega in the weight function.
weight (str) – Weight type: 'sin' for sin(omega*x) or 'cos' for cos(omega*x). Defaults to 'sin'.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If weight is not 'sin' or 'cos', or if abs_tol is not positive, or if limits are infinite.

integrationdifferentiation.integrate_weighted_fourier(f, a, omega, weight='sin', abs_tol=1.49e-08, max_fev=500)¶

Adaptive weighted Fourier quadrature on [a, ∞) (IMSL QDAWF).

Computes ∫_a^∞ f(x) * w(x) dx using QUADPACK DQAWFE via scipy.integrate.quad() with semi-infinite limits and trigonometric weight.

Parameters:

f (Callable) – Non-oscillatory part of the integrand: f(x) -> float. Must decay to zero as x → ∞.
a (float) – Lower integration limit (must be finite).
omega (float) – Oscillation frequency in the weight function.
weight (str) – Weight type: 'sin' for sin(omega*x) or 'cos' for cos(omega*x). Defaults to 'sin'.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If weight is not 'sin' or 'cos', abs_tol is not positive, or a is infinite.

integrationdifferentiation.integrate_weighted_algebraic(f, a, b, alpha, beta, weight_type=1, abs_tol=1.49e-08, max_fev=500)¶

Adaptive quadrature with algebraic-logarithmic endpoint weight (IMSL QDAWS).

Computes ∫_a^b f(x) * w(x) dx where w(x) has algebraic and/or logarithmic endpoint singularities. Uses QUADPACK DQAWS via scipy.integrate.quad().

Weight types (weight_type):

1: (x-a)^alpha * (b-x)^beta
2: (x-a)^alpha * (b-x)^beta * log(x-a)
3: (x-a)^alpha * (b-x)^beta * log(b-x)
4: (x-a)^alpha * (b-x)^beta * log(x-a) * log(b-x)

Parameters:

f (Callable) – Integrand function: f(x) -> float.
a (float) – Lower integration limit (finite).
b (float) – Upper integration limit (finite), must satisfy b > a.
alpha (float) – Exponent for the (x-a) factor. Must be > -1.
beta (float) – Exponent for the (b-x) factor. Must be > -1.
weight_type (int) – One of 1, 2, 3, 4. Defaults to 1.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Integral estimate with error bound and diagnostics.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If weight_type is not in {1, 2, 3, 4}, alpha <= -1, beta <= -1, abs_tol not positive, or limits are non-finite.

integrationdifferentiation.integrate_cauchy_principal_value(f, a, b, c, abs_tol=1.49e-08, max_fev=500)¶

Adaptive Cauchy principal value integration (IMSL QDAWC).

Computes the Cauchy principal value of ∫_a^b f(x)/(x-c) dx via scipy.integrate.quad() with weight='cauchy'.

Parameters:

f (Callable) – Numerator function: f(x) -> float. The integrand is f(x) / (x - c).
a (float) – Lower integration limit (finite).
b (float) – Upper integration limit (finite), must satisfy b > a.
c (float) – Cauchy singularity point; must satisfy a < c < b.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.
max_fev (int) – Maximum number of function evaluations. Defaults to 500.

Returns:

Principal-value integral estimate with error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If c is not strictly inside (a, b), abs_tol is not positive, or limits are non-finite.

integrationdifferentiation.integrate_nonadaptive(f, a, b, n=10)¶

Non-adaptive quadrature using a Gaussian rule (IMSL QDNG).

Uses fixed-order Gauss-Legendre quadrature via scipy.integrate.fixed_quad(). Suitable for smooth integrands where the number of evaluation points is known in advance.

Parameters:

f (Callable) – Integrand: f(x) -> float.
a (float) – Lower integration limit (finite).
b (float) – Upper integration limit (finite).
n (int) – Number of Gaussian quadrature points. Defaults to 10.

Returns:

Integral estimate. Error is set to 0.0 since no: error estimate is available from the fixed rule.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If n < 1 or limits are non-finite.

integrationdifferentiation.integrate_2d(f, ax, bx, g, h, abs_tol=1.49e-08)¶

2D iterated integral over a variable-limit domain (IMSL TWODQ).

Computes:

∫_{ax}^{bx} ∫_{g(x)}^{h(x)} f(y, x) dy dx

via scipy.integrate.dblquad().

Parameters:

f (Callable) – Integrand f(y, x) -> float. Note the argument order: inner variable first.
ax (float) – Lower limit for the outer integral in x.
bx (float) – Upper limit for the outer integral in x.
g (Callable) – Lower limit function for the inner integral: g(x) -> float.
h (Callable) – Upper limit function for the inner integral: h(x) -> float.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.

Returns:

Integral estimate with error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f, g, or h is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_2d_singular(f, ax, bx, ay, by, abs_tol=1.49e-08)¶

2D quadrature over a rectangular domain with possible singularity (IMSL QDAG2D).

Computes ∫_{ax}^{bx} ∫_{ay}^{by} f(y, x) dy dx via scipy.integrate.dblquad() with constant inner limits.

Parameters:

f (Callable) – Integrand f(y, x) -> float.
ax (float) – Lower outer limit in x.
bx (float) – Upper outer limit in x.
ay (float) – Lower inner limit in y.
by (float) – Upper inner limit in y.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.

Returns:

Integral estimate with error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_3d_singular(f, ax, bx, ay, by, az, bz, abs_tol=1.49e-08)¶

3D quadrature over a rectangular domain with possible singularity (IMSL QDAG3D).

Computes ∫_{ax}^{bx} ∫_{ay}^{by} ∫_{az}^{bz} f(z, y, x) dz dy dx via scipy.integrate.tplquad().

Parameters:

f (Callable) – Integrand f(z, y, x) -> float.
ax (float) – Lower outer limit in x.
bx (float) – Upper outer limit in x.
ay (float) – Lower limit in y.
by (float) – Upper limit in y.
az (float) – Lower inner limit in z.
bz (float) – Upper inner limit in z.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.

Returns:

Integral estimate with error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If abs_tol is not positive.

integrationdifferentiation.integrate_nd(f, ranges, abs_tol=1.49e-08)¶

N-dimensional quadrature over a hyperrectangle (IMSL QAND).

Computes the integral of f(x_1, ..., x_d) over the product domain ∏_i [ranges[i][0], ranges[i][1]] using scipy.integrate.nquad().

Parameters:

f (Callable) – Integrand: f(*args) -> float where args has len(ranges) elements.
ranges (list[tuple[float, float]]) – Integration bounds for each dimension as (lower, upper) pairs. Innermost dimension first.
abs_tol (float) – Absolute error tolerance. Defaults to 1.49e-8.

Returns:

Integral estimate with error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If ranges is empty or abs_tol is not positive.

integrationdifferentiation.integrate_monte_carlo(f, ranges, n_points=10000, seed=None)¶

Quasi-Monte Carlo integration over a hyperrectangle (IMSL QMC).

Uses a Halton quasi-random sequence via scipy.integrate.qmc_quad() for low-discrepancy sampling, which typically converges faster than plain Monte Carlo for smooth integrands.

Parameters:

f (Callable) – Integrand that accepts a 1D array of length len(ranges) and returns a scalar, i.e. f(x) -> float where x.shape == (d,).
ranges (list[tuple[float, float]]) – Integration bounds for each dimension as (lower, upper) pairs.
n_points (int) – Number of quasi-random sample points per estimate. Defaults to 10000.
seed (int | None) – Optional random seed for reproducibility. Defaults to None.

Returns:

Integral estimate with standard error bound.

Return type:

IntegrationResult

Raises:

TypeError – If f is not callable.
ValueError – If ranges is empty or n_points < 1.

integrationdifferentiation.gauss_quadrature_rule(n, weight_type='legendre')¶

Gauss quadrature rule for classical weight functions (IMSL GQRUL).

Returns the n-point Gaussian quadrature abscissas and weights for the specified classical weight function using scipy.special root-finding routines.

Supported weight types:

'legendre' — Gauss-Legendre on [-1, 1], weight 1
'chebyshev' — Gauss-Chebyshev (Type 1) on [-1, 1], weight (1-x²)^(-1/2)
'laguerre' — Gauss-Laguerre on [0, ∞), weight exp(-x)
'hermite' — Gauss-Hermite on (-∞, ∞), weight exp(-x²)
'jacobi' — Gauss-Jacobi on [-1, 1] with parameters alpha=0, beta=0 (equivalent to Legendre)

Parameters:

n (int) – Number of quadrature points. Must be >= 1.
weight_type (str) – Classical weight function identifier. Defaults to 'legendre'.

Returns:

Rule with points sorted in ascending order and: corresponding weights.

Return type:

GaussQuadratureRule

Raises:

ValueError – If n < 1 or weight_type is not one of the supported strings.

integrationdifferentiation.gauss_quadrature_from_recurrence(alpha, beta)¶

Gauss quadrature rule from three-term recurrence coefficients (IMSL GQRCF).

Constructs the quadrature rule using the Golub–Welsch algorithm: builds the symmetric tridiagonal Jacobi matrix from the recurrence coefficients, computes its eigendecomposition via numpy.linalg.eigh(), and derives the weights from the first components of the eigenvectors.

The recurrence relation is:

p_{k+1}(x) = (x - alpha_k) p_k(x) - beta_k p_{k-1}(x)

with p_{-1} = 0, p_0 = 1, and beta_0 = mu_0 (total weight integral).

Parameters:

alpha (np.ndarray) – Diagonal recurrence coefficients, shape (n,).
beta (np.ndarray) – Off-diagonal recurrence coefficients, shape (n,). beta[0] is the total weight mu_0; beta[k] for k >= 1 satisfies beta[k] = ||p_k||^2 / ||p_{k-1}||^2.

Returns:

n-point quadrature rule with nodes and weights.

Return type:

GaussQuadratureRule

Raises:

ValueError – If alpha and beta do not have the same length, or if n < 1.

integrationdifferentiation.recurrence_coefficients_classical(n, weight_type='legendre', a=0.0, b=1.0)¶

Three-term recurrence coefficients for classical weight functions (IMSL RECCF).

Returns arrays (alpha, beta) of length n for the recurrence

p_{k+1}(x) = (x - alpha_k) p_k(x) - beta_k p_{k-1}(x)

Supported weight types:

'legendre' — on [-1, 1]; a, b unused
'chebyshev' — Chebyshev T on [-1, 1]; a, b unused
'laguerre' — on [0, ∞); a unused, b is the shape parameter (default b=1 → standard Laguerre)
'hermite' — on (-∞, ∞); a, b unused
'jacobi' — on [-1, 1]; a = Jacobi alpha, b = Jacobi beta

Parameters:

n (int) – Number of recurrence coefficients. Must be >= 1.
weight_type (str) – Classical weight function. Defaults to 'legendre'.
a (float) – Parameter a (Jacobi alpha or Laguerre shape). Defaults to 0.0.
b (float) – Parameter b (Jacobi beta). Defaults to 1.0.

Returns:

(alpha, beta) arrays of length n.: beta[0] equals the total integral of the weight function mu_0; beta[k] for k >= 1 equals ||p_k||^2/||p_{k-1}||^2.

Return type:

tuple[np.ndarray, np.ndarray]

Raises:

ValueError – If n < 1 or weight_type is not supported.

integrationdifferentiation.recurrence_from_quadrature(points, weights)¶

Recurrence coefficients from a discrete quadrature rule (IMSL RECQR).

Uses the Stieltjes–Lanczos three-term recurrence procedure to compute (alpha, beta) from the given n-point quadrature rule. The computed coefficients satisfy:

p_{k+1}(x) = (x - alpha_k) p_k(x) - beta_k p_{k-1}(x)

with beta[0] = sum(weights) (total mass mu_0).

Parameters:

points (np.ndarray) – Quadrature abscissas, shape (n,).
weights (np.ndarray) – Corresponding positive quadrature weights, shape (n,).

Returns:

(alpha, beta) each of length n.

Return type:

tuple[np.ndarray, np.ndarray]

Raises:

ValueError – If points and weights have different lengths or are empty, or if any weight is non-positive.

integrationdifferentiation.fejer_quadrature_rule(n)¶

Fejer Type-1 quadrature rule (IMSL FQRUL).

Computes the n-point Fejér Type-1 rule on [-1, 1] with nodes at

x_k = cos((2k - 1) * π / (2n)),  k = 1, ..., n

and weights computed via the trigonometric sum formula. All weights are positive and the rule integrates polynomials of degree up to n - 1 exactly.

Parameters:: n (int) – Number of quadrature points. Must be >= 1.
Returns:: Rule with nodes and weights in ascending order.
Return type:: GaussQuadratureRule
Raises:: ValueError – If n < 1.

integrationdifferentiation.differentiate(f, x, order=1, h=None)¶

Approximate derivative using Richardson-extrapolated central differences (IMSL DERIV).

Computes the 1st, 2nd, or 3rd derivative of f at the point x using central-difference formulas combined with Richardson extrapolation to achieve \(O(h^4)\) accuracy.

Formulas:

order=1: D1(h) = (f(x+h) - f(x-h)) / (2h); Richardson: (4*D1(h/2) - D1(h)) / 3
order=2: D2(h) = (f(x+h) - 2f(x) + f(x-h)) / h²; Richardson: (4*D2(h/2) - D2(h)) / 3
order=3: D3(h) = (f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2h)) / (2h³); Richardson: (4*D3(h/2) - D3(h)) / 3

Parameters:

f (Callable) – Function to differentiate: f(x) -> float.
x (float) – Point at which to evaluate the derivative.
order (int) – Derivative order: 1, 2, or 3. Defaults to 1.
h (float | None) – Step size. If None, an adaptive step max(abs(x), 1.0) * 1e-4 is used. Defaults to None.

Returns:

Approximation to the order-th derivative of f at x.

Return type:

float

Raises:

TypeError – If f is not callable.
ValueError – If order is not 1, 2, or 3, or if h is provided but is not positive.

integrationdifferentiation.launch_documentation(opener=<function open>, package_dir=None)¶

Launch package documentation in a web browser.

The launcher opens local HTML documentation only.

Parameters:

opener (Callable[[str], bool]) – Browser opener function. Defaults to webbrowser.open.
package_dir (Path | None) – Optional package directory override used for tests. Defaults to the directory of this module.

Returns:

The file://... URI opened.

Return type:

str

Raises:

FileNotFoundError – If no local documentation index can be located.