quadax.quadcc

Global adaptive quadrature using Clenshaw-Curtis 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.

A good general purpose integrator for most reasonably well behaved functions over finite or infinite intervals, and a reasonable alternative to quadgk(). It’s main advantage is in allowing arbitrary high orders, which can be useful for smooth but highly oscillatory integrands in the absence of a specialized solver, in which case choosing order to have ~7-8 points per period is often the most efficient.

As with quadgk(), an interior jump or singularity is best passed as a breakpoint in interval rather than left for the subdivision to find.

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 a multiple of 4, or with closed=False any even order of at least 4; see ClenshawCurtisRule.
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.
closed (bool, optional) – Whether the interval endpoints are among the nodes of the local rule. The default closed rule is cheaper on smooth, peaked and oscillatory integrands. The open (Fejer-2) rule never evaluates the integrand at an interval endpoint, which is what to use for integrands that are singular or undefined there; it is markedly cheaper on infinite intervals whose integrand decays algebraically, and on endpoint singularities. See ClenshawCurtisRule.
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.