Re: [Help-gsl] Jacobi polynomials

[Francesco Abbate]

2010/5/28 Jonny Taylor <address@hidden>:
> Hi all,
>
> I hope this is not considered too off-topic: my question refers to a 
> third-party extension to GSL (for Jacobi polynomials), but I am hoping the 
> mathematical experts here may be able to point me in the right direction.
> Unfortunately, there are some circumstances when BOTH methods suffer from 
> cancelation e.g. P_34^(-16,16) for x~0.8055. I am therefore rather stuck as 
> to how I can obtain an accurate answer for cases such as this. I notice that 
> Mathematica is able to give an accurate answer in such cases, but I have no 
> idea how to go about calculating that answer myself (I am mathematically out 
> of my depth with this sort of thing). Should I be looking for yet another way 
> of generating the Jacobi function, or do I simply need to resort to higher 
> precision arithmetic in these pathological cases?

Hi,

it seems that your problem is intrinsically difficult to treat because
the Jacobi polynomial with alpha= -16 and beta= 16 is essentially 0
for x=1 and you get huge values for x= -1 (for a given n it will be
(-1)^n*choose(n+b,n)) with oscillatory behaviour. This is also clear
if you look to the weight function which is:

(1-x)^alpha * (1+x)^beta

you can note that this function is very small at one end of the
interval and very big at the other end.

I'm wondering if you really need to evaluate this kind of Jacobi
polynomial or if you could may be adopt another way of performing your
computations.

Having said so, I've been successful with some tests by using the simple formula

gamma(a+n+1)/(gamma(n+1)*gamma(a+b+n+1)) *
\sum_{m=0}^n choose(n,m) * gamma(a+b+n+m+1)/gamma(a+m+1) * ((z-1)/2)^m

you just need to take into account that 1/gamma(a+m+1) is equal to 0
for a+m = -1, -2, -3, ... .

I hope that helps but be aware that I'm not an expert of numerical computations.

Best regards,
Francesco