help-gsl
[Top][All Lists]

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

 From: Lowell Johnson Subject: Re: [Help-gsl] Different value for mathieu_ce in Mathematica and GSL Date: Sat, 18 Feb 2017 08:59:27 -0600 User-agent: KMail/5.2.3 (Linux/4.4.0-63-generic; KDE/5.28.0; x86_64; ; )

```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

```