[SciPy-dev] Numerical stability of special functions (here: chebyu)
Norbert Nemec
Norbert.Nemec.list at gmx.de
Mon Apr 24 03:24:03 CDT 2006
Hi there,
I just had my first encounter with Chebyshev polynomials of the second
kind. A lengthy calculation could be solved very elegantly using these
special functions. In the end, I wanted to plot the solution using
scipy&numpy&matplotlib. To my disappointment, the results showed strong
numerical errors, so I had to dig in deeper.
I quickly found out that the call chebyu(N)(x) itself returned bad
results for N>30 and complex x ~= +/- 1.0
I checked the implementation of chebyu and found that it basically first
calculates the coefficients from the generating function and then adds
up the polynomial. In the numerical range of question, this means adding
up huge numbers with alternating sign, producing bad numerical errors.
After a long odyssee of optimizations (I even considered using gmpy to
simply increase the internal precision) I found a ridiculously simple
solution:
----------------
def chebyu(N,x):
previous = 0.0
current = 1.0
for n in range(N):
new = 2*x*current - previous
previous = current
current = new
return current
----------------
proved to be perfectly stable for all my applications. For N>1000 and
|x|>1.1 it starts producing NaNs because of an internal overflow, but I
have not, yet, encountered error accumulation.
Now, it would, of course, make sense to contribute this implementation,
but this is, where I'm lost:
Currently, chebyu is defined in special/orthogonal.py, where it first
creates a orthopoly1d object which can then be evaluated. As mentioned,
this two-step process is intrisically problematic. I believe, the same
applies to many of the other polynomials defined in the same place. On
the other hand, in other applications, the coefficients themselves are
valuable as well and it would be a pity to throw them out completely.
The best idea that I had so far was to replace the function chebyu by a
class definition:
----------------------------
class chebyu(orthopoly1d):
def __init__(self,N):
self.N = N
... same initialization code as before ...
def __call__(self,x):
previous = 0.0
current = 1.0
for n in range(self.N):
new = 2*x*current - previous
previous = current
current = new
return current
---------------------------
however, I'm not sure whether this solution is unproblematic.
Furthermore, it looks like a rather dirty hack to me...
Greetings,
Norbert Nemec
More information about the Scipy-dev
mailing list