Gaussian Quadrature
Gauss-Legendre quadrature
Gauss quadrature for the weight function $w(x) = 1$.
- For $n ≤ 5$: Use an analytic expression.
- For $n ≤ 60$: Use Newton's method to solve $P_n(x)=0$. Evaluate Legendre polynomials $P_n$ and their derivatives $P_n'$ by 3-term recurrence. Weights are related to $P_n'$.
- For $n > 60$: Use asymptotic expansions for the Legendre nodes and weights [1].
FastGaussQuadrature.gausslegendre
— Methodgausslegendre(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Legendre quadrature.
\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausslegendre(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 2/5
true
Gauss-Hermite quadrature
Gauss quadrature for the weight function $w(x) = \exp(-x^2)$ on the real line.
- For $n < 200$: Use Newton's method to solve $H_n(x)=0$. Evaluate Hermite polynomials $H_n$ and their derivatives $H_n'$ by three-term recurrence.
- For $n ≥ 200$: Use Newton's method to solve $H_n(x)=0$. Evaluate $H_n$ and $H_n'$ by a uniform asymptotic expansion, see [7].
The paper [7] also derives an $O(n)$ algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form $\exp(-V(x))$, where $V(x)$ is a real polynomial.
FastGaussQuadrature.gausshermite
— Methodgausshermite(n::Integer; normalize = false) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Hermite quadrature.
\[\int_{-\infty}^{+\infty} f(x) \exp(-x^2) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
The option normalize=true
instead computes
\[\int_{-\infty}^{+\infty} f(x) \phi(x) dx \approx \sum_{i=1}^{n} w_i f(x_i),\]
where $\phi$ is the standard normal density function.
Examples
julia> x, w = gausshermite(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 3(√π)/4
true
Gauss-Laguerre quadrature
Gauss quadrature for the weight function $w(x) = \exp(-x)$ on $[0,+\infty)$
- For $n < 128$: Use the Golub-Welsch algorithm.
- For
method=GLR
: Use the Glaser-Lui-Rohklin algorithm. Evaluate Laguerre polynomials $L_n$ and their derivatives $L_n'$ by using Taylor series expansions near roots generated by solving the second-order differential equation that $L_n$ satisfies, see [2]. - For $n ≥ 128$: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight $w(x) = x^\alpha \exp(-q_m x^m)$ and this is $O(\sqrt{n})$ when allowed to stop when the weights are below the smallest positive floating point number.
FastGaussQuadrature.gausslaguerre
— Methodgausslaguerre(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Laguerre quadrature.
\[\int_{0}^{+\infty} f(x) e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausslaguerre(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 24
true
FastGaussQuadrature.gausslaguerre
— Methodgausslaguerre(n::Integer, α::Real) -> x, w # nodes, weights
Return nodes x
and weights w
of generalized Gauss-Laguerre quadrature.
\[\int_{0}^{+\infty} f(x) x^\alpha e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausslaguerre(3, 1.0);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 120
true
Optionally, a reduced quadrature rule can be computed. In that case, only those points and weights are computed for which the weight does not underflow in the floating point precision type. Supply the optional argument reduced = true
.
Though the code is generic, heuristical choices on the choice of the algorithm are based on achieving machine precision accuracy only for Float64
type. In case the default choice of algorithm does not achieve the desired accuracy, the user can manually invoke the following routines:
gausslaguerre_GW
: computation based on Golub-Welschgausslaguerre_rec
: computation based on Newton iterations applied to evaluation using the recurrence relationgausslaguerre_asy
: the asymptotic expansions
Gauss-Chebyshev quadrature
There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions:
- 1st kind, weight function $w(x) = 1/\sqrt{1-x^2}$
- 2nd kind, weight function $w(x) = \sqrt{1-x^2}$
- 3rd kind, weight function $w(x) = \sqrt{(1+x)/(1-x)}$
- 4th kind, weight function $w(x) = \sqrt{(1-x)/(1+x)}$
They are all have explicit simple formulas for the nodes and weights [4].
FastGaussQuadrature.gausschebyshevt
— Methodgausschebyshevt(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Chebyshev quadrature of the 1st kind.
\[\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausschebyshevt(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 3π/8
true
FastGaussQuadrature.gausschebyshevu
— Methodgausschebyshevu(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Chebyshev quadrature of the 2nd kind.
\[\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausschebyshevu(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ π/16
true
FastGaussQuadrature.gausschebyshevv
— Methodgausschebyshevv(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Chebyshev quadrature of the 3rd kind.
\[\int_{-1}^{1} f(x)\sqrt{\frac{1+x}{1-x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausschebyshevv(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 3π/8
true
FastGaussQuadrature.gausschebyshevw
— Methodgausschebyshevw(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Chebyshev quadrature of the 4th kind.
\[\int_{-1}^{1} f(x)\sqrt{\frac{1-x}{1+x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausschebyshevw(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 3π/8
true
Gauss-Jacobi quadrature
Gauss quadrature for the weight functions $w(x) = (1-x)^\alpha(1+x)^\beta$, $\alpha,\beta > -1$.
- For $n ≤ 100$: Use Newton's method to solve $P_n(x)=0$. Evaluate $P_n$ and $P_n'$ by three-term recurrence.
- For $n > 100$: Use Newton's method to solve $P_n(x)=0$. Evaluate $P_n$ and $P_n'$ by an asymptotic expansion (in the interior of $[-1,1]$) and the three-term recurrence $O(n^{-2})$ close to the endpoints. (This is a small modification to the algorithm described in [3].)
- For $\max(\alpha,\beta) > 5$: Use the Golub-Welsch algorithm requiring $O(n^2)$ operations.
FastGaussQuadrature.gaussjacobi
— Methodgaussjacobi(n::Integer, α::Real, β::Real) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Jacobi quadrature for exponents α
and β
.
\[\int_{-1}^{1} f(x) (1-x)^\alpha(1+x)^\beta dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gaussjacobi(3, 1/3, -1/3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 268π/729(√3)
true
Gauss-Radau quadrature
Gauss quadrature for the weight function $w(x)=1$, except the endpoint $-1$ is included as a quadrature node.
The Gauss-Radau nodes and weights can be computed via the $(0,1)$ Gauss-Jacobi nodes and weights[3].
FastGaussQuadrature.gaussradau
— Methodgaussradau(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Radau quadrature.
\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gaussradau(3);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 2/5
true
Note that the first node is fixed at -1.
julia> x, w = gaussradau(3);
julia> x[1]
-1.0
Gauss-Lobatto quadrature
Gauss quadrature for the weight function $w(x)=1$, except the endpoints $-1$ and $1$ are included as nodes.
The Gauss-Lobatto nodes and weights can be computed via the $(1,1)$ Gauss-Jacobi nodes and weights[3].
FastGaussQuadrature.gausslobatto
— Methodgausslobatto(n::Integer) -> x, w # nodes, weights
Return nodes x
and weights w
of Gauss-Lobatto quadrature.
\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]
Examples
julia> x, w = gausslobatto(4);
julia> f(x) = x^4;
julia> I = dot(w, f.(x));
julia> I ≈ 2/5
true
Note that the both ends of nodes are fixed at -1 and 1.
julia> x, w = gausslobatto(4);
julia> x[1], x[end]
(-1.0, 1.0)