GitHub

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 — Method
gausslegendre(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
source

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 — Method
gausshermite(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
source

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 — Method
gausslaguerre(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
source
FastGaussQuadrature.gausslaguerre — Method
gausslaguerre(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-Welsch
  • gausslaguerre_rec: computation based on Newton iterations applied to evaluation using the recurrence relation
  • gausslaguerre_asy: the asymptotic expansions
source

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 — Method
gausschebyshevt(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
source
FastGaussQuadrature.gausschebyshevu — Method
gausschebyshevu(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
source
FastGaussQuadrature.gausschebyshevv — Method
gausschebyshevv(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
source
FastGaussQuadrature.gausschebyshevw — Method
gausschebyshevw(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
source

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 — Method
gaussjacobi(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
source

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 — Method
gaussradau(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
source

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 — Method
gausslobatto(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)
source