`quadgk`

examples

The following are several examples illustrating the usage of the main `quadgk`

numerical-integration function of QuadGK, focusing on more complicated circumstances than the smooth scalar integral of the Quick start section.

## Improper integrals: Infinite limits

`quadgk`

supports "improper" integrals over infinite and semi-infinite intervals, simply by passing `±Inf`

for the endpoints.

For example, $\int_0^\infty e^{-x} dx = 1$ is computed by:

```
julia> quadgk(x -> exp(-x), 0, Inf)
(1.0, 4.507383379289404e-11)
```

which give gives the correct answer (1) exactly in this case. Note that the error estimate `≈ 4.5e-11`

is pessimistic, as is often the case.

The Gaussian integral $\int_{-\infty}^{+\infty} e^{-x^2} dx = \sqrt{\pi} = 1.772453850905516027298167483341\ldots$ is computed by:

```
julia> quadgk(x -> exp(-x^2), -Inf, Inf)
(1.7724538509055137, 6.4296367126505234e-9)
```

which is the correct answer to nearly machine precision, despite the pessimistic error estimate `≈ 6.4e-9`

.

Internally, `quadgk`

handles infinite limits by the changes of variables

\[\int_a^\infty f(x)dx = \int_0^1 f\left(a + \frac{t}{1-t}\right) \frac{1}{(1-t)^2} dt\]

and

\[\int_{-\infty}^\infty f(x)dx = \int_{-1}^1 f\left(\frac{t}{1-t^2}\right) \frac{1+t^2}{(1-t^2)^2} dt\]

respectively. Although the transformed integrands are singular at the endpoints $t = 1$ and $t = \pm 1$, respectively, because the singularities are integrable and `quadgk`

*never evaluates the integrand exactly at the endpoints*, it is able to perform the numerical integration successfully.

Important tip: This change-of-variables trick works best if your function decays with $|x|$ over a lengthscale of order $\sim 1$. If your decay length is much larger or shorter than that, it will perform poorly. For example, with $f(x) = e^{-x/10^6}$ the decay is over a lengthscale $\sim 10^6$ and `quadgk`

requires many more evaluations (705) than for $e^{-x}$ (135, as measured by `quadgk_count`

, described below):

```
julia> quadgk_count(x -> exp(-x), 0, Inf)
(1.0, 4.507382674563286e-11, 135)
julia> quadgk_count(x -> exp(-x/1e6), 0, Inf)
(1.000000000001407e6, 0.0014578728207218505, 705)
```

If your function decays over a lengthscale $\sim L$, it is a good idea to compute improper integrals using a change of variables $\int_a^b f(x)dx = \int_{a/L}^{b/L} f(uL) L \, du$, for example:

```
julia> L = 1e6 # decay lengthscale
1.0e6
julia> f(x) = exp(-x/L)
f (generic function with 1 method)
julia> quadgk_count(u -> f(u*L) * L, 0, Inf) # rescaled integration over u = x/L
(1.0e6, 4.507388807828416e-5, 135)
```

## Counting and printing integrand evaluations

Often, it is a good idea to count the number of times the integrand is evaluated, in order to have a sense of how efficiently `quadgk`

is performing the integral; this is especially useful with badly behaved integrand (e.g. with singularities, discontinuities, sharp spikes, and/or rapid oscillations) to see whether some transformation of the problem might be helpful (see below).

This is easy enough by simply incrementing a global counter in your integrand function and printing progress as desired. But it is such a common desire in pedagogy and debugging that we provide convenience functions `quadgk_count`

and `quadgk_print`

to automate this task.

For example, in our $\int_0^\infty e^{-x} dx$ example from above, we could do:

```
julia> quadgk_count(x -> exp(-x), 0, Inf)
(1.0, 4.507383379289404e-11, 135)
```

