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
