Abscissas and Weights of Classical Gaussian Quadrature Rules

by Pavel Holoborodko on May 30, 2012

The latest version of the Multiprecision Computing Toolbox introduces new routines to compute the nodes and weights of classical Gaussian quadrature rules with any required accuracy. A summary of the supported functionality is outlined below.

The Gaussian quadrature is targeted to approximate an integral by taking the weighted sum of integrand values sampled at special points (called abscissas). The abscissas and weights are calculated in a special way so that the rule provides a precise answer for all polynomials up to certain degree.

    \[ \int_a^b \omega(x)f(x)\,\mathrm{d}x \approx \sum_{i=1}^{n}{w_i\,f(x_i)} \]

The most common case, the Gauss-Legendre quadrature, occurs when the weight function \omega(x) = 1. However if the integrand has a factor (weight) of a special form, then a more efficient quadrature can be applied. This is beneficial for both accuracy and computational speed.

The table below lists the classical Gaussian quadrature rules with corresponding weight functions and parameters. Note that the formulas are provided in generalized form, already adjusted for the arbitrary intervals of integration.

NameWeight function, \omega(x)IntervalFunction
Chebyshev Type 1((b-x)(x-a))^{-1/2}(a,b)mp.GaussChebyshevType1
Chebyshev Type 2((b-x)(x-a))^{1/2}(a,b)mp.GaussChebyshevType2
Generalized Laguerre(x-a)^\alpha\exp(-b(x-a))(a,+\infty)mp.GaussLaguerre
Generalized Hermite|x-a|^\alpha\,\exp(-b(x-a)^2)(-\infty,+\infty)mp.GaussHermit

The last column in the table shows the functions provided by the Multiprecision Computing Toolbox for computing abscissas and weights for each corresponding quadrature with arbitrary precision. A detailed description of the functions and parameters is available using the standard MATLAB’s command help:

>> help mp.GaussLaguerre
 mp.GaussLaguerre - Arbitrary precision coefficients of Generalized Gauss-Laguerre Quadrature
     [x,w] = mp.GaussLaguerre(order, alpha, a, b) computes abscissas
     'x' and weights 'w' of Generalized Gauss-Laguerre quadrature.
     The generalized Gauss-Laguerre quadrature rule is used as follows:
         Integral ( a < x < +inf ) (x-a)^alpha * exp(-b*(x - a)) f(x) dx
     is to be approximated using 'x' and 'w' by
          sum ( 1 <= i <= order ) w(i) * f(x(i))      
     'order' is the number of points in the quadrature rule.
     'alpha' is the exponent of |x| in the weight function. 
            The value of alpha may be any real value greater than -1.0. 
            Specifying alpha = 0.0 (default) results in the basic (non-generalized) rule.
     'a' is the left point (default 0);
     'b' is the scale factor in the exponential (default 1);        
     Precision is controlled by mp.Digits and by precision of the arguments

For example, let’s compute the 10th order Gauss-Laguerre quadrature rule to an accuracy of 50 decimal digits with the default parameters (\alpha=0, a=0, b=1):

>> mp.Digits(50);
>> [x,w] = mp.GaussLaguerre(10);
>> x
>> w

Example 1. The Runge function integral approximation on [-1,1] using the Gauss-Legendre quadrature:

    \[ \int_{-1}^{1} \frac{1}{1+25x^2}\,\mathrm{d}x \]

>> mp.Digits(50);
>> f = @(x) 1./(1+25.*x.^2);
>> [x,w] = mp.GaussLegendre(50,-1,1);        % GaussLegendre(order, a, b)
>> sum(f(x).*w)

Example 2. The Gauss-Laguerre rule can be used to compute values of Gamma function:

    \[  \Gamma(z) = \int_0^\infty  e^{-t} t^{z-1}\,\mathrm{d}t \]

The function under the integral corresponds to the Laguerre weight with \alpha=z-1, a=0, b=1, and the approximation can be calculated as (z=\pi):

>> mp.Digits(50);
>> z = mp('pi');
>> [t,w] = mp.GaussLaguerre(30, z-1);
>> sum(w)
>> gamma(z)


The toolbox uses the Golub-Welsch method[1],[2] to compute all Gaussian quadrature abscissas and weights.

At first we construct a tri-diagonal Jacobi matrix using coefficients of the recurrence relation of orthogonal polynomials built with the \omega(x) weight function on the required interval. The quadrature’s nodes are eigenvalues of the Jacobi matrix, and weights are computed from the first components of the corresponding eigenvector. Efficient, implicit QL decomposition is used to find the eigenvalues and eigenvectors[3].

The code of the procedure:

    % Build tri-diagonal Jacobi matrix 
    [aj, bj, mu] = jacobi_matrix(kind, N, alpha, beta);
    J = diag(aj,0) + diag(bj(1:N-1),1) + diag(bj(1:N-1),-1);
    % Find eigenvectors & eigenvalues of tridiagonal Jacobi matrix
    % This gives us quadrature rules for default intervals
    [V,D] = eig(J);
    [x,i] = sort(diag(D));
    w = mu * V(1,i)'.^2;

Toolbox’s implementation was inspired by the MATLAB code of John Burkardt[4], which is based on the FORTRAN library of Sylvan Elhay and Jaroslav Kautsky[5],[6].

Additional information on the Gaussian quadrature and numerical integration can be found at [7] and [8].


  1. Golub G., Welsch J., Calculation of Gauss quadrature rules, Math. Соmр. 23 1969.^
  2. Gautschi W., Orthogonal Polynomials: Computation and Approximation. Oxford University Press. 2004.^
  3. Martin R., Wilkinson J., The Implicit QL Algorithm, Numerische Mathematik, Volume 12, Number 5, December 1968, pages 377-383.^
  4. Burkardt, J., Weights for Interpolatory Quadrature. ^
  5. Elhay S., Kautsky J., Calculation of the Weights of Interpolatory Quadratures, Numerische Mathematik, Volume 40, 1982, pages 407-422.^
  6. Elhay S., Kautsky J., Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of Interpolatory Quadrature, ACM Transactions on Mathematical Software, Volume 13, Number 4, December 1987, pages 399-415.^
  7. Davis P., Rabinowitz P., Methods of Numerical Integration, Second Edition, Dover, 2007, ISBN: 0486453391, LC: QA299.3.D28.^
  8. Stroud A., Secrest D., Gaussian Quadrature Formulas, Prentice Hall, 1966, LC: QA299.4G3S7.^

{ 0 comments… add one now }

Leave a Comment

Previous post:

Next post: