The following is from my notes I did a while back, whilst evaluating different numerical quadrature methods. Quite a lot of numerical quadrature in software used seems to be Gauss-Kronrod, in particular derived from QUADPACK.
Any thoughts or corrections? Drop me a mail.
Gaussian Quadrature
Gaussian quadrature approximates integrals in the form of: $$ \int^{b}_{a}f(x)\ dx $$ For some range \(a\) to \(b\) inclusive or exclusive depending upon the type. With a sum in the form of: $$ \int^b_af(x)\ dx\approx\sum_{i=1}^nw_if(x_i) $$ Where \(w_i\) is some weight and the node \(x_i\) being the root of an orthogonal polynomial.
Gauss-Chebyshev Type 1
Gauss-Chebyshev type 1 quadrature numerically approximates integrals in the form of: $$ \int^{+1}_{-1}\frac{f(x)}{\sqrt{1-x^2}}\ dx \approx \sum^n_{i=1}w_if(x_i) $$ Where the abscissas are: $$ x_i=\cos\left(\frac{2i-1}{2n}\pi\right) $$ and the weight \(w_i=\pi/n\). One thing to mention is that if $$ n\not\equiv 0\mod 2 $$ the sample points will be biased towards one side, which may be undesirable if high precision is required.
This is a special case of Gauss-Jacobi quadrature where \(\alpha=\beta=-\frac{1}{2}\) and with a static weight as above.
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Gauss-Chebyshev Type 2
The type 2 Gauss-Cherbychev quadrature approximates the integral $$ \int^{+1}_{-1}\sqrt{1-x^2}g(x)\ dx\approx \sum^n_{i=1}w_if(x_i) $$ with the abscissas being defined as: $$ x_i=\cos\left(\frac{i}{n+1}\pi\right) $$ and with the weights of: $$ w_i=\frac{\pi}{n+1}\sin^2\left(\frac{i}{n+1}\pi\right) $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Gauss-Hermite
Gauss-Hermite qudrature approximates integrals in the form of: $$ \int_{-\infty}^{\infty}e^{-x^2}f(x)\ dx\approx \sum^n_{i=1}w_if(x_i) $$ Where the abscissas, \(x_i\) are roots of the Hermite Polynomial \(H_n\) and the weights \(w_i\) are defined as: $$ w_i=\frac{2^{n-1}n!\sqrt{\pi}}{n^2\left[H_{n-1}(x_i)\right]^2} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Gauss-Jacobi
Gauss-Jacobi quadrature approximates integrals of the form: $$ \int_{-1}^1f(x)(1-x)^\alpha(1+x)^\beta\ dx \approx \sum^n_{i=1}w_if(x_i) $$ Where the abscissas \(x_i\) are the roots of the Jacobi Polynomial \(P_n^{(\alpha,\beta)}\) and the weights are: $$ w_i=-\frac{2n+\alpha+\beta+2}{n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+\beta+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)(n+1)!}\frac{2^{\alpha+\beta}}{\left(\frac{d}{dx}P_n^{(\alpha,\beta)}(x_i)\right)P_{n+1}^{(\alpha,\beta)}(x_i)} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\), \(\alpha=10\) and \(\beta=50\)
(more for show, than for applicability with values as large as this):
Gauss-Gegenbauer
This is a special case of Gauss-Jacobi quadrature where \(\alpha=\beta\). In this case the weight simplifies to: $$ w_i=\frac{2(n+\alpha+1)}{n+2\alpha+1}\frac{\Gamma(n+\alpha+1)^2}{\Gamma(n+2a+1)(n+1)!}\frac{2^{2\alpha}}{\left(\frac{d}{dx}P_n^{(\alpha,\beta)}(x_i)\right)P_{n+1}^{(\alpha,\beta)}(x_i)} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\) and \(\alpha=10\):
Gauss-Laguerre
Gauss-Laguerre quadrature approximates integrals in the form of: $$ \int_0^{\infty}e^{-x}f(x)\ dx\approx \sum^n_{i=1}w_if(x_i) $$ The abscissas \(x_i\) are the i-th root of the Laguarre Polynomial \(L_n\) and the weights are computed using: $$ w_i=\frac{x_i}{(n+1)^2\left[L_{n+1}(x_i)\right]^2} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Gauss-Laguerre (Generalized)
Generalized Gauss-Laguerre quadrature approximates integrals in the form of: $$ \int_0^{\infty}x^\alpha e^{-x}f(x)\ dx\approx \sum^n_{i=1}w_if(x_i) $$ The abscissas \(x_i\) are the i-th root of the Generalized Laguerre Polynomial \(L^{(\alpha)}_n\) and the weights \(w_i\) are computed using: $$ w_i=\frac{\Gamma(n+\alpha+1)x_i}{n!(n+1)^2\left[L^{(\alpha)}_{n+1}(x_i)\right]^2} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\) and \(\alpha=1.2\):
Gauss-Legendre
Gauss-Legendre quadrature approximates an integral in the form of: $$ \int_{-1}^1f(x)\ dx\approx \sum^n_{i=1}w_if(x_i) $$ With the abscissas \(x_i\) being the roots of the Legendre Polynomial and the weights being: $$ w_i=\frac{2}{(1-x_i^2)\left[\frac{d}{dx}P_n(x_i)\right]^2} $$
This is a special case of Gauss-Jacobi quadrature where \(\alpha=\beta=0\).
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Extensions
Gauss-Lobatto
Gauss-Lobatto is an extension of the Gauss-Legendre quadrature, with the inclusion of the integration limits \(\pm1\) to the abscissas and the addition of a remainer \(R_n\) as defined below. These normally result in a division by zero when evaluated on the first derivative of the \(P_n\) polynomial. Gauss-Lobatto solves integrals in the form of: $$ \int_{-1}^1f(x)\ dx \approx \frac{2}{n(n-1)}\left[f(1)+f(-1)\right]+\sum^{n-1}{i=2}w_if(x_i)+R_n $$ With the weights: $$ w_i=\frac{2}{n(n-1)\left[P{n-1}(x_i)\right]^2} $$ and the remainder: $$ R_n=\frac{-n(n-1)^32^{2n-1}\left[(n-2)!\right]^4}{(2n-1)\left[(2n-2)!\right]^3}f^{(2n-2)}(\xi) $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\):
Gauss-Radau
Gauss-Radau is an extension of Gauss-Legendre, where one of the end points of the integration interval \(x_1\) is included in the evaluation. This is usually \(x_1=-1\). $$ \int_{-1}^1f(x)\ dx \approx w_1f(x_1)+\sum_{i=2}^nw_if(x_i) $$ With the abscissas \(x_i\) being the roots of: $$ \frac{P_{n-1}(x)+P_n(x)}{1+x} $$ and the weights being: $$ w_i = \begin{cases} \frac{2}{n^2} & \text{if } i=1 \\ \frac{1-x_i}{n^2\left[P_{n-1}(x_i)\right]^2} &\text{otherwise} \end{cases} $$
The following shows a plot of \(w_i\) against \(x_i\) for \(n=26\) with the choosen point being \(-1\):
Gauss-Kronrod
Gauss-Kronrod probably warrants it’s own article, as they are not as straight forward to calculate. So I won’t cover it here. For now here some references:
- Laurie, Dirk “Calculation of Gauss-Kronrod quadrature rules” (1997)
- Piessens et. al “QUADPACK: A subroutine package for automatic integration.” ISBN 978-3-540-12553-2 (1983)
MAXIMA Code
I wrote the following code for MAXIMA, which evaluates the \(x_i\) and \(w_i\) coefficents for all of the above described quadrature rules.
It is important that orthopoly_returns_intervals : false is set, as otherwise realroots() for some
of the polynomials will return a bound interval due to inaccuracies in constants and arithmetic within MAXIMA,
rather than a definite value. For the
illustrative purposes of this, an approximated value within that bound interval is fine however.
Function Reference
These functions return a list of lists, where the latter contains \(x_i\) and \(w_i\) for \(n\).
gausschebyshev1(n)- Gauss-Chebyshev Type 1gausschebyshev2(n)- Gauss-Chebyshev Type 2gausshermite(n)- Gauss-Hermitegausslaguerre(n)- Gauss-Laguerregaussgenlaguerre(n,a)- Gauss-Laguerre (Generalized)gausslegendre(n)- Gauss-Legendregaussjacobi(n,a,b)- Gauss-Jacobigaussgegenbauer(n,a)- Gauss-Gegenbauergausslobatto(n)- Gauss-Lobattogaussradau(n,p)- Gauss-Radau
These functions return a list of the nodes of the respective method.
gausschebyshev1_xi(n)- Gauss-Chebyshev Type 1gausschebyshev2_xi(n)- Gauss-Chebyshev Type 2gausshermite_xi(n)- Gauss-Hermitegausslaguerre_xi(n)- Gauss-Laguerregaussgenlaguerre_xi(n,a)- Gauss-Laguerre (Generalized)gausslegendre_xi(n)- Gauss-Legendregaussjacobi_xi(n,a,b)- Gauss-Jacobigaussgegenbauer_xi(n,a)- Gauss-Gegenbauergausslobatto_xi(n)- Gauss-Lobattogaussradau_xi(n)- Gauss-Radau (Note: Does not include the choosen endpoint)
These functions return the weight \(w_i\) associated with a given node \(x_i\):
gausschebyshev1_wi(n,x)- Gauss-Chebyshev Type 1gausschebyshev2_wi(n,x)- Gauss-Chebyshev Type 2gausshermite_wi(n,x)- Gauss-Hermitegausslaguerre_wi(n,x)- Gauss-Laguerregaussgenlaguerre_wi(n,x,a)- Gauss-Laguerre (Generalized)gausslegendre_wi(n,x)- Gauss-Legendregaussjacobi_wi(n,x,a,b)- Gauss-Jacobigaussgegenbauer_wi(n,x,a)- Gauss-Gegenbauergausslobatto_wi(n,x)- Gauss-Lobattogaussradau_wi(n,x)- Gauss-Radau (Note: Does handle the weight of the choosen endpoint)
To obtain diagrams like the above, use:
plot2d([discrete,gausschebyshev1(26)],
[title, simplode(["Gauss-Chebyshev Type 1 n=", 26])],
[xlabel, "x_i"],[ylabel, "w_i"],
[style, points],[point_type, asterisk],[yx_ratio, 0.25]);
References
Gaussian Quadrature
- Kahaner, David et al. “Numerical Methods and Software” ISBN 013-627258-4 (1989)
- Kincaid, David, Cheney, Ward “Numerical Analysis” 2nd Ed. ISBN 0534-33892-5 (1996)
- Piessens et. al “QUADPACK: A subroutine package for automatic integration.” ISBN 978-3-540-12553-2 (1983)
- Gautschi, Walter, Li, Shikang “Gauss-Radau and Gauss-Lobatto quadratures with double end points” (1990)
- Golub, Gene H., Welsch, John H. “Calculation of Gauss Quadrature Rules” (1967)
Roots of Polynomials
Non-exhaustive list, just the ones I wrote down in my notes:
- Saizer, Herbert E. et al. “Table of the Zeros and Weight Factors of the First Twenty Hermite Polynomials.” (1952).
- Masoero, Davide, Roffelsen Pieter “Roots of generalised Hermite polynomials when both parameters are large” (2019).
- Townsend, Alex et al. “Fast computation of Gauss quadrature node and weights on the whole real line” (2014)
- Salzer, Herbert E. et al. “Table of zeros and weight factors of the first fifteen Laguerre polynomials” (1949)
- Shao, T. S., Chen, T. C “Table of zeros and Gaussian Weights of certain Associated Laguerre Polynomials and the related Hermite Polynomials” (1964)
- Lowan, Arnold N. et al. “Table of the zeros of the Legendre polynomials of order 1-16 and the weight coefficients for Gauss’ mechanical quadrature formula” (1942)
- Moak, D. S. et al. “On the Zeros of Jacobi Polynomials \(P^{(|alpha_n,|beta_n)}_n(x)\)” (1979)