bugGNU Scientific Library - Bugs: bug #25075, blips in mathieu function

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #25075: blips in mathieu function

Submitter:  -Deleted Account- <bjg>
Submitted:  Fri 12 Dec 2008 10:47:52 AM UTC
   
 
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  Fixed
Assigned to:  bjg Open/Closed:  Closed
Release:  1.12

Wed 14 Jan 2009 03:12:33 PM UTC, comment #2: 

applied in 1c51c985f4981fa5571eecb05d424e52f6a51e69


-Deleted Account- <bjg>
Tue 30 Dec 2008 02:04:21 AM UTC, comment #1: 

From a fairly thorough testing, the attached patch (mathieu.patch) seems to fix the "blips" for 0 <= |order| < 70 for 0 <= q <= 5000.  For orders greater than 69, the initial approximation is set to order^2, which starts to fail for q > 1000.

Unless anyone discovers further problems within the "fixed" ranges listed above, I suggest applying the patch for the next release.

(file #17161)

Lowell D. Johnson <ldj>
Fri 12 Dec 2008 10:47:52 AM UTC, original submission:  

From: Lowell Johnson <ldj00@sio.midco.net>
To: "Chad Mitchell" <chad.cemitch@gmail.com>
Cc: bug-gsl@gnu.org, "Alex J. Dragt" <dragt@physics.umd.edu>
Subject: [Bug-gsl] Re: GSL Mathieu Function Bug
Date: Tue, 9 Dec 2008 22:04:35 -0600

On Tuesday 09 December 2008 8:42:25 pm Lowell Johnson wrote:

>

[snip]

> 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.


Following up on my earlier response, I just tried recreating your problem case
of order=29, 0 < q < 1000.  The array functions do indeed avoid the "blips" at
q~480, 630, and 850.

But the array is truncated too short for accurate results beyond q=50 or so.
So I rebuilt, adding another 20 rows to the recurrence matrix.  That appears
to give accurate results all the way to q=1000 (as far as I computed).  If you
want to give it a shot, apply the following patch to your build:

-------------------------------------------------
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


_____________________________________________
Bug-gsl mailing list
-email is unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl


-Deleted Account- <bjg>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #17161:  mathieu.patch added by ldj (5KiB - text/x-diff - Patch to correct characteristic values associated with wrong orders.)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by ldj (Updated the item)
  •  

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2009-01-14 bjg StatusConfirmed Fixed
        Assigned toNone bjg
        Open/ClosedOpen Closed
    2008-12-30 ldj Attached File- Added mathieu.patch, #17161

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code