Hi all,
This RFC proposes adding support for efficiently computing nodes and weights of high-order Gauss–Legendre quadrature rules in SciPy. Currently, nodes are computed by a Golub-Welsch algorithm: finding the eigenvalues of a tri-diagonal matrix and performing one application of Newton’s method which can become computationally intensive as the quadrature order increases. The weights are then calculated from the nodes through a simple formula.
I propose augmenting this method to use an asymptotic expansion to provide an initial guess for the roots of high-order Legendre polynomials, then conduct 1-2 rounds of Newton to converge to the values of the roots. This replaces the eigenvalue computation with a closed-form asymptotic estimate and improves the runtime for large quadrature order while still maintaining accuracy in the nodes and weights. The paper I am pulling from is cited below.
A similar approach is already used in SciPy for Hermite polynomials by automatically applying asymptotic approximations based on parabolic cylinder functions for degrees n \geq 150. I believe it would not be outside the scope of the package to add this feature for Legendre polynomials.
Please let me know what you think and if this would be a good addition to SciPy. Thank you for your time!
N. Hale and A. Townsend. Fast and accurate computation of gauss–legendre and gauss–jacobi quadrature nodes and weights. SIAM Journal on Scientific Computing, 35(2):A652–A674, 2013. [SIAM Link]