If for some weighting function it is known not only the system of orthogonal polynomials, but also good initial approximations for the polynomial roots, nodes of N-point quadrature Gauss formula could be found as roots of polynomial p_N (x) using Newton method. There is also formula which calculates weights if nodes of quadrature formula are known. It is this method that is used by many subroutine from this section.

Algorithm on the basis of three-term recurrence

If good initial approximations for the polynomial roots are unknown, the method described above is inapplicable: Newton method iteration may not found all the roots. In this case it is reasonable to use method which requires only coefficients of three-term recurrence which generates a system of orthogonal polynomials:

Such recurrence exists for every orthogonal polynomial system, that is, every weighting function W(x) has its own set of coefficients a_j  and b_j  (at that, b₀  isn't used and could be arbitrary). It should be noted that such a recurrence generates a system of orthogonal polynomials having leading coefficient equals 1. But then, it is often used another normalization method which makes an orthogonal system to be an orthonormal (polynomial scalar square equals 1). These two methods and corresponding coefficients a and b should be distinguished.

If coefficients a_i  and b_j  for weighting function W(x) are known, nodes of Gauss quadrature formula could be found as eigenvalues of tridiagonal symmetric matrix

and corresponding weights are calculated by formula , where v_j  is a j-th eigenvector and μ₀  equals . At that, nodes of quadrature formula are separated and located in [a, b], and weights are strictly positive.

Implementation of algorithm

In fact, the algorithm given here is only the wrapper to the algorithm of finding eigenvalues and eigenvectors of a symmetric tridiagonal matrix. On the basis of passed arrays Alpha and Beta (they correspond to a_i  and b_j  coefficients) algorithm forms a tridiagonal matrix J_n , finds its eigenvalues and eigenvectors and then calculates the result - nodes and weights of quadrature formula. At that, there is a special operating mode in the subroutine which finds eigenvectors. The algorithm can calculate only the first row of the eigenvectors matrix. It will take only O(N ²) instead of O(N ³).

Report a bug

C#

C# 1.0 source.
gqgengauss.csharp.zip - Computation of Gauss quadrature rule nodes and weights

C++

C++ source.
gqgengauss.cpp.zip - Computation of Gauss quadrature rule nodes and weights
ablas.zip - optimized basic linear algebra subroutines with SSE2 support (for C++ sources only)

C++, multiple precision arithmetic

C++ source. MPFR/GMP is used.
GMP source is available from gmplib.org. MPFR source is available from www.mpfr.org.
gqgengauss.mpfr.zip - Computation of Gauss quadrature rule nodes and weights
mpfr.zip - precompiled Win32 MPFR/GMP binaries

Delphi

Delphi source.
Can be compiled under FPC (in Delphi compatibility mode).
gqgengauss.delphi.zip - Computation of Gauss quadrature rule nodes and weights

Visual Basic 6

Visual Basic 6 source.
gqgengauss.vb6.zip - Computation of Gauss quadrature rule nodes and weights

Zonnon ^beta

Zonnon source.
Zonnon is an experimental language developed at ETH Zurich.
See www.zonnon.ethz.ch for more information.
gqgengauss.zonnon.zip - Computation of Gauss quadrature rule nodes and weights