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

[Phyks]

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