repository·master·Indexed 43 days ago
https://github.com/juliamath/quadgk.jlA Julia package for one-dimensional numerical integration using adaptive Gauss-Kronrod quadrature. It supports arbitrary numeric types (including BigFloat), vector-valued integrands, and improper integrals with infinite limits. Key functions include `quadgk` for integration, `gauss` for Gaussian quadrature points and weights, and `kronrod` for Kronrod rules. Features include h-adaptive quadrature, support for weighted integrals, and performance optimizations via `BatchIntegrand` and internal buffer reuse with `segbuf`.
QuadGK.jl provides support for one-dimensional numerical integration using adaptive Gauss-Kronrod quadrature. Key features include:
BigFloat).quadgk for integration, gauss for Gaussian quadrature points/weights, and kronrod for Kronrod points/weights.If your integrand is vectorized (can evaluate multiple points simultaneously), you can significantly improve performance by using the BatchIntegrand wrapper. This wrapper takes your integrand along with pre-allocated input and output buffers to avoid repeated allocations during the integration process.
# Concept usage for vectorized integrands
# Requires pre-allocated buffers for input and output
batch_f = BatchIntegrand(your_vectorized_function, input_buffer, output_buffer)
integral, error = quadgk(batch_f, a, b)The quadgk function uses Gauss–Kronrod quadrature to provide both an integral estimate and an error estimate.
quadgk, the order parameter refers to the size $n'$ of the embedded Gauss rule. The total number of points used is $n = 2n' + 1$.quadgk performs $h$-adaptive quadrature. It splits the integration interval into two halves and recursively applies the rule to each half. This is particularly robust for functions with localized sharp peaks or discontinuities.quadgk defaults to order=7 (which uses $n=15$ points).# Example of changing the order in quadgk
quadgk(f, a, b, order=10)When the integration interval is symmetric (e.g., $[-b, b]$) and the weight function is symmetric ($W(x) = W(-x)$), the Jacobi matrix is "hollow" (diagonal entries $\alpha_k = 0$).
Using QuadGK.HollowSymTridiagonal is more efficient. Its constructor only requires the off-diagonal elements $\sqrt{\beta_k}$.
Note the difference in return behavior:
gauss returns all points.kronrod returns only the $x_i \le 0$ points.using QuadGK
# Constructor only needs the off-diagonal elements
J_hollow = QuadGK.HollowSymTridiagonal([sqrt_beta_1, sqrt_beta_2, ...])
# Returns all points
x, w = gauss(J_hollow, unitintegral)
# Returns only x_i <= 0 points
x, w, wg = kronrod(J_hollow, n, unitintegral)Weighted Gaussian quadrature is used to approximate integrals of the form $\int_a^b W(x) f(x) dx$ by a weighted sum $\sum_{i=1}^n w_i f(x_i)$.
This technique is particularly useful for handling 'poorly behaved' integrands (such as those that are singular, discontinuous, or highly oscillatory) by factoring the problematic behavior into a known weight function $W(x)$. As long as $W(x) \ge 0$ on the interval $(a,b)$ and the remaining factor $f(x)$ is smooth, the quadrature rule will provide fast convergence. Once a rule is constructed for a specific $W(x)$, it can be reused for many different $f(x)$ functions.
You can compute quadrature rules for weighted integrals of the form $\int_a^b W(x) f(x) dx$.
By incorporating the weight function $W(x)$ into the quadrature weights $w_i$ and points $x_i$, you can handle poorly behaved integrands (singular, discontinuous, or highly oscillatory) more efficiently. As long as the remaining factor $f(x)$ is smooth, the quadrature rule will converge quickly.
QuadGK supports arbitrary weight functions $W(x)$ and arbitrary precision.
The quadgk integrand f(x) can return complex numbers, vectors, matrices, or any type implementing ±, scalar multiplication, and norm.
complex(x) for sqrt).norm of the error. By default, this is the Euclidean (L2) norm, but you can specify a custom norm using the norm keyword argument.StaticArrays.jl (e.g., SVector) to significantly improve efficiency. quadgk is compatible with StaticArrays via multiple dispatch.# Complex integration
quadgk(x -> 1/sqrt(complex(x)), -1, 0, 1)
# Vector-valued integration with custom norm
quadgk(x -> [1,x,x^2,x^3], 0, 1, norm=v->maximum(abs, v))
# Efficient integration with StaticArrays.jl
using StaticArrays
quadgk(x -> @SVector[1,x,x^2,x^3], 0, 1)QuadGK allows you to construct quadrature points $x_i$ and weights $w_i$ for nearly any desired weight function $W(x) \ge 0$. You can specify the weight function $W(x)$ and the interval $(a,b)$ using one of two methods:
FastGaussQuadrature.jl).quadgk. Note that this method can be computationally expensive to achieve high accuracy, but it is beneficial if the resulting rule is reused many times or if $f(x)$ is very expensive to evaluate.For integrands with singularities of the form $(1-x)^\alpha (1+x)^\beta$ at the endpoints, standard adaptive quadrature (like quadgk) can be slow. You can accelerate convergence by using Gauss-Jacobi quadrature, which incorporates these singularities into the weight function.
Workflow:
FastGaussQuadrature.jacobi_jacobimatrix(n, α, β) to get the Jacobi matrix $J_n$.FastGaussQuadrature.jacobimoment(α, β) to get the unit integral $I[1]$.gauss or kronrod in QuadGK.This method allows Gaussian quadrature to converge exponentially fast even when the integrand is not smooth at the boundaries, provided the singularities are captured by $\alpha$ and $\beta$.
using FastGaussQuadrature, QuadGK
α, β, n = 0.5, -0.1, 5
J_n = FastGaussQuadrature.jacobi_jacobimatrix(n, α, β)
I_1 = FastGaussQuadrature.jacobimoment(α, β)
# Compute points and weights for the weighted integrand
x, w = gauss(J_n, I_1)
# To integrate a smooth function f(x) with these weights:
I = sum(@. f(x) * w)While quadgk can handle integrable singularities at the endpoints (because it never evaluates the integrand exactly at the endpoints), singularities or discontinuities in the interior of the integration domain can cause errors (like DomainError) or extremely slow convergence.
To resolve this, use the syntax quadgk(f, a, b, c, ...) to subdivide the integral into multiple segments: $\int_a^b f(x)dx + \int_b^c f(x)dx + \dots$. By adding the location of the singularity or discontinuity as an extra endpoint, you ensure the algorithm handles the non-smooth behavior more efficiently.
You can also reduce the number of evaluations for subdivided intervals by lowering the order parameter.
# Problem: Singularity at x=1 in the interval (0, 2)
# This might fail or be slow:
quadgk_count(x -> 1/sqrt(abs(x-1)), 0, 2)
# Solution: Subdivide at the singularity
quadgk_count(x -> 1/sqrt(abs(x-1)), 0, 1, 2)
# Problem: Discontinuity at x=1
# Solution: Subdivide and optionally lower the order
quadgk_count(x -> x > 1, 0, 1, 3, order=1)You can compute integrals over infinite or semi-infinite intervals by passing Inf or -Inf as the integration endpoints.
Important Performance Note: quadgk uses a change-of-variables transformation for infinite limits. This works best if your function decays over a lengthscale of approximately $\sim 1$. If your function decays over a much larger or smaller lengthscale $L$, the integration may be inefficient. In such cases, you should manually rescale the integral using the substitution $u = x/L$ to improve convergence.
Example of rescaling an integral $\int_a^b f(x)dx$ where $f$ decays over lengthscale $L$: $$\int_a^b f(x)dx = \int_{a/L}^{b/L} f(uL) L , du$$
# Standard infinite integral
quadgk(x -> exp(-x), 0, Inf)
# Rescaling for a large decay lengthscale L
L = 1e6
f(x) = exp(-x/L)
quadgk(u -> f(u*L) * L, 0, Inf)quadgk is an h-adaptive quadrature method that subdivides the integration domain into pieces and sums them in an arbitrary order. This makes it unsuitable for conditionally convergent integrals (where the integral of the absolute value of the integrand diverges, such as the Dirichlet integral $\int_0^\infty \frac{\sin x}{x} dx$).
Applying quadgk naively to conditionally convergent integrals often results in extremely large error estimates and incorrect results.
Recommended Workaround: Richardson Extrapolation
You can compute these integrals by taking the limit of a Laplace transform. For example, to compute $\int_0^\infty \frac{\sin x}{x} dx$, you can integrate $e^{-sx} \frac{\sin x}{x}$ (which is absolutely convergent for $s > 0$) and then use the Richardson.jl package to extrapolate the limit as $s \to 0^+$.
using Richardson
using QuadGK
# Extrapolating the limit s -> 0+
extrapolate(1, rtol=1e-8) do s
quadgk(x -> exp(-s*x) * sin(x)/x, 0, Inf, rtol=1e-10)[1]
end