to return the number of function evaluations (`135`

) along with the integral (`1.0`

) and error estimate (`≈ 4.5e-11`

). A relatively large number of function evaluations are required, even though the function $e^{-x}$ is very smooth, because the infinite endpoint implicitly introduces a singularity (via the change of variables discussed above).

We can also print the evaluation points, setting a lower requested relative accuracy of `rtol=1e-2`

so that we don't get so much output, by:

```
julia> quadgk_print(x -> exp(-x), 0, Inf, rtol=1e-2)
f(1.0) = 0.36787944117144233
f(0.655923922306948) = 0.5189623601162878
f(1.52456705113438) = 0.21771529605384055
f(0.026110451522445993) = 0.9742274787577301
f(38.29883980138536) = 2.3282264081806294e-17
f(0.00429064542600238) = 0.9957185462423211
f(233.06516868994527) = 6.04064500678147e-102
f(0.14841469193194298) = 0.8620735458624529
f(6.737877409458626) = 0.0011851600983136456
f(0.07246402202084404) = 0.9300992091482783
f(13.799951646519888) = 1.015680581505947e-6
f(0.42263178703605514) = 0.6553198859704107
f(2.3661258586654506) = 0.09384358625528809
f(0.26096469051417465) = 0.7703081183184751
f(3.831936029467112) = 0.021667625834125973
(0.9999887201849575, 0.0009180738585039538, 15)
```

to see that (for a relatively low requested relative accuracy of `rtol=1e-2`

) it evaluates the integrand only 15 times at points from $x \approx 0.00429$ to $x \approx 233.1$, and still managed to get about 5 significant digits correct. (Note that `quadgk_print`

again returns a 3-element tuple, like `quadgk_count`

, where the third element is the number of integrand evaluations.)

## Integrands with singularities and discontinuities

The integral $\int_0^1 x^{-1/2} dx = \left. 2 \sqrt{x} \right|_0^1 = 2$ is perfectly finite even though the integrand $1/\sqrt{x}$ blows up at $x=0$. This is an example of an *integrable singularity*, and `quadgk`

can compute this integral:

```
julia> quadgk_count(x -> 1/sqrt(x), 0, 1)
(1.9999999845983916, 2.3762511924588765e-8, 1305)
```

Notice the large number (`1305`

) of integrand evaluations returned by `quadgk_count`

: this is an indication of how much more work it is to evaluate an integral with a singularity or any other form of non-smoothness. At its heart, the Gauss–Kronrod algorithm employed by `quadgk`

works by interpolating the integrand with polynomials over segments of the domain, and polynomials are bad at representing non-analytic functions like $1/\sqrt{x}$.

The good news is that `quadgk`

**never evaluates functions exactly at the endpoints**, so it is okay if your function blows up or errors at those points. (However, you may have to relax your error tolerance because of the slow convergence, and floating-point limitations may prevent `quadgk`

from reaching very low error tolerances for singular integrands.) Of course, it is always better to remove the singularity by some analytical transformation if you can. For example, if you need $\int_0^a f(x) x^{-1/2} dx$, you can do a change of variables $x = y^2$ to obtain an equivalent integral $\int_0^\sqrt{a} f(y^2) 2 dy$ that has no singularity and will therefore converge *much* more quickly.

If your integrand blows up (or has *any* singularity or discontinuity) in the *interior* of the integration domain, you should *add an extra "endpoint"* at that point to make sure we never evaluate it. (Also, `quadgk`

can often converge more quickly if you tell it where your singularities are via the endpoints.) For example, suppose we are integrating $\int_0^2 |x-1|^{-1/2} dx = 4$, which has an (integrable) singularity at $x=1$. If we don't tell `quadgk`

about the singularity, it gets "unlucky" and evaluates the integrand exactly at $x=1$, which ends up throwing an error:

```
julia> quadgk_count(x -> 1/sqrt(abs(x-1)), 0, 2)
ERROR: DomainError with 1.0:
integrand produced NaN in the interval (0, 2)
...
```

