quadax.tanhsinh

Tanh-sinh (aka double exponential) quadrature on a uniformly refined mesh.

Returns the integral of fun (a function of one variable) over interval.

Substitutes x = tanh(pi/2 sinh(t)), which flattens an endpoint singularity into a doubly exponentially decaying tail, and applies the trapezoidal rule to the result on a mesh that is halved until the requested tolerance is met.

Performs well for functions with mild singularities at the endpoints or integration over infinite intervals.

Interior breaks are the one thing it cannot do at all: its mesh is uniform, it accepts no breakpoints, and a singularity or a jump away from the limits should go to quadgk() with that point passed in interval.

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. 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. Breakpoints are not supported; pass exactly two limits.
args (tuple, optional) – Extra arguments passed to fun.
full_output (bool, optional) – If True, return the estimate from every refinement level. See below.
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.
divmax (int, optional) – Maximum number of divisions, ie the most times the mesh may be halved.
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.
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 divmax 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. Defaults to 2**divmin, which is one batch for the whole starting grid and exactly one for the refinement level after it. Each refinement level doubles the number of new points, so raising this together with divmin is usually worth a lot on GPU/TPU, at the cost of peak memory scaling with it. Levels with fewer new points than one batch are padded up to a full batch, so a level costs batch_size evaluations however few points it places; that padding is what keeps a single batch shape traced for every level rather than one per level. Clipped to the largest number of points any one level places.
divmin (int, optional) – Number of refinement levels the first pass places at once, rather than arriving at one level at a time. Costs nothing over starting from a single interval: the coarser levels are filled from the same evaluations.
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) – Approximation to the integral
info (QuadratureInfo) – Named tuple with the following fields:
- err : (float) Estimate of the error in the approximation. Built from how far the estimate moved over the last few refinement levels rather than over the last one alone, plus the tail that movement’s own contraction rate implies, the mass the map leaves outside the range integrated over, and floored at the precision the integrand can be summed to.
- 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 : (ndarray or None) Only present if full_output is True: the trapezoidal estimate at each refinement level, of shape (divmax + 1, ...). Levels beyond the one the routine stopped at are zero.

Notes

The number of new points a refinement level places is only known at run time, so there is no single shape to vectorize over. Integrand evaluations are made in fixed size batches of batch_size, defaulting to 16; raise it to get more parallelism on GPU/TPU.