It is not possible to calculate the recurrence coefficients accurately with OrthogonalPolynomial

[nvtroshkin]

OrthogonalPolynomial shows not stable behaviour when increasing ncap. As a result it is not possible to calculate the recurrence coefficients accurately (around 10^-2 at most in my case), because there is a maximum ncap from which the numerical error starts dominating.

Presumably, the cause is the lancz method used by default _orthogonal_polynomial.py:102, which is not that stable as Gautschi believes. Probably, it is better to replace it with the sti method (Stieltjes).

Example

Let us consider example 4.1 from Gautschi, 1994:

import orthpol as orth
import numpy as np
import math

ncap=10000

p = orth.OrthogonalPolynomial(20, left=-1, right=1, wf=np.vectorize(lambda x: 1/ncap), ncap=ncap)

for beta in p.beta:
    print(beta)

for ncap=10000 it gives good output:

0.014142135623730923
0.5773502691896286
0.5163977794943228
0.5070925528371096
0.5039526306789716
0.5025189076296065
0.5017452060042544
0.5012804118276035
0.5009794328681206
0.5007733956671929
0.5006261743217585
0.500517330712618
0.5004345937369759
0.5003702332976762
0.5003191829243054
0.5002780094738031
0.5002443195845779
0.500216403386026
0.5001930129390563
0.5001732201680237
0.5001563232803555

but if we increase ncap to ncap=100000, the results are completely different:

0.004472135954986121
0.6613645116574611
0.4416171241289155
0.6373692004945412
0.3178815453321839
0.7808459694960573
nan
2226.7477418947597
3276.4070621703863
1045.0788877186317
1193.8779228463982
nan
243.83832924472938
nan
nan
0.0033354296191723124
nan
nan
2.654608019718409e-06
nan
2.0914153814406423e-08

Comparison of _orthpol.dlancz and _orthpol.dssti

Let us compare Gautschi's dlancz and dsti (the same example).

_orthpol.dlancz:

import numpy as np
import orthpol._orthpol as orth

n=20
ncap = 100000
quad = orth.dfejer(ncap)

for beta in orth.dlancz(n,quad[0], np.array(quad[1]))[1]:
    print(beta)

Output:

1.9999999999879625
0.437403017279912
0.19502568432389394
0.40623949773905066
0.10104867686277746
0.6097204280781962
-0.05588880446006792
4958362.216046399
10734776.222417397
1092180.8409676065
1425321.9460764646
-341783.1687876356
29920.67015771357
-130235.64486376994
-223691960.9528409
1.1114319445532937e-05
-3.464481670371827e-05
-1.4102366473178456e-09
7.053745347898197e-12
-8.787640567426355e-17

_orthpol.dsti:

for beta in orth.dsti(n,quad[0], np.array(quad[1]))[1]:
    print(beta)

Output:

1.9999999999879625
0.3333333333302986
0.2666666666636578
0.2571428571405676
0.25396825396539363
0.2525252525230131
0.25174825174540916
0.2512820512798144
0.25098039215404
0.25077399380583304
0.25062656641321385
0.25051759834144605
0.2504347826058571
0.2503703703681423
0.2503192847992098
0.2502780867608271
0.2502443792737972
0.25021645021421135
0.2501930501902173
0.25017325017100217

Conclusion

The Stieltjes method is more stable in contrast to the Lanczos method used by OrthogonalPolynomial.