fix quadrature bug (increased error w/quad_pts)

[nestor98]

Hi!
running some experiments (with the current version) I noticed that sometimes increasing quad_pts in iad.Sample.rt made the quadrature diverge and then settle on some different (wrong) value at higher quad_pts. I fixed the problem (with some help from an LLM...)

The issue

An example small reproduction:

import numpy as np, iadpython as iad
for q in range(2, 70, 2):
    s = iad.Sample(a=0.5, b=1.0, d=1.0, g=0.9, n=1, n_above=1, n_below=1, quad_pts=q)
    print(f"q={q:3d}  T_L={s.rt()[1]}")

This prints:

q=  8  T_L=0.5910431400369376
q= 16  T_L=0.5907416292385642
q= 28  T_L=0.5907126796822936
q= 32  T_L=0.5897418927090614
q= 34  T_L=0.6037212508561103     # starts diverging
q= 36  T_L=0.4949748616350780
q= 38  T_L=0.9387339150820570   
q= 40  T_L=0.4139662134049849
q= 50  T_L=0.3885445116257321
q= 64  T_L=0.3794837664452013     # different wrong convergence
q= 80  T_L=0.3767384295008461

Then, around 100 quad_pts it starts raising ValueError: array must not contain infs or NaNs

Issue 2

A related error happens with radau quadrature (also for higher quad_pts), which returns weights that do not add to one as they should:

from iadpython import quadrature as Q
for n in (16, 32, 34, 40, 64, 80):
    print(n, Q.gauss(n, 0, 1)[1].sum(), Q.radau(n, 0, 1)[1].sum())

n	gauss sum(w)	radau sum(w)
16	1.00000000	1.00000000
32	1.00000000	0.99996484
34	1.00000000	1.00145222
40	1.00000000	0.90883375
64	1.00000000	0.49334245
80	1.00000000	0.45410315

The fix

The problem is pp = scipy.special.legendre(n - 1).deriv(1), where the deriv function evaluates the polynomial, which seems to be very ill conditioned in floating point for high degrees.

The fix is to switch that derivative to an analytical version (Bonnet's) that doesn't have the numerical issue:

$$ P'_m(x) = \frac{m(x P_m(x) - P_{m-1}(x))}{x^2 - 1} $$

This has been patched in quadrature.py, also changing legendre for eval_legendre, which should evaluate faster (didn't actually measure this).

Results

Now the integration from issue 1 converges correctly:

quad_pts	stock T_L	patched T_L
8	0.5910431	0.59104314
16	0.5907416	0.59074163
32	0.5897419	0.59072054
34	0.6037213	0.59071865
40	0.4139662	0.59071979
64	0.3794838	0.59071947
80	0.3767384	0.59071946
160	error	0.59071701