help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Phyks
Subject: Re: [Help-gsl] Different value for mathieu_ce in Mathematica and GSL
Date: Mon, 20 Feb 2017 10:47:12 +0100
User-agent: Roundcube Webmail/1.2.3

Hi,

Thanks a lot for all the insights and explanations!

I had thought about the different branch cuts between Mathematica and GSL, but MathieuCharacteristicA[0, -1] gives the same value as GSL, -0.455139. I could not yet find values for which MathieuCharacteristicA and the GSL implementation differ.


By going again through the Mathematica documentation, I found out it references Abramovitz and Stegun as well. So the computation should be in agreement.

Finally, I noticed something in their MathieuC documentation that I missed (http://reference.wolfram.com/language/ref/MathieuC.html), under the "Properties and Relations" paragraph. Doing N[MathieuC[MathieuCharacteristicA[0, -1], -1, 2*Pi/180], 50] indeed gives the same result as GSL. So it was definitely an issue on the Mathematica side.


I am not sure if this is doable, but adding some comment in https://www.gnu.org/software/gsl/manual/html_node/Angular-Mathieu-Functions.html#Angular-Mathieu-Functions docstrings, indicating that the `gsl_mathieu_ce` functions are computed acording to the definition in Abramovitz and Stegun (basically a summarized version of Lowell's message below) might be really helpful!


Thanks again for the help!

Le 2017-02-18 15:59, Lowell Johnson a écrit :
On Saturday, February 18, 2017 2:55:00 PM CST maxgacode wrote:
Il 17/02/2017 23:17, Patrick Alken ha scritto:
>>> 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.

Looking at Abramovitz and Stegun I found the following power serie for
Ce(0,q,z) ( for small |q| ).

Ce(0,q,z) = ( 1/sqrt(2) ) * [ 1 - q * cos(2 z)/2 + q^2 * ((cos(4 z)/32)
- 1/16) +........


for  q= -1 , z = 2 pi / 180

Ce(0,q,z) =~ 1.04 + ....

That is not proving anything but my guess is that GSL implementation
agrees with Abramovitz and Stegun.

Moreover Scilab (using the Mathieu Toolbox from R.Coisson & G. Vernizzi,
Parma University, 2001-2002.)

-->mathieu_ang_ce(0,-1, 2 * %pi / 180 ,1)
  ans  =

     0.9975194

again in agreement with GSL, Specfun and Abramovitz.

The Wolfram site says

"For nonzero q, the Mathieu functions are only periodic in z for certain
values of a. Such characteristic values are given by the Wolfram
Language functions MathieuCharacteristicA[r, q] and
MathieuCharacteristicB[r, q] with r an integer or rational number. These
values are often denoted a_r and b_r. In general, both a_r and b_r are
multivalued functions with very complicated branch cut structures.
Unfortunately,

there is no general agreement on how to define the branch cuts.

As a result, the Wolfram Language's implementation simply picks a
convenient sheet. "


What are the values returned by

MathieuCharacteristicA[0, -1]

Hope this helps

Hi all,

I'm the original author of the GSL Mathieu functions, having converted them from the Fortran library I wrote as part of my graduate thesis back around 1990 (which solves for complex order and argument). The goal was to match the results in Abromowitz & Stegun: the primary source of information for this work. I haven't worked on or with these functions for many years, so I may have some of the following statements a bit off. I'll apologize for any
inaccuracies right up front.

It's difficult to compute the correct solution for a given order and large q, because the continued fraction root-finding process is happy to jump to an adjacent solution. The key seems to be having an accurate-enough initial
guess for the solution.

The best way (that I'm aware of) to ensure the solution is the one you want is to find the eigenvalues of the recurrence relation sparse matrix. This will find all of the solutions up to the size of the matrix; they just need to be sorted by value to get them in order. Since this infinite matrix has to be truncated, the computed eigenvalues aren't necessarily as accurate as we'd like. So we use these computed eigenvalues as the initial guesses to the
root- finding process, where we specify the level of accuracy required.

The vast majority of the time spent on developing the GSL Mathieu functions was applied to tweaking and testing various techniques to improve the initial guesses to the solver and to adjusting the root-finding algorithm to throw away steps that jump too far (which means it has hopped to a different branch).

I hope this information is helpful. Because of the above, I wouldn't be surprised (but would be disappointed) to find that the GSL solutions are incorrect for certain regions of the order/argument space. With the help of Brian Gladman and others (whose names I don't recall), we tested for both orders and arguments up to 5000. But it's time-consuming to test all of the
regions that includes in detail.

Regards,

Lowell

--
Phyks



reply via email to

[Prev in Thread] Current Thread [Next in Thread]