[SciPy-user] Numerical stability and Legendre polynomials
yaroslavvb at gmail.com
Fri Sep 9 13:29:18 CDT 2005
Thanks for the insight.
The 3-term recursion indeed does much better.
I did some experiments to measure exactly how much better.
Here's the graph of the number of bits needed to represent maximum
error when you project f(x)=0.5 onto Legendre polynomial basis.
Green graph is the old method using scipy's functions, and blue graph
uses 3-term recursion
And here's just the graph of 3-term's recursion until 250'th order
For some reason the error explodes around 210th order
Here's the code used to generate error values
k = 30
vals = zeros(200, Float)
errors = zeros(k, Float)
for j in range(k):
ff = lambda x: eval_poly_3term(j, x)
coef = (2.*j+1)/2*(integrate.quad(ff, -1, 1))
errors[j] = log(max(max(abs(vals-1)),2**-53))/log(2)+53
> In : for n in range(30):
> ....: pn = scipy.special.legendre(n)
> ....: print pn(-1.0)-((-1)**n)*pn(1.0)
For some reason when I execute the same code I get
I'm using 64 bit arithmetic (floats), scipy version 0.3.2, why do you
get higher accuracy?
bulatov at cs.oregonstate.edu
Oregon State University
More information about the SciPy-user