Instead, if we *tell* it to subdivide the integral at $x=1$, we get the correct answer(`≈ 4`

):

```
julia> quadgk_count(x -> 1/sqrt(abs(x-1)), 0, 1, 2)
(3.9999999643041515, 5.8392038954259235e-8, 2580)
```

In general, the syntax `quadgk(f, a, b, c, ...)`

denotes the integral $\int_a^b f(x)dx + \int_b^c f(x)dx + \cdots$, and `quadgk`

never evaluates the integrand $f(x)$ exactly at the endpoints $a, b, c, \ldots$.

As another example, consider an integral $\int_0^3 H(x-1) = 2$ of the discontinuous Heaviside step function $H(x)$, which $=1$ when $x > 0$ and $=0$ when $x \le 0$:

```
julia> quadgk_count(x -> x > 1, 0, 3)
(2.0000000043200235, 1.7916158219741817e-8, 705)
```

Even though $H(x-1)$ is nearly constant, `quadgk`

struggles to integrate it (`705`

function evaluations to get about 8 digits), thanks to the discontinuity at $x=1$. (Note that `true`

and `false`

in Julia are equal to numeric `0`

and `1`

, which is why we could implement $H(x-1)$ as simply `x > 1`

.) On the other hand, if we *tell it* the location of the discontinuity:

```
julia> quadgk_count(x -> x > 1, 0, 1, 3)
(2.0, 0.0, 30)
```

then it gives the *exact* answer in only 30 evaluations. The reason it takes 30 evaluations is because `quadgk`

defaults to 7th-order Gauss–Kronrod integration rule, which uses 15 points to interpolate with a high-degree polynomial. Once we subdivide the integral, we could actually get away with a lower-order rule by setting the `order`

parameter, e.g.:

```
julia> quadgk_count(x -> x > 1, 0, 1, 3, order=1)
(2.0, 0.0, 6)
```

## Complex and vector-valued integrands

The integrand `f(x)`

can return not just real numbers, but also complex numbers, vectors, matrices, or any Julia type supporting `±`

, multiplication by scalars, and `norm`

(i.e. implementing any Banach space).

For example, we can integrate $1/\sqrt{x}$ from $x=-1$ to $x=1$, where we tell the `sqrt`

function to return a complex result for negative arguments:

```
julia> quadgk(x -> 1/sqrt(complex(x)), -1, 0, 1)
(1.9999999891094182 - 1.9999999845983916im, 4.056765398346683e-8)
```

which correctly gives $\approx 2 - 2i$. Note that we explicitly put an endpoint at $x=0$ to tell `quadgk`

about the singularity at that point, as described above.

Or let's integrate the vector-valued function $f(x) = [1, x, x^2, x^3]$ for $x \in (0,1)$:

```
julia> quadgk(x -> [1,x,x^2,x^3], 0, 1)
([1.0, 0.5, 0.3333333333333333, 0.25], 6.206335383118183e-17)
```

which correctly returns $\approx [1, \frac{1}{2}, \frac{1}{3}, \frac{1}{4}]$. Note that the error estimate in this case is an approximate bound on the norm of the error, as computed by the `LinearAlgebra.norm`

function in Julia. It defaults to the Euclidean (L2) norm, but you can change this with the `norm`

argument:

```
julia> quadgk(x -> [1,x,x^2,x^3], 0, 1, norm=v->maximum(abs, v))
([1.0, 0.5, 0.3333333333333333, 0.25], 5.551115123125783e-17)
```

e.g. to use the maximum norm or some other norm (e.g. a weighted norm if the different components have different units or have unequal error tolerances).

For integrands whose values are *small* arrays whose length is known at compile time, it is usually most efficient to modify your integrand to return an `SVector`

from the StaticArrays.jl package. For the example above:

```
julia> using StaticArrays
julia> integral, error = quadgk(x -> @SVector[1,x,x^2,x^3], 0, 1)
([1.0, 0.5, 0.3333333333333333, 0.25], 6.206335383118183e-17)
julia> typeof(integral)
SVector{4, Float64} (alias for SArray{Tuple{4}, Float64, 1, 4})
```

