[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Re: GSL Mathieu Function Bug
|
From: |
Lowell Johnson |
|
Subject: |
[Bug-gsl] Re: GSL Mathieu Function Bug |
|
Date: |
Thu, 11 Dec 2008 21:57:54 -0600 |
|
User-agent: |
KMail/1.10.3 (Linux/2.6.27-10-generic; KDE/4.1.3; x86_64; ; ) |
On Tuesday 09 December 2008 8:42:25 pm Lowell Johnson wrote:
> 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.
[snip]
> 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.
Ok, I've done a little more work on this. Turns out the initial
approximations get a little too far off the mark for the regions associated
with the "blips" you're seeing. Normally, my "hack" (given above) nudges the
approximation in the right direction. But for this case of order = 29 (and
others that are similar), the check |a - approx| > 3 + 0.01*order*|a - approx|
is false because the approx value is 1010 and the solution it's finding is
around 1084. This means it thinks it has the correct solution.
So I've implemented a secondary check: |a - approx| > 60. This, combined
with allowing a few more search iterations, has removed the "blips" in the
plot you submitted. I've attached a plot of the results. Note that there's
a chance that adding this secondary check may cause unnecessary iterations in
cases where the approximation is greater than 60 away from the correct
solution (which is found anyway). But at this point, I don't know how likely
that is.
The problem is that there are so many "fringe" cases where the existing
approximations and my secondary checks might fail, but I don't have the time
right now to extensively test a large number of order/argument combinations.
So if someone is willing and able to help in that capacity, I think we could
knock out most of these cases (or at least get a feel for where things still
break down).
An even better alternative would be if someone can improve on the existing
approximations for certain ranges of order and argument. I'm not aware of any
improvements or extensions to those found in Alhargan's paper (given above),
but then I no longer work in a field that keeps me up on the latest papers. So
if someone comes across better approximations for the characteristic values in
these "middle ground" a/q ranges, that would be great.
In the meantime, here's an updated patch for the testing changes I've made:
---------------------------------------------------------------------
diff --git a/specfunc/mathieu_charv.c b/specfunc/mathieu_charv.c
index d059cdf..131280e 100644
--- a/specfunc/mathieu_charv.c
+++ b/specfunc/mathieu_charv.c
@@ -369,7 +369,7 @@ static double approx_s(int order, double qq)
int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result *result)
{
- int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
+ int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;
double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
@@ -437,7 +437,7 @@ int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result
*result)
result->err = GSL_DBL_EPSILON;
break;
}
- if (ii > 20)
+ if (ii > 40)
{
result->err = dela;
break;
@@ -448,7 +448,8 @@ int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result
*result)
/* If the solution found is not near the original approximation,
tweak the approximate value, and try again. */
- if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
+ if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||
+ fabs(aa - aa_orig) > 60)
{
counter++;
if (counter == maxcount)
@@ -482,7 +483,7 @@ int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result
*result)
int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result *result)
{
- int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
+ int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;
double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
@@ -556,7 +557,7 @@ int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result
*result)
result->err = GSL_DBL_EPSILON;
break;
}
- if (ii > 20)
+ if (ii > 40)
{
result->err = dela;
break;
@@ -567,7 +568,8 @@ int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result
*result)
/* If the solution found is not near the original approximation,
tweak the approximate value, and try again. */
- if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
+ if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||
+ fabs(aa - aa_orig) > 60)
{
counter++;
if (counter == maxcount)
diff --git a/specfunc/mathieu_workspace.c b/specfunc/mathieu_workspace.c
index 782c7dd..e6a4615 100644
--- a/specfunc/mathieu_workspace.c
+++ b/specfunc/mathieu_workspace.c
@@ -36,6 +36,7 @@ gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t
nn,
/* Compute the maximum number of extra terms required for 10^-18 root
accuracy for a given value of q (contributed by Brian Gladman). */
extra_values = (int)(2.1*pow(fabs(qq), 0.37)) + 9;
+ extra_values += 20; /* additional fudge */
if (nn + 1 == 0)
{
---------------------------------------------------------------------
Regards,
--
Lowell
a_29.eps
Description: image/eps