Package for auxiliary integration, i.e. integrating multiple functions at the same time while ensuring that each converges to its own tolerance. This has a few advantages over vector-valued integrands with custom norms in that the errors from different integrands can be treated separated and the adaptive algorithm can decide which integrand to prioritize based on whether others have already converged. This results in conceptually simpler algorithms, especially when the various integrands may differ in order of magnitude.

This module heavily reuses the source code of QuadGK.jl

Statement of need

Calculating integrals of the form a^2/(f(x)^2+a^2)^2 is challenging in the a -> 0 limit because they become extremely localized while also having vanishingly small tails. I.e. the tails are O(a^2) however the integral is O(a^-1). Thus, setting an absolute tolerance is impractical, since the tolerance also needs to be O(a^2) to resolve the tails (otherwise the peaks will be missed) and that may be asking for many more digits than desired. Another option is to specify a relative tolerance, but a failure mode is that if there is more than one peak to integrate, the algorithm may only resolve the first one because the errors in the tails to find the other peaks become eclipsed by the first peak error magnitudes. When the peak positions are known a priori, the convential solution is to pass in several breakpoints to the integration interval so that each interval has at most one peak, but often finding breakpoints can be an expensive precomputation that is better avoided. Instead, an integrand related to the original may more reliably find the peaks without requiring excessive integrand evaluations or being expensive to compute. Returning to the original example, an ideal auxiliary integrand would be 1/(f(x)+im*a)^2, which has O(1) tails and a O(1) integral. Thus the tails will be resolved in order to find the peaks, which don't need to be resolved to many digits of accuracy. However, since one wants to find the original integral to a certain number of digits, it may be necessary to adapt further after the auxiliary integrand has converged. This is the problem the package aims to solve.

Example

f(x)    = sin(x)/(cos(x)+im*1e-5)   # peaked "nice" integrand
imf(x)  = imag(f(x))                # peaked difficult integrand
f2(x)   = f(x)^2                    # even more peaked
imf2(x) = imf(x)^2                  # even more peaked!

x0 = 0.1    # arbitrary offset of between peak

function integrand(x)
    re, im = reim(f2(x) + f2(x-x0))
    AuxValue(imf2(x) + imf2(x-x0), re)
end

using QuadGK    # plain adaptive integration

quadgk(x -> imf2(x) + imf2(x-x0), 0, 2pi, atol = 1e-5)   # 1.4271103714584847e-7
quadgk(x -> imf2(x) + imf2(x-x0), 0, 2pi, rtol = 1e-5)   # 235619.45750214785

quadgk(x -> imf2(x), 0, 2pi, rtol = 1e-5)   # 78539.81901117883

quadgk(x -> imf2(x-x0), 0, 2pi, rtol = 1e-5)   # 157079.63263294287

using AuxQuadGK # auxiliary integration

auxquadgk(integrand, 0, 2pi, atol=1e-2) # 628318.5306881254
auxquadgk(integrand, 0, 2pi, rtol=1e-2) # 628318.5306867635

As can be seen from the example, plain integration can completely fail to capture the integral despite using stringent tolerances. With a well-chosen auxiliary integrand, often arising naturally from the structure of the integrand, the integration is much more robust to error because it can resolve the regions of interest with the more-easily adaptively integrable problem.

BatchIntegrand(f!, y::Type, x::Type=Nothing; max_batch=typemax(Int))

Constructor for a BatchIntegrand whose range type is known. The domain type is optional. Array buffers for those types are allocated internally.

BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int))

Constructor for a BatchIntegrand accepting an integrand of the form f!(y,x) = y .= f.(x) that can evaluate the integrand at multiple quadrature nodes using, for example, threads, the GPU, or distributed-memory. The max_batch keyword limits the number of nodes passed to the integrand, and it must be at least 4*order+2 to evaluate two GK rules simultaneously. The buffers y,x must both be resize!-able since the number of evaluation points may vary between calls to f!.

BatchIntegrand(f!, y, x; max_batch=typemax(Int))

Constructor for a BatchIntegrand with pre-allocated buffers.

auxquadgk!(f!, result, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm)

Like auxquadgk, but make use of in-place operations for array-valued integrands (or other mutable types supporting in-place operations). In particular, there are two differences from quadgk:

1. The function f! should be of the form f!(y, x) = y .= f(x). That is, it writes the return value of the integand f(x) in-place into its first argument y. (The return value of f! is ignored.)

2. Like auxquadgk, the return value is a tuple (I,E) of the estimated integral I and the estimated error E. However, in auxquadgk! the estimated integral is written in-place into the result argument, so that I === result.

Otherwise, the behavior is identical to auxquadgk.