Note that the return value also gives the `integral`

as an `SVector`

(a statically sized array).

The QuadGK package did not need any code specific to StaticArrays, and was written long before that package even existed. The fact that unrelated packages like this can be composed is part of the beauty of multiple dispatch and duck typing for generic programming.

## Batched integrand evaluation

User-side parallelization of integrand evaluations is also possible by providing an in-place function of the form `f!(y,x) = y .= f.(x)`

, which evaluates the integrand at multiple points simultaneously. To use this API, `quadgk`

dispatches on a `BatchIntegrand`

type containing `f!`

and buffers for `y`

and `x`

. These buffers may be pre-allocated and reused for multiple `BatchIntegrand`

s with the same domain and range types.

For example, we can perform multi-threaded integration of a highly oscillatory function that needs to be refined globally:

```
julia> f(x) = sin(100x)
f (generic function with 1 method)
julia> function f!(y, x)
n = Threads.nthreads()
Threads.@threads for i in 1:n
y[i:n:end] .= f.(@view(x[i:n:end]))
end
end
f! (generic function with 1 method)
julia> quadgk(BatchIntegrand{Float64}(f!), 0, 1)
(0.0013768112771231598, 8.493080824940099e-12)
```

Batching also changes how the adaptive refinement is done, which typically leads to slightly different results and sometimes more integrand evaluations. You can limit the maximum batch size by setting the `max_batch`

parameter of the `BatchIntegrand`

, which can be useful in order to set an upper bound on the size of the buffers allocated by `quadgk`

.

## Arbitrary-precision integrals

`quadgk`

also supports arbitrary-precision arithmetic using Julia's `BigFloat`

type to compute integrals to arbitrary accuracy (albeit at increased computational cost).

For example, we can compute the error function $\frac{\sqrt{\pi}}{2} \text{erf}(1) = \int_0^1 e^{-x^2} dx$ to 50 digits by:

```
julia> setprecision(60, base=10) # use 60-digit arithmetic
60
julia> quadgk_count(x -> exp(-x^2), big"0.0", big"1.0", rtol=1e-50)
(0.74682413281242702539946743613185300535449968681260632902766195, 6.8956257635323481758755998484087241330474674891762053644928492e-51, 15345)
```

The correct answer is `≈ 0.746824132812427025399467436131853005354499686812606329027654498958605…`

, and we are matching that to nearly the full precision (≈ 60 digits). (As usual, the error estimate of `quadgk`

is very conservative for smooth functions.)

Unfortunately, it took 15345 function evaluations to obtain such an accurate answer. Since this a smooth integrand, for high-accuracy calculations it is often advisable to increase the "order" of the quadrature algorithm (which is related to the degree of polynomials used for interpolation). The default is `order=7`

, but let's try tripling it to `order=21`

:

```
julia> quadgk_count(x -> exp(-x^2), big"0.0", big"1.0", rtol=1e-50, order=21)
(0.74682413281242702539946743613185300535449968681260632902765324, 2.1873898701681913100611385149037136705674736373054902472850425e-58, 129)
```

It got the same accuracy (≈ 60 digits) with only 129 integrand evaluations!

## Contour integration

You can specify a sequence of points in the complex plane to perform a contour integrals with `quadgk`

along a piecewise-linear contour.

For example, consider the function $f(z) = \cos(z)/z$. By the residue theorem of complex analysis, if we integrate counter-clockwise in a "loop" around the pole at $z=0$, we should get exactly $2\pi i \cos(0) = 2\pi i$.

One way to do this integral is to parameterize a contour, say a circle $|z|=1$ parameterized by $z= e^{i\phi}$ (`= cis(ϕ)`

in Julia), which gives $dz = i z\, d\phi$, to obtain an ordinary integral $\int_0^{2\pi} \frac{\cos(e^{i\phi})}{e^{i\phi}} ie^{i\phi} d\phi$ over the *real* parameter $\phi \in (0,2\pi)$:

