API Documentation

Adaptive integration of a callable function or method

quadgk(fun, interval[, args, full_output, ...])	Global adaptive quadrature using Gauss-Kronrod rule.
quadcc(fun, interval[, args, full_output, ...])	Global adaptive quadrature using Clenshaw-Curtis rule.
quadts(fun, interval[, args, full_output, ...])	Global adaptive quadrature using trapezoidal tanh-sinh rule.
romberg(fun, interval[, args, full_output, ...])	Romberg integration of a callable function or method.
rombergts(fun, interval[, args, ...])	Romberg integration with tanh-sinh (aka double exponential) transformation.

Quadrature Rules

AbstractQuadratureRule()	Abstract base class for 1D quadrature rules.
GaussKronrodRule([order, norm, batch_size])	Integrate a function from a to b using a fixed order Gauss-Kronrod rule.
ClenshawCurtisRule([order, norm, batch_size])	Integrate a function from a to b using a fixed order Clenshaw-Curtis rule.
TanhSinhRule([order, norm, batch_size])	Integrate a function from a to b using a fixed order Tanh-Sinh trapezoidal rule.

Adjoints

Adjoints control how derivatives of a quadrature are computed, without changing what the quadrature itself returns. Pass one as the adjoint argument.

AbstractAdjoint()	Abstract base class for adjoint methods.
DirectAdjoint([checkpoint, chunk_size])	Differentiate the quadrature exactly, reusing the converged subdivision.
LeibnizAdjoint()	Differentiate by the Leibniz rule, either mode, with its own error control.

Choosing an adjoint

There are two, both supporting forward and reverse mode.

DirectAdjoint is the default. It differentiates the discretization the primal solve settled on, so the derivative costs no error control of its own, and for a cheap integrand it is usually the cheaper option in either mode.

LeibnizAdjoint instead evaluates the derivative with a second adaptive solve, giving it its own error control rather than inheriting the subdivision chosen for the integral. This buys:

• Accuracy. When the derivative of the integrand is sharply peaked somewhere the integrand itself is smooth, the subdivision that resolves the integral need not resolve its derivative. On such a problem at a loose tolerance the difference can be several orders of magnitude; at a tolerance tight enough that the integral’s subdivision resolves the derivative anyway, it may be less of an issue.

• Speed when the integrand is expensive. The derivative solve stops as soon as the derivative has converged, rather than covering the subdivision the integral needed, so the more one evaluation of the integrand costs, the more there is to save. This applies to both forward and reverse modes.

Against that, its reverse pass carries the workspace of the second solve, so on a scalar integrand it generally costs more memory than the default; on a vector or matrix valued integrand, where the stored subdivision is what dominates, it can cost less.

Both pick up the jump term from differentiating with respect to a breakpoint that sits on a discontinuity. Neither can see a discontinuity that has no breakpoint at it - declare one there.

romberg() and rombergts() have no subdivision to reuse - DirectAdjoint freezes the number of Richardson levels instead - and there the two cost about the same, so the choice is about accuracy alone.

These are rules of thumb, not laws. The balance shifts with how expensive the integrand is relative to the quadrature around it, how hard its derivative is to integrate compared to the integrand, and how many parameters are involved. Time both on your own problem before caring much about the difference.

Both take two options controlling the memory a derivative needs. checkpoint (off by default) recomputes each block of sub-intervals during the backward pass instead of storing it, and is where nearly all of the saving is for DirectAdjoint, whose backward pass replays the frozen subdivision. LeibnizAdjoint solves for the adjoint instead of replaying anything, so there both options touch only the fixed-mesh pass that supplies the derivative with respect to the interval, and checkpoint does nothing at all unless the limits are being differentiated. chunk_size sets how many sub-intervals of the frozen subdivision are evaluated at once, and multiplies with the batch_size of the routine itself: a gradient evaluates the integrand at up to chunk_size * batch_size points at a time.

Integrating function from sampled values

trapezoid(y, *[, x, dx, axis])	Integrate along the given axis using the composite trapezoidal rule.
cumulative_trapezoid(y, *[, x, dx, axis, ...])	Cumulatively integrate y(x) using the composite trapezoidal rule.
simpson(y, *[, x, dx, axis])	Integrate y(x) from samples using the composite Simpson's rule.
cumulative_simpson(y, *[, x, dx, axis, initial])	Cumulatively integrate y(x) using the composite Simpson's 1/3 rule.

Low level routines and wrappers

adaptive_quadrature(rule, fun, interval[, ...])

Global adaptive quadrature.