6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 1function [x, w] = lp_gauss(p)
2% LP_GAUSS Gauss-Legendre nodes and weights on [-1, 1].
3%
4% Newton's method on the Legendre polynomial P_p, starting from the Chebyshev
5% estimate of the roots; three or four iterations reach roundoff for the small
6% p used here. The weights are 2 / ((1 - x^2) P_p'(x)^2).
7x = cos(pi * ((1:p)' - 0.25) / (p + 0.5));
8dp = ones(p, 1);
9for it = 1:50
10 p0 = ones(p, 1);
11 p1 = x;
12 for k = 2:p
13 p2 = ((2*k - 1) * x .* p1 - (k - 1) * p0) / k;
14 p0 = p1;
15 p1 = p2;
16 end
17 dp = p * (x .* p1 - p0) ./ (x.^2 - 1);
18 dx = p1 ./ dp;
19 x = x - dx;
20 if max(abs(dx)) < 1e-15
21 break;
22 end
23end
24w = 2 ./ ((1 - x.^2) .* dp.^2);
25end