Computing the coefficients for Gaussian integral

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?

Is it a programming practice or you just need the Gaussian integral weights? If the latter, search GaussianQuadratureWeights in the document. — xzczd, May 04 '15 at 03:27
Apropos: here is an answer that uses some heavy machinery to generate Gauss-Legendre quadrature rules. — J. M.'s missing motivation, May 04 '15 at 03:50

score 3 · Accepted Answer · edited Apr 13 '17 at 12:19

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]

Computing the coefficients for Gaussian integral

1 Answers1