Re: [Help-gsl] Different value for mathieu_ce in Mathematica and GSL

[Phyks]

Hello,

Thanks! I was not sure it was appropriate to open a bug in the bug
tracker until I was convinced it was indeed a bug in GSL.

-- 
Phyks

Le 17/02/2017 à 23:17, Patrick Alken a écrit :
> Hello,
> 
>   I'm not familiar with the mathieu functions, but I entered your report
> into the bug tracker so hopefully someone with more knowledge than I can
> take a look and fix it.
> 
> On 02/17/2017 08:53 AM, Phyks wrote:
>> Hi,
>>
>> I have some code that I prototyped in Mathematica and am now writing
>> in C using GSL, that makes use of Mathieu functions. I have different
>> results between the two of them, and I cannot figure out whether this
>> is a bug in GSL, Mathematica or simply some misunderstanding from my
>> part.
>>
>> I am using `MathieuC` function in latest Mathematica
>> (http://reference.wolfram.com/language/ref/MathieuC.html) which should
>> be the same function as `gsl_sf_mathieu_ce`
>> (https://www.gnu.org/software/gsl/manual/html_node/Angular-Mathieu-Functions.html#Angular-Mathieu-Functions)
>> except that the former one takes a single `a` argument being the
>> characteristic value whereas the GSL 2.3 implementation takes the
>> order `n` and the `q` parameter directly.
>>
>> So, I guess,
>> ```
>>> N[MathieuC[MathieuCharacteristicA[0, -1], -1, 2*Pi/180]]
>> 1.41071
>> ```
>>
>> should be equivalent to
>> ```
>> gsl_sf_mathieu_ce(0, -1.0, 2.0 * M_PI / 180.0)
>> ```
>> which gives a totally different value: 0.99751942347886335.
>>
>>
>> I tried to debug with different values, and the discrepancies between
>> Mathematica and GSL seems to appear only when the `q` parameter (-1.0
>> here) is negative. If I take 1.0 instead, I get values in agreement. I
>> tried to find yet another implementation to debug it, and found Scipy
>> (https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.mathieu_cem.html#scipy.special.mathieu_cem)
>> which relies on Fortran SPECFUN library apparently, and is in
>> agreement with GSL.
>>
>> I am missing something? Thanks!
> 
> 
>