API reference

quadgk

The most commonly used function from the QuadGK package is the quadgk function, used to compute numerical integrals (by h-adaptive Gauss–Kronrod quadrature):

QuadGK.quadgkFunction
quadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm, segbuf=nothing, eval_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.

Note

It is advisable to increase the integration order in rough proportion to the precision, for smooth integrands.

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

Instead of passing the integration domain as separate arguments a,b,c..., you can alternatively pass the domain as a single argument: an array of endpoints [a,b,c...] or an array of interval tuples [(a,b), (b,c)]. (The latter enables you to integrate over a disjoint domain if you want.)

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

Note

Only one-dimensional integrals are provided by this function. For multi-dimensional integration (cubature), there are many different algorithms (often much better than simple nested 1d integrals) and the optimal choice tends to be very problem-dependent. See the Julia external-package listing for available algorithms for multidimensional integration or other specialized tasks (such as integrals of highly oscillatory or singular functions).

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. Alternatively, one can replace the first quadgk(...) call with quadgk_segbuf(...) to return the segment buffer from a given call. This buffer can be used across multiple calls to avoid repeated allocation. Upon return from quadgk, the segbuf array contains an array of subintervals that were used for the final quadrature evaluation.

By passing eval_segbuf=segbuf to a subsequent call to quadgk, these subintervals can be re-used as the starting point for the next integrand evaluation (over the same domain), even for an integrand of a different result type; by also passing maxevals=0, further refinement of these subintervals is prohibited, so that it forces the same quadrature rule to be used (which is useful for evaluating e.g. derivatives of the approximate integral).

source
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.

source

The quadgk function also has variants quadgk_count (which also returns a count of the integrand evaluations), quadgk_print (which also prints each integrand evaluation), quadgk! (which implements an in-place API for array-valued functions).

QuadGK.quadgk_countFunction
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.

source
QuadGK.quadgk_printFunction
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.

source
QuadGK.quadgk!Function
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.

source

For a vectorized integrand that can evaluate the integrand at multiple points simultaneously, quadgk accepts a BatchIntegrand wrapper around the user's integrand and pre-allocated input and output buffers.

QuadGK.BatchIntegrandType
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.

source

The quadgk functions accept a segbuf=segbuf keyword to re-use internal buffers from one call to the next, and an eval_segbuf=segbuf to re-use the adaptive subdivision from one call to the next. To return the internal segbuf from an integration in order to re-use it, you can use the quadgk_segbuf function and its variants. You can also use the alloc_segbuf function to pre-allocate a segment buffer if you know the correct numeric types for the domain endpoints, integral, and error estimate:

QuadGK.quadgk_segbufFunction
quadgk_segbuf(args...; kws...)

Identical to quadgk(args...; kws...), but returns a tuple (I, E, segbuf) where segbuf is a segment buffer (storing the subintervals used for the final integral evaluation) that can be passed as a segbuf and/or eval_segbuf argument on subsequent calls to quadgk and related functions.

source
QuadGK.quadgk_segbuf_countFunction
quadgk_segbuf_count(args...; kws...)

Identical to quadgk_count(args...; kws...), but returns a tuple (I, E, segbuf, count) where segbuf is a segment buffer (storing the subintervals used for the final integral evaluation) that can be passed as a segbuf and/or eval_segbuf argument on subsequent calls to quadgk and related functions.

source
QuadGK.quadgk_segbuf_printFunction
quadgk_segbuf_print(args...; kws...)

Identical to quadgk_print(args...; kws...), but returns a tuple (I, E, segbuf, count) where segbuf is a segment buffer (storing the subintervals used for the final integral evaluation) that can be passed as a segbuf and/or eval_segbuf argument on subsequent calls to quadgk and related functions.

source
QuadGK.quadgk_segbuf!Function
quadgk_segbuf!(args...; kws...)

Identical to quadgk!(args...; kws...), but returns a tuple (I, E, segbuf) where segbuf is a segment buffer (storing the subintervals used for the final integral evaluation) that can be passed as a segbuf and/or eval_segbuf argument on subsequent calls to quadgk and related functions.

source
QuadGK.alloc_segbufFunction
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.

Alternatively, you can replace your first call to quadgk(...) with a call toquadgksegbuf(...), which returns the computed segment buffer from your first integration. This saves you the trouble of figuring outdomaintype` etc., which may not be obvious if the integrand is a variable.

source

Gauss and Gauss–Kronrod rules

For more specialized applications, you may wish to construct your own Gauss or Gauss–Kronrod quadrature rules, as described in Gauss and Gauss–Kronrod quadrature rules. To compute rules for $\int_{-1}^{+1} f(x) dx$ and $\int_a^b f(x) dx$ (unweighted integrals), use:

QuadGK.gaussMethod
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.

source
QuadGK.kronrodMethod
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 wg (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
source

More generally, to compute rules for $\int_a^b W(x) f(x) dx$ (weighted integrals, as described in Gaussian quadrature and arbitrary weight functions), use the following methods if you know the Jacobi matrix for the orthogonal polynomials associated with your weight function:

QuadGK.gaussMethod
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]$.

source
QuadGK.kronrodMethod
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, wg) of $n$ quadrature points x[i] and weights w[i] to integrate functions, i.e. sum(w .* f.(x)) approximates the integral $I[f]$. wg are the weights of the embedded Gauss rule corresponding to the points x[2:2:end], which can be used for error estimation.

source
QuadGK.HollowSymTridiagonalType
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.

source

Most generally, if you know only the weight function $W(x)$ and the interval $(a,b)$, you can construct Gauss and Gauss–Kronrod rules completely numerically using:

QuadGK.gaussMethod
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.

source
QuadGK.kronrodMethod
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.

source