For integrands whose values are small arrays whose length is known at compile-time, it is usually more efficient to use quadgk and modify your integrand to return an SVector from the StaticArrays.jl package.

auxquadgk!(f::BatchIntegrand, result, a,b,c...; kws...)

Like auxquadgk!, but batches evaluation points for an in-place integrand to evaluate simultaneously. In particular, there are two differences from using quadgk with a BatchIntegrand:

1. f.y must be an array of dimension ndims(result)+1 whose first axes match those of result. The last dimension of y should be reserved for different Kronrod points, and the function f.f! should be of the form f!(y,x) = foreach((v,z) -> v .= f(z), eachslice(y, dims=ndims(y)), x) or

   function f!(y, x) idx = CartesianIndices(axes(y)[begin:end-1]) for (j,i) in zip(axes(y)[end], eachindex(x)) y[idx,j] .= f(x[i]) end end

2. f.y must be resize!-able in the last dimension. Consider using ElasticArrays.jl for this. Otherwise specialize QuadGK.resizelastdim!(A::T, n) for your array type T.

auxquadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm, segbuf=nothing)

Numerically integrate the function f(x) from a to b, and optionally over additional intervals b to c and so on. Keyword options include a relative error tolerance rtol (if atol==0, defaults to sqrt(eps) in the precision of the endpoints), an absolute error tolerance atol (defaults to 0), a maximum number of function evaluations maxevals (defaults to 10^7), and the order of the integration rule (defaults to 7).

Returns a pair (I,E) of the estimated integral I and an estimated upper bound on the absolute error E. If maxevals is not exceeded then E <= max(atol, rtol*norm(I)) will hold. (Note that it is useful to specify a positive atol in cases where norm(I) may be zero.)

Compared to quadgk from QuadGK.jl, auxquadgk implements a few more safeguards for integration of difficult functions. It changes how adaptive refinement is done when using a relative tolerance to refine all segments with an error above the tolerance (instead of just the segment with the largest error). Additionally, if an integrand returns an AuxValue then the heap first integrates the auxiliary value followed by the primary by resorting the heap of segments.

The endpoints a et cetera can also be complex (in which case the integral is performed over straight-line segments in the complex plane). If the endpoints are BigFloat, then the integration will be performed in BigFloat precision as well.

More generally, the precision is set by the precision of the integration endpoints (promoted to floating-point types).

The integrand f(x) can return any numeric scalar, vector, or matrix type, or in fact any type supporting +, -, multiplication by real values, and a norm (i.e., any normed vector space). Alternatively, a different norm can be specified by passing a norm-like function as the norm keyword argument (which defaults to norm).

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each interval is estimated using a Kronrod rule (2*order+1 points) and the error is estimated using an embedded Gauss rule (order points). The interval with the largest error is then subdivided into two intervals and the process is repeated until the desired error tolerance is achieved.

These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. For example, if f has a discontinuity at x=0.7 and you want to integrate from 0 to 1, you should use auxquadgk(f, 0,0.7,1) to subdivide the interval at the point of discontinuity. The integrand is never evaluated exactly at the endpoints of the intervals, so it is possible to integrate functions that diverge at the endpoints as long as the singularity is integrable (for example, a log(x) or 1/sqrt(x) singularity).

For real-valued endpoints, the starting and/or ending points may be infinite. (A coordinate transformation is performed internally to map the infinite interval to a finite one.)

In normal usage, auxquadgk(...) will allocate a buffer for segments. You can instead pass a preallocated buffer allocated using alloc_segbuf(...) as the segbuf argument. This buffer can be used across multiple calls to avoid repeated allocation.

auxquadgk(f::BatchIntegrand, a,b,c...; kws...)

Like auxquadgk, but batches evaluation points for an in-place integrand to evaluate simultaneously. In particular, there are two differences from quadgk

1. The function f.f! should be of the form f!(y, x) = y .= f.(x). That is, it writes the return values of the integand f(x) in-place into its first argument y. (The return value of f! is ignored.) See BatchIntegrand for how to define the integrand.

2. f.max_batch must be large enough to contain 4*order+2 points to evaluate two Kronrod rules simultaneously. Choosing max_batch=4*order+2 will reproduce the result of quadgk, however if max_batch=n*(4*order+2) up to 2n Kronrod rules will be evaluated together, which can produce different results for integrands with multiple peaks when used together with relative tolerances. For an example see the manual

auxquadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm, segbuf=nothing)

Numerically integrate the function f(x) from a to b, and optionally over additional intervals b to c and so on. Keyword options include a relative error tolerance rtol (if atol==0, defaults to sqrt(eps) in the precision of the endpoints), an absolute error tolerance atol (defaults to 0), a maximum number of function evaluations maxevals (defaults to 10^7), and the order of the integration rule (defaults to 7).

