1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Is it a programming practice or you just need the Gaussian integral weights? If the latter, search GaussianQuadratureWeights in the document. $\endgroup$ – xzczd May 4 '15 at 3:27
  • 2
    $\begingroup$ @xzczd just for fun. $\endgroup$ – Kattern May 4 '15 at 3:34
  • $\begingroup$ You are aware that LegendreP[] is built-in? $\endgroup$ – J. M. will be back soon May 4 '15 at 3:47
  • $\begingroup$ Apropos: here is an answer that uses some heavy machinery to generate Gauss-Legendre quadrature rules. $\endgroup$ – J. M. will be back soon May 4 '15 at 3:50
  • 1
    $\begingroup$ FWIW: here is an even more elaborate procedure. $\endgroup$ – J. M. will be back soon May 4 '15 at 4:01
2
$\begingroup$

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.

enter image description here

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]
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.