# 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.quadgk`

— Function`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.

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`

).

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

`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`

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

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_count`

— Function`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.quadgk_print`

— Function`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.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`

:

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

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.BatchIntegrand`

— Type```
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.

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

`QuadGK.quadgk_segbuf_count`

— Function`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.

`QuadGK.quadgk_segbuf_print`

— Function`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.

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

`QuadGK.alloc_segbuf`

— Function`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 to`

quadgk*segbuf(...) , which returns the computed segment buffer from your first integration. This saves you the trouble of figuring outdomain*type` etc., which may not be obvious if the integrand is a variable.

## 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.gauss`

— Method```
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`

.

`QuadGK.kronrod`

— Method```
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
```

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.gauss`

— Method`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]$.

`QuadGK.kronrod`

— Method`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.

`QuadGK.HollowSymTridiagonal`

— Type`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.

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.gauss`

— Method`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`

.

`QuadGK.kronrod`

— Method`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.