```
julia> quadgk(ϕ -> cos(cis(ϕ)) * im, 0, 2π)
(0.0 + 6.283185307179586im, 1.8649646913725044e-8)
```

which indeed gives us $\approx 2\pi i$ (to machine precision).

As an alternative, however, you can directly supply a sequence of *complex "endpoints"* to `quadgk`

and it will perform the contour integral along a sequence of line segments connecting these points. For example, instead of integrating around a circular contour, we can integrate around the diamond (rotated square) connecting the corners $\pm 1$ and $\pm i$:

```
julia> quadgk(z -> cos(z)/z, 1, im, -1, -im, 1)
(0.0 + 6.283185307179587im, 5.369976662961913e-9)
```

which again gives $\approx 2\pi i$ (to machine precision). Note that it is critically important to have a *closed* contour (a loop): the final endpoint must be the same as the starting point ($z=1$).

## Cauchy principal values

Integrands $f(x) = g(x)/x$ that diverge $\sim 1/x$ cannot be integrated through $x=0$ in the usual way (the singularity is not integrable). However, if you integrate *around* $x=0$, for both signs of $x$, then you can define a kind of integral that is the "difference" of the divergence on the two sides. This definition is called a Cauchy principal value, and is usually presented as a limit:

\[\text{p.v.}\int_a^b \frac{g(x)}{x} dx = \lim_{\varepsilon\to 0^+} \left[ \int_a^{-\varepsilon} \frac{g(x)}{x} dx + \int_{+\varepsilon}^b \frac{g(x)}{x} dx \right] \, .\]

That is, you subtract a "ball" of radius $\varepsilon$ from the integration domain $a < 0 < b$ to eliminate the singularity at $x=0$, and take the limit of the resulting integral as the $\varepsilon$ goes to zero.

In principle, you might imagine taking this limit numerically by extrapolation of numerical integrals for a sequence of $\varepsilon > 0$ values, perhaps using Richardson extrapolation via the Richardson.jl package. However, it is mathematically equivalent and *much* more efficient to use a simple singularity-subtraction procedure:

\[\frac{g(x)}{x} = \frac{g(x)-g(0)}{x} + \frac{g(0)}{x}\]

where the first term is *not singular* if $g(x)$ is differentiable at $x=0$, and the latter term can be integrated analytically, giving:

\[\text{p.v.}\int_a^b \frac{g(x)}{x} dx = \int_a^b \frac{g(x)-g(0)}{x} dx + g(0) \log|b/a|\]

Since the remaining integral has no singularity, we can do it numerically directly. There are two tricks to help us a bit further:

- As in Integrands with singularities and discontinuities above, we'll want to put an extra endpoint at $x=0$ to make sure
`quadgk`

doesn't evaluate the integrand exactly at that point (which would give`NaN`

from $0/0$). - We should be careful with the integration tolerances, to make sure that any relative tolerance
`rtol`

is applied with respect to the whole principal part and not just to the $g(x)-g(0)$ integral. An easy way to do this is to add $g(0) \frac{\log|b/a|}{b-a}$ to the*integrand*, so that we no longer compute the two pieces separately.

Putting it all together, here is a function `cauchy_quadgk(g, a, b)`

that computes our Cauchy principal part $\text{p.v.}\int_a^b g(x)/x$:

```
function cauchy_quadgk(g, a, b; kws...)
a < 0 < b || throw(ArgumentError("domain must include 0"))
g₀ = g(0)
g₀int = b == -a ? zero(g₀) : g₀ * log(abs(b/a)) / (b - a)
return quadgk_count(x -> (g(x)-g₀)/x + g₀int, a, 0, b; kws...)
end
```

For example, Mathematica tells us that

