Severe instability - it is not possible to calculate more than ~10 polynomials

[nvtroshkin]

Example:

def pdf(x):
    return math.exp(-x)

pp = OrthoPoly(pdf, intlims=(0,1))
pp.gen_poly(n)

# a new class field 'betas' has been added to store the beta recurrence coefficients for demonstration
print(pp.betas)

Output:

[1.0,
0.07932640579220766,
0.06711669572587618,
0.06433876970545459,
0.06350596398737984,
0.06313659912953457,
0.06293951687236084,
0.06282180581953442,
0.0627458612867561,
0.06269478868364445,
0.06283407650853606,
0.06089373421606997,
0.21595648578197688,
0.669882586586399,
-1.3444323914957315,
4.962392290085074,
0.4906144952293512,
0.20512208494961423,
-686.035257250946,
15925.79931665891,
56662.12826844703]

The last elements demonstrate instability.

[j-jith]

Sorry for the late response. I just noticed this issue.

Thank you for catching this, @nvtroshkin . I had not tried generating more than 5-6 polynomials. Unfortunately, this is an issue with the Stieltjes procedure used for calculating alpha and beta (the coefficients in the recurrence relations). This procedure is very numerically unstable.

There are other methods like Chebyshev procedure which are more stable. You can refer to [1] for more details. I may implement these methods in the future, But I can't guarantee anything at the moment.

[1] W. Gautschi. On Generating Orthogonal Polynomials. SIAM J. Sci. Stat. Comput., Vol. 3, No. 3, Sep 1982.

PS: You should normalize whatever function you're using so that its integral equals 1. Otherwise, beta_0 should be assigned the value of the integral of your function. This is not going to affect the orthogonality of the polynomials. It's just a convention. I have assumed beta_0 to be 1 because I developed this class for working with probability density functions.