quadax.quadgk

Global adaptive quadrature using Gauss-Kronrod 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.

Basically the same algorithm as scipy.integrate.quad, including the convergence acceleration. The general purpose integrator to reach for first, over finite or infinite intervals. It is generally the most robust and often also the most efficient, on smooth and non-smooth integrands alike.

Where an integrand has a jump or a singularity at a known interior point, passing that point as a breakpoint in interval is worth more than any change of method, since the subdivision no longer has to find it.

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, one of 15, 21, 31, 41, 51, 61.
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, on by default. Not needed for smooth integrands on finite domains, but can help significantly if there are algebraic singularities or infinite intervals. The additional cost is small and constant, so it is only worth switching off for a very cheap integrand where performance is critical.
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.