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

[Lowell Johnson]

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