Returns a pair (I,E) of the estimated integral I and an estimated upper bound on the absolute error E. If maxevals is not exceeded then E <= max(atol, rtol*norm(I)) will hold. (Note that it is useful to specify a positive atol in cases where norm(I) may be zero.)

Compared to quadgk from QuadGK.jl, auxquadgk implements a few more safeguards for integration of difficult functions. It changes how adaptive refinement is done when using a relative tolerance to refine all segments with an error above the tolerance (instead of just the segment with the largest error). Additionally, if an integrand returns an AuxValue then the heap first integrates the auxiliary value followed by the primary by resorting the heap of segments.

The endpoints a et cetera can also be complex (in which case the integral is performed over straight-line segments in the complex plane). If the endpoints are BigFloat, then the integration will be performed in BigFloat precision as well.

More generally, the precision is set by the precision of the integration endpoints (promoted to floating-point types).

The integrand f(x) can return any numeric scalar, vector, or matrix type, or in fact any type supporting +, -, multiplication by real values, and a norm (i.e., any normed vector space). Alternatively, a different norm can be specified by passing a norm-like function as the norm keyword argument (which defaults to norm).

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each interval is estimated using a Kronrod rule (2*order+1 points) and the error is estimated using an embedded Gauss rule (order points). The interval with the largest error is then subdivided into two intervals and the process is repeated until the desired error tolerance is achieved.

These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. For example, if f has a discontinuity at x=0.7 and you want to integrate from 0 to 1, you should use auxquadgk(f, 0,0.7,1) to subdivide the interval at the point of discontinuity. The integrand is never evaluated exactly at the endpoints of the intervals, so it is possible to integrate functions that diverge at the endpoints as long as the singularity is integrable (for example, a log(x) or 1/sqrt(x) singularity).

For real-valued endpoints, the starting and/or ending points may be infinite. (A coordinate transformation is performed internally to map the infinite interval to a finite one.)

In normal usage, auxquadgk(...) will allocate a buffer for segments. You can instead pass a preallocated buffer allocated using alloc_segbuf(...) as the segbuf argument. This buffer can be used across multiple calls to avoid repeated allocation.

auxquadgk(f::BatchIntegrand, a,b,c...; kws...)

Like auxquadgk, but batches evaluation points for an in-place integrand to evaluate simultaneously. In particular, there are two differences from quadgk

1. The function f.f! should be of the form f!(y, x) = y .= f.(x). That is, it writes the return values of the integand f(x) in-place into its first argument y. (The return value of f! is ignored.) See BatchIntegrand for how to define the integrand.

2. f.max_batch must be large enough to contain 4*order+2 points to evaluate two Kronrod rules simultaneously. Choosing max_batch=4*order+2 will reproduce the result of quadgk, however if max_batch=n*(4*order+2) up to 2n Kronrod rules will be evaluated together, which can produce different results for integrands with multiple peaks when used together with relative tolerances. For an example see the manual

auxquadgk!(f!, result, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm)

Like auxquadgk, but make use of in-place operations for array-valued integrands (or other mutable types supporting in-place operations). In particular, there are two differences from quadgk:

1. The function f! should be of the form f!(y, x) = y .= f(x). That is, it writes the return value of the integand f(x) in-place into its first argument y. (The return value of f! is ignored.)

2. Like auxquadgk, the return value is a tuple (I,E) of the estimated integral I and the estimated error E. However, in auxquadgk! the estimated integral is written in-place into the result argument, so that I === result.

Otherwise, the behavior is identical to auxquadgk.

For integrands whose values are small arrays whose length is known at compile-time, it is usually more efficient to use quadgk and modify your integrand to return an SVector from the StaticArrays.jl package.

auxquadgk!(f::BatchIntegrand, result, a,b,c...; kws...)

Like auxquadgk!, but batches evaluation points for an in-place integrand to evaluate simultaneously. In particular, there are two differences from using quadgk with a BatchIntegrand:

1. f.y must be an array of dimension ndims(result)+1 whose first axes match those of result. The last dimension of y should be reserved for different Kronrod points, and the function f.f! should be of the form f!(y,x) = foreach((v,z) -> v .= f(z), eachslice(y, dims=ndims(y)), x) or

   function f!(y, x) idx = CartesianIndices(axes(y)[begin:end-1]) for (j,i) in zip(axes(y)[end], eachindex(x)) y[idx,j] .= f(x[i]) end end

2. f.y must be resize!-able in the last dimension. Consider using ElasticArrays.jl for this. Otherwise specialize QuadGK.resizelastdim!(A::T, n) for your array type T.