[SciPy-user] Orthogonal polynomial evaluation

Andrew York Andrew.G.York+scipy@gmail....
Mon Apr 13 01:39:05 CDT 2009

 I am running a simulation with Scipy that needs to numerically
evaluate Jacobi polynomials. I have been using the function
scipy.special.orthogonal.jacobi() heavily. Evaluations of this
function are >90% of my simulation running time. Much more troubling,
I begin to see unphysical results when I need to evaluate high-order
Jacobi polynomials, and I suspect I am running into the problem warned
against here:


"Warning Evaluating large-order polynomials using these functions can
be numerically unstable."

Another Scipy user seems to have run into a similar problem:


 So, it seems the correct course of action is clear: I should identify
a stable algorithm to numerically evaluate Jacobi polynomials, I
should implement this algorithm in something fast like C, and then
call it from within python. I found some promising information on a
possible algorithm:


But I do not speak FORTRAN. I am also new to python and not yet
familiar with how to integrate a fast low-level language with python,
nor am I especially proficient in any low-level language. I suspect if
numerical evaluation of these functions was faster or more robust,
other scipy users would also benefit. Should I attempt this myself (in
C with weave?)? If I do, would others be interested in the results?

Is this an easy task for anyone else? I would greatly appreciate any
help or guidance.

Thank you for your attention.

More information about the SciPy-user mailing list