From: kovarik@mcmail.cis.McMaster.CA (Zdislav V. Kovarik)
Newsgroups: sci.math
Subject: Re: Gaussian quadrature abscissae and weights
Date: 4 Jun 1998 03:18:26 -0400

In article <6l2ama$36g$1@cdn-news.telecom.com.au>,
Jim Maguire  <jmaguire@nibunetw.telstra.com.au> wrote:
>Is there an algorithm which will give these, if you tell
>it what order you want? Or are they always copied from tables?

Yes, there is a procedure. Let us fix the terminology: suppose I want a
Gauss formula

  integral[a to b] (w(x)*f(x)) dx approximated by
  sum[j=1 to n] W_j * f(x_j)

so that the error is 0 as long as f is a polynomial of degree less than
2*n. The most published situation (called Gauss-Legendre integration) is

 a = -1, b = 1, w(x) = 1 for all x

In general, it is assumed that w(x) is measurable, defined and
non-negative almost everywhere on (a,b), positive on a set of positive
measure within (a,b), and such that all integrals over (a,b) of
w(x)*abs(x^n) dx are finite. 

Stage 1: Find an orthogonal polynomial p_n(x) of degree n corresponding to
weight w on the interval (a,b).

 (The scalar product is defined as 
   (p,q)=integral[a to b] w(x)*p(x)*q(x) dx, 
 and if nothing else is known, a way to find p_n  is Gram-Schmidt
 orthogonalization.)

   [Piece of theory: It can be proved that p_n has n distinct (real)
    roots, all in the interval (a,b) .]

Stage 2: Find the roots x_1, ... , x_n  of  p_n. They will be the
abscissas.

Stage 3. Construct auxiliary polynomials

 L_j(x) = p(x)/((x-x_j)*p'(x_j)) , L(x_j)=1

and the weights will be

 W_j = integral[a to b] w(x) * L_j(x) dx

   [Piece of theory: All W_j come out positive. Why? Because also

     W_j = integral[a to b] w(x) * (L_j(x))^2 dx .]

Reference: a good Numerical Analysis (as opposed to Numerical Methods)
textbook. Also, a Special Functions textbook with a chapter on orthogonal
polynomials.

Hope it helps, ZVK (Slavek).