
From:  Phyks 
Subject:  Re: [Helpgsl] Different value for mathieu_ce in Mathematica and GSL 
Date:  Mon, 20 Feb 2017 10:47:12 +0100 
Useragent:  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/AngularMathieuFunctions.html#AngularMathieuFunctions 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 20170218 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, 20012002.) >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 certainvalues of a. Such characteristic values are given by the Wolfram Language functions MathieuCharacteristicA[r, q] andMathieuCharacteristicB[r, q] with r an integer or rational number. Thesevalues 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 helpsHi 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 anyinaccuracies right up front.It's difficult to compute the correct solution for a given order and large q, because the continued fraction rootfinding process is happy to jump to an adjacent solution. The key seems to be having an accurateenough initialguess 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 theroot 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 rootfinding 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 timeconsuming to test all of theregions that includes in detail. Regards, Lowell
 Phyks
[Prev in Thread]  Current Thread  [Next in Thread] 