[SciPy-user] Orthogonal polynomial evaluation

[Andrew York]

 I found an inelegant but satisfactory solution to my problem,
perhaps. I'm more comfortable with MATLAB than FORTRAN, and the MATLAB
algorithm here:

https://people.sc.fsu.edu/~burkardt/m_src/polpak/jacobi_poly.m

was very straightforward to translate into python. Perhaps this
algorithm is similar or equivalent to what Paul Virtanen suggested?
This algorithm seems to run faster and resist explosion longer than
scipy.special.orthogonal.jacobi(). I suppose that's not surprising if
the scipy algorithm is meant to find zeros and only evaluates the
polynomial unstably as an afterthought.

If it's helpful, there are many similar algorithms in MATLAB (and
other languages) for special polynomials here:

https://people.sc.fsu.edu/~burkardt/m_src/polpak/polpak.html


On Mon, Apr 13, 2009 at 8:21 PM, Pauli Virtanen <pav at iki.fi> wrote:
> Mon, 13 Apr 2009 20:09:58 +0000, Pauli Virtanen wrote:
> [clip]
>> I wonder if simply making orthopoly1d evaluate the polynomials using the
>> Newton scheme would fix stability issues... Needs some thinking.
>
> This seems to actually help quite a lot: see the branch at
>
>        http://github.com/pv/scipy-work/tree/ticket/921-orthogonal
>
> For example,
>
> Before (Horner):
>
>>>> import scipy.special as sc
>>>> sc.chebyt(300)(0.999) - np.cos(300*np.arccos(0.999))
> -2.9807351165313966e+117
>
> After (using roots):
>
>>>> import scipy.special as sc
>>>> sc.chebyt(300)(0.999) - np.cos(300*np.arccos(0.999))
> 1.9129142714291447e-13
>
> Specialized evaluation routines could of course be faster; finding the
> zeros of the polynomial only for evaluating it is overkill...
>
> --
> Pauli Virtanen
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>