Gaussian Quadrature

We are trying to approximate the integral using the following formula:

  Integrate[ f[x], {x,-1,1} ]
  
     =  (Subscript[c, 1]) f[(Subscript[x, 1])]  +  (Subscript[c, 2]) f[(Subscript[x, 2])]  +  . . .  +  (Subscript[c, n]) f[(Subscript[x, n])].
 
 
 What we will do is pick these numbers xi and ci so that the formula is exactly true when f[x] is a polynomial of highest degree.  Since there are 2n unknowns, we may suppose that it is true for polynomials of degree 2n-1.  Thus we need only check it for f[x] = (x^k) for k between 1 and 2n-1.
 
 Lets try this out for n = 2.

n = 2; eqn = 
      Table[(∫_(-1)^1) (x^k) \[DifferentialD]x == 
          (∑_(i = 1)^n) c[i] x[i]^k, {k, 0, 2 n - 1}]

solutions = Solve[eqn]

Lets compare with NIntegrate.

Clear[f]GIntegrate[f_] = 
      (∑_(i = 1)^n) c[i] f[x[i]] /. 
        solutions[[1]]

f[x_] = Exp[-(x^2)]; NIntegrate[f[x], {x, -1, 1}]N[GIntegrate[f]]

Lets try this out for larger n.  Compare this with the table on page 209.

n = 3; eqn = 
      Table[(∫_(-1)^1) (x^k) \[DifferentialD]x == 
          (∑_(i = 1)^n) c[i] x[i]^k, {k, 0, 2 n - 1}]; solutions = NSolve[eqn]; TableForm[
      Chop[Table[{k, x[k], c[k]} /. 
            solutions[[1]], {k, 1, n}]]]

Clear[f]GIntegrate[f_] = 
      (∑_(i = 1)^n) c[i] f[x[i]] /. 
        solutions[[1]]; f[x_] = Exp[-(x^2)]; NIntegrate[f[x], {x, -1, 1}]N[GIntegrate[f]]

If you try this for larger n, the computer takes a very long time.  However, it is possible to figure out the values of (Subscript[x, i]) differently, as they are the zeroes of the so called Legendre Polynomials.  These have the property that LegP[n,x] is of degree n, and
  
  \(TraditionalForm
      \`∫\_\(-1\)\%1\ \ 
        LegP[n, x]\ \ \(x\^m\) dx\)      =     0

whenever m<n.

Actually these polynomials are built in to Mathematica.  (They differ from the ones in the book by a constant.)

Information["LegendreP", LongForm -> False]

Clear[n]TableForm[Table[{n, LegendreP[n, x]}, {n, 0, 10}]]

This table illustrates the defining property of the Legendre Polynomials.

MatrixForm[
      Table[(∫_(-1)^1) (x^i) LegendreP[j, x] 
            \[DifferentialD]x, {i, 0, 5}, {j, 0, 5}]]

Let's calculate the values of xi and ci.  The values of ci are given by a formula that is also in the book.

n = 10; sss = NSolve[LegendreP[n, x] == 0, x]; TableForm[
      Table[x[i] = x /. sss[[i]], {i, 1, 
          n}]]

Table[
      c[i] = (∫_(-1)^1)
            ((∏_(j = 1)^(i - 1))((x - x[j])/(x[i] - x[j])))
               ((∏_(j = i + 1)^n)
                  (( x - x[j])/(x[i] - x[j]))) \[DifferentialD]x, {i, 1, 
        n}]

And let's try it out:

Clear[f]GIntegrate[f_] = (∑_(i = 1)^n) c[i] f[x[i]]

f[x_] = Exp[-(x^2)]; NIntegrate[f[x], {x, -1, 1}]N[GIntegrate[f]]