Computing the coefficients for Gaussian integral [duplicate]

Question

I wrote a small script to find the Gaussian integral points and weights:

FindInstance[{Table[x[i]^j, {j, 0, 2 # + 1}, {i, 0, #}].Table[
       w[i], {i, 0, #}] == 
     Table[Integrate[x^i, {x, -1, 1}], {i, 0, 2 # + 1}], 
    Less @@ Table[x[i], {i, 0, #}]} , 
   Flatten[{Table[x[i], {i, 0, #}], Table[w[i], {i, 0, #}]}]] &[n]

For $n$ is small the script is working:

n=1
{{x[0] -> -(1/Sqrt[3]), x[1] -> 1/Sqrt[3], w[0] -> 1, w[1] -> 1}}

n=2
{{x[0] -> -Sqrt[(3/5)], x[1] -> 0, x[2] -> Sqrt[3/5], w[0] -> 5/9, 
  w[1] -> 8/9, w[2] -> 5/9}}

but for $n\ge5$, the script is extremely slow.

Is it possible to speed up the script? If not, is it possible if I only want numerical results?

Answer 1

If you don't need to many quadrature points (say, less than 1000), then the easiest way to get them is simply by finding the eigenvalues of the tridiagonal Jacobi matrix, the Golub Welsch algorithm.

Here's a post describing the rules for the matrix elements.

pointsAndWeightsGaussLegendre[n_] := Module[{jacobi},
  jacobi = ConstantArray[0, {n, n}];
  Do[jacobi[[i, i + 1]] = 
    jacobi[[i + 1, i]] = i/Sqrt[4.0 i^2 - 1], {i, n - 1}];
  jacobi = Eigensystem@jacobi;
  jacobi[[2]] = (2 #[[1]]^2 &) /@ jacobi[[2]];
  Chop@Transpose@Sort@Transpose@jacobi]

This gives the same results as the OP's code, but much faster. Even 1000 points can be calculated in a fraction of a second.

I know there are plenty of codes out there to give the quadrature points and weights, like here and here, I just find this one to be intuitive and easy to implement. I've used Gauss-Hermite quadrature quite often, and so I had this function already laying around:

pointsAndWeightsGaussHermite[n_] := Module[{jacobi},
  jacobi = ConstantArray[0, {n, n}];
  Do[jacobi[[i, i + 1]] = jacobi[[i + 1, i]] = Sqrt[i/2.0], {i, 
    n - 1}];
  jacobi = Eigensystem@jacobi;
  jacobi[[2]] = (Sqrt[\[Pi]] #[[1]]^2 &) /@ jacobi[[2]];
  Chop@Transpose@Sort@Transpose@jacobi]