2
$\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$
6
  • $\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, 2015 at 3:27
  • 2
    $\begingroup$ @xzczd just for fun. $\endgroup$
    – Kattern
    May 4, 2015 at 3:34
  • $\begingroup$ You are aware that LegendreP[] is built-in? $\endgroup$ May 4, 2015 at 3:47
  • $\begingroup$ Apropos: here is an answer that uses some heavy machinery to generate Gauss-Legendre quadrature rules. $\endgroup$ May 4, 2015 at 3:50
  • 1
    $\begingroup$ FWIW: here is an even more elaborate procedure. $\endgroup$ May 4, 2015 at 4:01

1 Answer 1

3
$\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.