[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## Re: [Help-gsl] legendre poly deriv

**From**: |
Patrick Alken |

**Subject**: |
Re: [Help-gsl] legendre poly deriv |

**Date**: |
Fri, 10 Nov 2017 09:24:29 +0100 |

**User-agent**: |
Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.4.0 |

I think the older routine from 1.16 has the same issue. For Legendre
derivatives, I believe the current algorithm divides by sin(theta), so
at the poles there is a singularity, therefore the routine fails for x=1
or x=-1 as you found.
The following paper discusses the problem and proposes an algorithm to
fix it:
http://www.sciencedirect.com/science/article/pii/S1464189500001010
It is not currently implemented in GSL but if you're interested in
working on it I would be glad to accept a patch for this
Patrick
On 11/10/2017 07:09 AM, Li Dong wrote:
>* Hi,*
>
>* First of all, thanks for this great library! This is my first post in this*
>* group. Sorry if there is any rule I fail to follow.*
>
>* I was using gsl 1.16 and recently upgraded to 2.3. I found*
>* gsl_sf_legendre_Plm_deriv_array is replaced by gsl_sf_legendre_deriv_array.*
>* Maybe in most cases, this replacement is just fine. However, I find that in*
>* certain cases there is no way to get a result from current*
>* gsl_sf_legendre_deriv_array.*
>
>* In 1.16, the following code works.*
>* double L[5], DL[5]; // here 5 is just an arbitrary large number to*
>* initialize the array*
>* int lmax = 3, m = 2;*
>* double x=-1.0;*
>* gsl_sf_legendre_Plm_deriv_array (lmax, m, x, L, DL); // If m=1, this line*
>* doesn't work since x=-1.*
>
>* In 2.3, it needs to be written as:*
>* double L[5], DL[5]; \\ 5 is again arbitrary*
>* int lmax = 3;*
>* double x = -1.0;*
>* gsl_sf_legendre_deriv_array (GSL_SF_LEGENDRE_NONE, lmax, x, L, DL); // this*
>* line won't work since we calculate all 0<=m<=lmax, it breaks when*
>* calculates m=1 and x=-1*
>
>* I wonder if there is a way around to implement this?*
>
>* Thanks,*
>* Li*