Integration pkg algo replacement suggestion

[Fotios Kasolis]

In the integration 1.0.7 pkg exists a routine that generates gaussian points and weights (in the file grule.m) using the algorithm given by Davis and Rabinowitz in Methods of numerical integration.

 

I suggest that the algorithm is replaced by Trefethen Golub - Welsch algorithm for the obvious reason! (+ the beauty of the theory!)

 

function [ x, w ] = grule (n)
b = 0.5 ./ sqrt (1 - (2 * (1:n - 1)) .^ - 2);
L = diag (b, 1) + diag (b, - 1);
[ V, D ] = eig (L);
x = diag (D);
[ x, i ] = sort (x);
w = (2 * V(1, i) .^ 2)';

 

/Fotios