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

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

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

applied in 1c51c985f4981fa5571eecb05d424e52f6a51e69

-Deleted Account- <bjg>In charge of this item.
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
-unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl

-Deleted Account- <bjg>In charge of this item.

 

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
  • -unavailable- added by ldj (Updated the item)
  •  

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 4 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 14 Jan 2009 03:12:33 PM UTCbjgStatusConfirmed=>Fixed
      Assigned toNone=>bjg
      Open/ClosedOpen=>Closed
    Tue 30 Dec 2008 02:04:21 AM UTCldjAttached File-=>Added mathieu.patch, #17161

    Back to the top


    Powered by Savane 3.1-cleanup1