[Bug-gsl] Re: GSL Mathieu Function Bug

[Lowell Johnson]

On Tuesday 09 December 2008 10:38:18 am Chad Mitchell wrote:
> Dear Lowell Johnson,
>
> I'm writing to inform you of a small bug appearing in the Special Functions
> segment of the Gnu Scientific Library.  The bug appears in the routines
> "gsl_sf_mathieu_a" and "gsl_sf_mathieu_b" of the file "gsl_sf_mathieu.h",
> which are used to compute the characteristic values of the Mathieu
> functions ce_n and se_n, respectively.
>
> Our group has used Mathieu functions and their characteristic values
> heavily, and the problem is one that we have noticed in other software as
> well.  The attached figure BlipGSL.eps is a plot of the characteristic
> value a_29(q) (associated with the function ce_29) versus the parameter q. 
> Values were computed using the routine gsl_sf_mathieu_a, over the interval
> [0,1200].  The "blips" that appear near q=450, q=625, and q=850 are a
> problematic feature we have seen often.  For example, I have included in
> BlipMathematica.pdf the same characteristic value as computed using
> Mathematica 6.0.  In this case, the "blips" are quite narrow in q, and the
> plot is therefore shown on the subdomain [300,700].  In another code that
> we have investigated, the characteristic values are computed as the zeros
> of a transcendental function, and the appearance of these blips is
> associated with a convergence problem in the root-finding algorithm near
> asymptotes of this function.

Hi Chad,

Thanks for the report.  As you're probably aware, the largest effort involved 
with computing the characteristic values is getting convergence to the 
*correct* characteristic value.  When first implementing the Mathieu functions 
back in the late '80s, I'm sure I spent well over 90% of my time working on 
getting accurate enough initial approximations to keep the root-finding 
algorithm from converging to an adjacent solution.

> Is there a reference or documentation available for the algorithm used by
> gsl_sf_mathieu_a and gsl_sf_mathieu_b?  I am curious to investigate this
> problem further.  We understand the computational challenges posed by the
> Mathieu functions, and I would be eager to communicate more about this
> issue.

No, there isn't much for specific documentation.  I basically pulled the 
algorithms together based on Abramowitz & Stegun, Chapter 20.  The recurrence 
relations are used to create the continued fractions, truncated at order+50 
terms.  This truncated continued fraction equation is solved for the 
characteristic value.

Initial approximations for the char. values are computed using the equations 
given in F. A. Alhargan's paper, "Algorithms for the Computation of All 
Mathieu Functions of Integer Orders," ACM Transactions on Mathematical 
Software, Vol. 26, No. 3, September 2000, pp. 390-407.

But the approximations all too often are not accurate enough to keep the root-
finding process from locking onto the wrong solution.  To overcome this, the 
solution is checked to determine whether it is thought to be the correct 
solution.  The following check is performed:

  if (|a - approx| > 3 + 0.01*order*|a - approx|)
      counter = counter + 1
      if (a > approx)
          approx_new = approx - 0.025*counter
      else
        approx_new = approx + 0.025*counter

where a is the root found, approx is the original approximation, order is the 
order of the Mathieu function, counter is an integer starting at 1 and 
incrementing with each failed iteration and approx_new is the new 
approximation to be used in the next root-finding attempt.

The above is basically a hack to improve the chances of finding the correct 
solution -- and for a hack, it's surprisingly successful. ;)  I'd say that the 
original approximation alone is only successful for a relatively limited range 
of values, whereas with the hack, I wasn't aware of any failures until your 
report (although we haven't tried extremely hard to push the limits).

I'm guessing that the "blips" you've found are associated with holes in the 
above algorithm, which empirical to a large extent.

If you don't mind reading C source code, all of the above can be found in the 
specfunc/mathieu_charv.c GSL file.

Having said all that, have you tried using the functions that compute ranges 
of characteristic values in one shot?  These functions are 
gsl_sf_mathieu_a_array and gsl_sf_mathieu_b_array.  They use the matrix, 
eigenvector approach, finding the eigenvalues of the truncated tri-diagonal 
array created from the recurrence relations for a given value of q. A 
workspace needs to be allocated first, using the function gsl_sf_mathieu_alloc. 
 
The array functions have two features that make them preferred for many (but 
not all) uses:

1) They solve for an array of characteristic values in one shot.  Oftentimes, 
the user wants a series of characteristic values for a given argument, q.  
Rather than calling the single solution function multiple times, one call to 
the array function is all that's needed.

2) Because it solves for a series of characteristic values, the problem of 
jumping to the wrong solution doesn't exist.  I've never known the array 
functions to fail, or have the "blip" issue we're discussing.

So I suggest trying the array functions.  Hopefully, they'll work for your 
needs.  But I do plan to use your report to improve the root-finding procedure. 
 
I'm fairly confident that I can patch the "blips" you've reported (hopefully 
fixing others in the process).  But the equations are pretty fickle -- even 
very 
minor perturbations can throw them off the desired convergence.

Thanks again for your report.  I'll send out an update when I've made some 
progress.  But I can't guarantee a rapid response, since I have my "paying 
job" and family that absorb most of my time. ;)  This is just done for the 
pure enjoyment of contributing something back.

Regards,

    Lowell
-- 
Lowell D. Johnson     address@hidden
Free and Open Source Software:  Of the people, by the people,
    for the people.
Support open standards:  Use software designed to improve communication,
    not hinder it.
Do you Ubuntu?  (www.ubuntu.com)