ENH: `special.gauss_legendre`: improve accuracy of weights

[susilehtola]

Describe your issue.

The Gauss-Legendre weights generated by scipy are inaccurate. A comparison to SymPy's analytic roots shows that while the accuracy is good for small numbers of nodes, the discrepancies become noticeable for large numbers of nodes.

Our C++ implementation of Gauss-Legendre quadrature in https://github.com/wavefunction91/IntegratorXX/blob/master/include/integratorxx/quadratures/primitive/gausslegendre.hpp maintains an accuracy of 10*eps in the nodes and weights at least up to 920 nodes. In contrast, the discrepancy in scipy grows with the number of nodes, and reaches an error for the weights of 1.9785e-14 with 500 nodes, which is 89 times machine epsilon, and 1.3362e-13 with 1000 nodes or 601 times epsilon.

Reproducing Code Example

from sympy.integrals.quadrature import gauss_legendre
from sympy.core import S
from sympy import sqrt
import numpy
import os
import time
import scipy

# Generate tests with 20 digit precision
ndigits = 20

for order in [10, 30, 50, 100, 300, 500, 1000]:
    x, w = gauss_legendre(order, ndigits)
    xi, wi = scipy.special.roots_legendre(order)
    dx = xi-x
    dw = wi-w

    print(f'{order=} {max(abs(dx))=} {max(abs(dw))=}')

Error message

order=10 max(abs(dx))=8.5408873170492866089e-17 max(abs(dw))=9.9880261667743136522e-16
order=30 max(abs(dx))=1.4661039174140508279e-16 max(abs(dw))=4.0621141449071046232e-15
order=50 max(abs(dx))=1.5321894279588313403e-16 max(abs(dw))=4.9165183798412732599e-15
order=100 max(abs(dx))=1.6209754214900271019e-16 max(abs(dw))=7.0182654149157263719e-15
order=300 max(abs(dx))=1.6611925458256612820e-16 max(abs(dw))=8.4826851919822337760e-15
order=500 max(abs(dx))=1.6198403973407063394e-16 max(abs(dw))=1.9785190888976244499e-14
order=1000 max(abs(dx))=1.6572623129504013284e-16 max(abs(dw))=1.3361781919090518247e-13

SciPy/NumPy/Python version and system information

>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.10.1 1.24.4 sys.version_info(major=3, minor=11, micro=5, releaselevel='final', serial=0)
lapack_armpl_info:
  NOT AVAILABLE
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['flexiblas', 'flexiblas']
    library_dirs = ['/usr/lib64']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['flexiblas', 'flexiblas']
    library_dirs = ['/usr/lib64']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_armpl_info:
  NOT AVAILABLE
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['flexiblas', 'flexiblas']
    library_dirs = ['/usr/lib64']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['flexiblas', 'flexiblas']
    library_dirs = ['/usr/lib64']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

[ev-br]

Pull requests welcome!

[mdhaber]

I don't think we need to pull in another library. The closed-form approximations in this paper get more accurate as the degree increases, so we can just switch algorithms once the degree gets beyond a threshold. I only coded up the first three terms of a series, and it looks like the worst error is 618 ULP with $n=100$, 41 ULP with $n=150$, and 9 ULP with $n=200$. The paper claims that with all the terms, the approach gives roots and weights accurate to machine precision for $n \geq 30$.

```python3 import numpy as np from scipy import special from sympy.integrals.quadrature import gauss_legendre def gauss_legendre2(n): # as written, for even `n` only. # Just need to insert the middle node for odd `n` # Here, I've used SciPy functions for # j0 (roots of Bessel function of order 0) and # J1 (evaluation of Bessel function of order 1) # but asymptotic expansions for these are given in the paper, too. # See section 4. def F0(x, u): # 2.25 return x def F1(x, u): # 2.26 return 1/8 * (u*x - 1)/x def F2(x, u): # 2.27 return 1/384 * (6*x**2 * (1 + u**2) + 25 - u*(31*u**2 + 33)*x**3)/x**3 def W0(x, u): # 3.9 return np.ones_like(x) def W1(x, u): # 3.10 return 1/8 * (u*x + x**2 - 1) / x**2 def W2(x, u): # 3.11 return 1/384 * (81 - 31*u*x - (3 - 6*u**2)*x**2 + 6*u*x**3 - (27 + 84*u**2 + 56*u**4)*x**4) / (x**4) F = [F0, F1, F2] W = [W0, W1, W2] vn = 1 / (n + 0.5) # 2.5 j0 = special.jn_zeros(0, n/2) ank = vn * j0 # 2.12 # 2.24 x = ank u = np.cos(ank)/np.sin(ank) thnk = 0 for m in range(3): thnk += F[m](x, u)*vn**(2*m) xnk = np.cos(thnk) # 2.3 roots = np.concatenate((-xnk, xnk[::-1])) # 3.8 sm = 0 for m in range(3): sm += W[m](x, u)*vn**(2*m) J = special.j1(j0) # just before 2.16 wnk = 2 / (J**2/vn**2 * ank/np.sin(ank) * sm) weights = np.concatenate((wnk, wnk[::-1])) return roots, weights ```

n = 200
x, w = gauss_legendre2(n)
x0, w0 = gauss_legendre(n, 20)
x0, w0 = np.asarray(x0, dtype=np.float64), np.asarray(w0, dtype=np.float64)
np.testing.assert_array_almost_equal_nulp(w, w0)
# AssertionError: Arrays are not equal to 1 ULP (max is 9)

[j-bowhay]

Just a comment to suggest 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 (2013), A652-A672. This is used in Chebfun and is impressively fast

[mdhaber]

Faster than the code above? With everything implemented, it would require no iteration under the hood. The paper behind it cites Hale and Townsend and notes several other advantages.

In any case, this is about the accuracy and what can get done. I would be willing to review a PR based on the Bogaert paper since I'm now familiar with it and there is already code.

[j-bowhay]

That is a good question, I'm away from a computer for some time so can't investigate. Do of course feel free to proceed with what you have proposed, it was just that I remembered that one of Chebfun's party tricks was the computation of guass Legendre weights so thought it was worth throwing this paper into the ring