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
_______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
