API reference

quadgk(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.)

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 quadgk(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, quadgk(...) 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.

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

Like quadgk, 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 changes how the adaptive refinement is done by batching multiple segments together. Choosing max_batch<=4*order+2 will reproduce the result of non-batched quadgk, however if max_batch=n*(4*order+2) up to 2n Kronrod rules will be evaluated together, which can produce slightly different results and sometimes require more integrand evaluations when using relative tolerances.

quadgk_count(f, args...; kws...)

Identical to quadgk but returns a triple (I, E, count) of the estimated integral I, the estimated error bound E, and a count of the number of times the integrand f was evaluated.

The count of integrand evaluations is a useful performance metric: a large number typically indicates a badly behaved integrand (with singularities, discontinuities, sharp peaks, and/or rapid oscillations), in which case it may be possible to mathematically transform the problem in some way to improve the convergence rate.

quadgk_print([io], f, args...; kws...)

Identical to quadgk, but prints each integrand evaluation to the stream io (defaults to stdout) in the form:

f(x1) = y1
f(x2) = y2
...

which is useful for pedagogy and debugging.

Also, like quadgk_count, it returns a triple (I, E, count) of the estimated integral I, the estimated error bound E, and a count of the number of times the integrand f was evaluated.

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

Like quadgk, 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 quadgk, the return value is a tuple (I,E) of the estimated integral I and the estimated error E. However, in quadgk! the estimated integral is written in-place into the result argument, so that I === result.

Otherwise, the behavior is identical to quadgk.

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.

function alloc_segbuf(domain_type=Float64, range_type=Float64, error_type=Float64; size=1)

This helper will allocate a segment buffer for segments to a quadgk(...) call with the given domain_type, which is the same as the type of the integration limits, range_type i.e. the range of the function being integrated and error_type, the type returned by the norm given to quadgk(...) and starting with the given size. The buffer can then be reused across multiple compatible calls to quadgk(...) to avoid repeated allocation.

BatchIntegrand(f!, y::AbstractVector, [x::AbstractVector]; max_batch=typemax(Int))
BatchIntegrand{Y,X}(f!; max_batch=typemax(Int)) where {Y,X}
BatchIntegrand{Y}(f!; max_batch=typemax(Int)) where {Y}

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 arguments y, and optionally x, are pre-allocated buffers of the correct element type for the domain and range of f!. Additionally, they must be resize!-able since the number of evaluation points may vary between calls to f!. Alternatively, the element types of those buffers, Y and optionally X, may be passed as type parameters to the constructor. The max_batch keyword roughly limits the number of nodes passed to the integrand, though at least 4*order+2 nodes will be used by the GK rule.

Internal allocations

If x or X are not specified, quadgk internally creates a new BatchIntegrand with the user-supplied y buffer and a freshly-allocated x buffer based on the domain types. So, if you want to reuse the x buffer between calls, supply {Y,X} or pass y,x explicitly.

gauss([T,] n)
gauss([T,] n, a, b)

Return a pair (x, w) of n quadrature points x[i] and weights w[i] to integrate functions on the interval $(a, b)$, which defaults to $(-1,1)$, i.e. sum(w .* f.(x)) approximates the integral $\int_a^b f(x) dx$.

Uses the Golub–Welch method described in Trefethen & Bau, Numerical Linear Algebra, to find the n-point Gaussian quadrature rule in O(n²) operations.

T is an optional parameter specifying the floating-point type, defaulting to Float64. Arbitrary precision (BigFloat) is also supported. If T is not supplied, but the interval (a, b) is passed, then the floating-point type is determined from the types of a and b.

kronrod([T,] n)
kronrod([T,] n, a, b)

Compute $2n+1$ Kronrod points x[i] and weights w[i] based on the description in Laurie (1997), appendix A, for integrating on the interval $(a,b)$ (defaulting to $[-1,1]$).

If a and b are not passed, since the rule is symmetric, this only returns the n+1 points with x <= 0. The function Also computes the embedded n-point Gauss quadrature weights gw (again for x <= 0 if a and b are not passed), corresponding to the points x[2:2:end]. Returns (x,w,wg) in O(n²) operations.

T is an optional parameter specifying the floating-point type, defaulting to Float64. Arbitrary precision (BigFloat) is also supported.

Given these points and weights, the estimated integral I and error E can be computed for an integrand f(x) as follows:

x, w, wg = kronrod(n)
fx⁰ = f(x[end])                # f(0)
x⁻ = x[1:end-1]                # the x < 0 Kronrod points
fx = f.(x⁻) .+ f.((-).(x⁻))    # f(x < 0) + f(x > 0)
I = sum(fx .* w[1:end-1]) + fx⁰ * w[end]
if isodd(n)
    E = abs(sum(fx[2:2:end] .* wg[1:end-1]) + fx⁰*wg[end] - I)
else
    E = abs(sum(fx[2:2:end] .* wg[1:end])- I)
end

gauss(J::AbstractMatrix, unitintegral::Real=1, [ (a₀,b₀) => (a,b) ])

Construct the $n$-point Gaussian quadrature rule for $I[f] = \int_a^b w(x) f(x) dx$ from the $n \times n$ symmetric tridiagonal Jacobi matrix J corresponding to the orthogonal polynomials for that weighted integral. The value of unitintegral should be $I[1]$, the integral of the weight function.

An optional argument (a₀,b₀) => (a,b) allows you to specify that J was originally defined for a different interval $(a_0, b_0)$, which you want to rescale to a given $(a, b)$. (gauss will rescale the points and weights for you.)

Returns a pair (x, w) of $n$ quadrature points x[i] and weights w[i] to integrate functions, i.e. sum(w .* f.(x)) approximates the integral $I[f]$.

kronrod(J::AbstractMatrix, n::Integer, unitintegral::Real=1, [ (a₀,b₀) => (a,b) ])

Construct the $2n+1$-point Gauss–Kronrod quadrature rule for $I[f] = \int_a^b w(x) f(x) dx$ from the $m \times m$ symmetric tridiagonal Jacobi matrix J corresponding to the orthogonal polynomials for that weighted integral, where m ≥ (3n+3)÷2. The value of unitintegral should be $I[1]$, the integral of the weight function.

An optional argument (a₀,b₀) => (a,b) allows you to specify that J was originally defined for a different interval $(a_0, b_0)$, which you want to rescale to a given $(a, b)$. (gauss will rescale the points and weights for you.)

Returns a tuple (x, w, gw) of $n$ quadrature points x[i] and weights w[i] to integrate functions, i.e. sum(w .* f.(x)) approximates the integral $I[f]$. gw are the weights of the embedded Gauss rule corresponding to the points x[2:2:end], which can be used for error estimation.

QuadGK.HollowSymTridiagonal(ev::AbstractVector)

Construct a "hollow" symmetric tridiagonal matrix, whose diagonal entries are zero and whose first sub/super-diagonal is ev.

The HollowSymTridiagonal type can be passed to gauss or kronrod for Jacobi matrices to dispatch to specialized methods that exploit the special "hollow" structure arising for symmetric weight functions, in order to generate symmetric quadrature rules more efficiently.

gauss(W, N, a, b; rtol=sqrt(eps), quad=quadgk)

Return a pair (x, w) of N quadrature points x[i] and weights w[i] to integrate functions on the interval (a, b) multiplied by the weight function W(x). That is, sum(w .* f.(x)) approximates the integral ∫ W(x)f(x)dx from a to b.

This function performs 2N numerical integrals of polynomials against W(x) using the integration function quad (defaults to quadgk) with relative tolerance rtol (which defaults to half of the precision eps of the endpoints). This is followed by an O(N²) calculations. So, using a large order N is expensive.

If W has lots of singularities that make it hard to integrate numerically, you may need to decrease rtol. You can also pass in a specialized quadrature routine via the quad keyword argument, which should accept arguments quad(f,a,b,rtol=_,atol=_) similar to quadgk. (This is useful if your weight function has discontinuities, in which case you might want to break up the integration interval at the discontinuities.)

The precision of the calculations and return value is determined from the types of a and b.

kronrod(W, N, a, b; rtol=sqrt(eps), quad=quadgk)

Return a tuple (x, w, wg) of N quadrature points x[i] and weights w[i] to integrate functions on the interval (a, b) multiplied by the weight function W(x), along with the weights wg of an embedded Gauss rule corresponding to x[2:2:end], similar to the gauss(W, N, a, b) function and analogous to kronrod(N) (which only returns the x ≤ 0 points for a constant weight function).

That is, I = sum(w .* f.(x)) approximates the integral ∫ W(x)f(x)dx from a to b. And an error estimate is abs(I - Ig), where Ig is the result Ig = sum(wg .* f.(x[2:2:end])) of the embedded Gauss rule.

This function performs ≈ 3N+3 numerical integrals of polynomials against W(x) using the integration function quad (defaults to quadgk) with relative tolerance rtol (which defaults to half of the precision eps of the endpoints). This is followed by an O(N²) calculations. So, using a large order N is expensive.