\[\text{p.v.} \int_{-1}^2 \frac{\cos(x^2-1)}{x} dx \approx 0.212451309942989788929352736695\ldots ,\]

and we can reproduce this with `cauchy_quadgk`

:

```
julia> cauchy_quadgk(x -> cos(x^2-1), -1, 2)
(0.21245130994298977, 1.8366794196644776e-11, 60)
```

which is correct to about 16 digits.

This approach and other approaches to computing Cauchy principal values are discussed in Keller and Wróbel (2016). This kind of "singularity subtraction" is a powerful approach to efficient computation of integrals with singularities or near singularities. A huge variety of related techniques have been developed for boundary element methods, where a vast number of singular integrals must be computed and efficiency is at a premium. See, for example, Reid *et al.* (2014) and references therein.

## Nearly singular integrands

Even if the integrand is only *nearly* singular, so that there is a sharp but *finite* peak within the integration domain, it can greatly increase the efficiency of numerical integration if you can separate the sharp peak analytically.

For example, suppose that you are integrating:

\[I = \int_a^b \frac{g(x)}{x - i\alpha} dx\]

for a small $0 < \alpha \ll 1$. For $\alpha \to 0^+$, it approaches $i\pi g(0)$ plus a Cauchy principal part (the latter being zero if $a = -b$ and $g(x)=g(-x)$), but for small $\alpha > 0$ you have to numerically integrate (for a general function $g(x)$) a function with a sharp spike at $x=0$, which will require a large number of quadrature points. But you can subtract out the singularity analytically:

\[I = \int_a^b \left[ \frac{g(x)-g(0)}{x - i\alpha} + \frac{g(0)}{x - i\alpha} \right] dx \\ = \int_a^b \frac{g(x)-g(0)}{x - i\alpha}dx + \underbrace{g(0) \left[\frac{1}{2}\log(x^2 + \alpha^2) + i\tan^{-1}(x/\alpha) \right]_a^b}_{I_0}\]

and then you only need to numerically integrate $I - I_0$, which has the spike subtracted.

As for Cauchy principal values above, we want to include a $I_0 / (b-a)$ term directly in the integrand so that the error tolerances are computed correctly, and include $x=0$ as an explicit endpoint to let `quadgk`

know that the integral is badly behaved there.

In code:

```
using QuadGK
function int_slow(g, α, a, b; kws...)
if a < 0 < b
# put an explicit endpoint at x=0 since we know it is badly behaved there
return quadgk_count(x -> g(x) / (x - im*α), a, 0, b; kws...)
else
return quadgk_count(x -> g(x) / (x - im*α), a, b; kws...)
end
end
function int_fast(g, α, a, b; kws...)
g₀ = g(0)
denom_int(x) = log(x^2 + α^2)/2 + im * atan(x/α)
I₀ = g₀ * (denom_int(b) - denom_int(a))
if a < 0 < b
# put an explicit endpoint at x=0 since we know it is badly behaved there
(I,E,c) = quadgk_count(x -> I₀/(b-a) + (g(x) - g₀) / (x - im*α), a, 0, b; kws...)
else
(I,E,c) = quadgk_count(x -> I₀/(b-a) + (g(x) - g₀) / (x - im*α), a, b; kws...)
end
return (I,E,c+1) # add 1 for g(0) evaluation
end
```

This gives:

```
julia> int_slow(cos, 1e-6, -1, 1)
(1.1102230246251565e-16 + 3.1415896808206125im, 1.4895091715264936e-9, 1230)
julia> int_fast(cos, 1e-6, -1, 1)
(3.3306690738754696e-16 + 3.1415896808190418im, 1.8459683038047577e-12, 31)
```

which agree to about 13 digits, but the slow brute-force method requires 1230 function evaluations while the fast singularity-subtracted method requires only 31 function evaluations.

As an added bonus, `int_fast`

works even for `α = 0`

, where it gives you $i\pi g(0)$ (for $0 \in (a,b)$) plus the Cauchy principal part as above.