[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Re: GSL Mathieu Function Bug
[Bug-gsl] Re: GSL Mathieu Function Bug
Tue, 9 Dec 2008 20:42:25 -0600
KMail/1.10.3 (Linux/2.6.27-10-generic; KDE/4.1.3; x86_64; ; )
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.
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
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
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
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
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.
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)