[Bug-gsl] GSL Mathieu Function Bug

[Chad Mitchell]

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.

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.

Thanks and best wishes,


Chad Mitchell
DSAT Group, Physics
University of Maryland
College Park, MD

BlipGSL.eps
Description: PostScript document

BlipMathematica.pdf
Description: Adobe PDF document