quadax.quadts

Global adaptive quadrature using trapezoidal tanh-sinh rule.

Integrate fun from interval[0] to interval[-1] using a h-adaptive scheme with error estimate. Breakpoints can be specified in interval where integration difficulty may occur.

This can often be the most efficient method for smooth integrands or weak endpoint singularities (up to around x**-0.5, including those induced by an algebraically decaying integrand on an infinite interval). Beyond that the truncation floor in the map (limited by working precision) dominates and no amount of refinement can do better. In those cases the extrapolating method quadgk() is the more reliable choice.

Parameters:

fun (callable) – Function to integrate, should have a signature of the form fun(x, *args) -> float, Array. Should be JAX transformable.
interval (array-like) – Lower and upper limits of integration with possible breakpoints. Use np.inf to denote infinite intervals. Its dtype sets the working precision: the integrand is called with an x of this dtype, and the result follows it unless the integrand upcasts. Integer types or python floats fall back to the JAX default. Must be real; complex integrands are supported, complex limits are not.
args (tuple, optional) – Extra arguments passed to fun.
full_output (bool, optional) – If True, return the full state of the integrator. See below for more information.
epsabs (float, optional) – Absolute and relative error tolerance. Default is the square root of the machine precision of the working dtype, ie of interval, or of the integrand’s own dtype if that is the coarser of the two. Algorithm tries to obtain an accuracy of abs(i-result) <= max(epsabs, epsrel*abs(i)) where i = integral of fun over interval, and result is the numerical approximation.
epsrel (float, optional) – Absolute and relative error tolerance. Default is the square root of the machine precision of the working dtype, ie of interval, or of the integrand’s own dtype if that is the coarser of the two. Algorithm tries to obtain an accuracy of abs(i-result) <= max(epsabs, epsrel*abs(i)) where i = integral of fun over interval, and result is the numerical approximation.
max_ninter (int, optional) – An upper bound on the number of sub-intervals used in the adaptive algorithm.
order (int) – Order of local integration rule. Must be odd.
norm (int, callable) – Norm to use for measuring error for vector valued integrands. No effect if the integrand is scalar valued. If an int, uses p-norm of the given order, otherwise should be callable.
extrapolate (bool, optional) – Whether to accelerate convergence by applying Wynn’s epsilon algorithm to the sequence of running totals, off by default. Unlike the other adaptive routines this rarely helps here: the tanh-sinh rule converges doubly exponentially, so the running totals have no geometric tail for the epsilon algorithm to sum. Where a tanh-sinh integration is inaccurate the limit is generally the resolution of the abscissas near the endpoints, which acceleration cannot recover.
adjoint (AbstractAdjoint, optional) – How to compute derivatives of the quadrature. Default is DirectAdjoint(), which gives the exact derivative of the discretized problem, and is the cheaper option for a cheap integrand. LeibnizAdjoint gives the derivative its own error control (ie, can better approximate the true continuous derivative), and is faster when the integrand is expensive or max_ninter is generous; see Adjoints for when that is worth paying for.
batch_size (int, optional) – Maximum number of points at which to evaluate the integrand in parallel. Default is all of the local rule’s nodes at once, which is fastest but makes peak memory scale with the order. Lower it to reduce memory on an expensive integrand. Values larger than the number of nodes are clipped to it.
throw (bool, optional) – Whether to raise an error if the routine does not converge. If True, a run that terminates for any reason other than reaching the requested tolerance raises with the message its status carries. If False, the default, that status is reported on the returned info and left to the caller to act on.

Returns:

y (float, Array) – The integral of fun from a to b.
info (QuadratureInfo) – Named tuple with the following fields:
- err : (float) Estimate of the error in the approximation.
- neval : (int) Total number of function evaluations.
- status : (int) Code for why the routine terminated, one of quadax.STATUS. STATUS.normal (0) means the requested tolerances were reached; every other code names a difficulty, whose message is print(quadax.STATUS[status]). Where a run meets more than one condition the most severe is reported.
- info : (dict or None) Other information returned by the algorithm. Only present if full_output is True. Contains the following:
  - ’ninter’ : (int) The number, K, of sub-intervals produced in the subdivision process.
  - ’a_arr’ : (ndarray) rank-1 array of length max_ninter, the first K elements of which are the left end points of the (remapped) sub-intervals in the partition of the integration range.
  - ’b_arr’ : (ndarray) rank-1 array of length max_ninter, the first K elements of which are the right end points of the (remapped) sub-intervals.
  - ’r_arr’ : (ndarray) rank-1 array of length max_ninter, the first K elements of which are the integral approximations on the sub-intervals.
  - ’e_arr’ : (ndarray) rank-1 array of length max_ninter, the first K elements of which are the moduli of the absolute error estimates on the sub-intervals.

Notes

Adaptive algorithms are inherently somewhat sequential, so perfect parallelism is generally not achievable. The local quadrature rule evaluates the integrand at all of its nodes at once, so using higher order methods will generally be more efficient on GPU/TPU. batch_size splits that evaluation up where the memory it needs is the binding constraint instead.