[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] [bug #25075] blips in mathieu function

From: Brian Gough
Subject: [Bug-gsl] [bug #25075] blips in mathieu function
Date: Fri, 12 Dec 2008 10:47:53 +0000
User-agent: Emacs-w3m/1.4.4 w3m/0.5.2


                 Summary: blips in mathieu function
                 Project: GNU Scientific Library
            Submitted by: bjg
            Submitted on: Fri 12 Dec 2008 10:47:52 AM GMT
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: 
                  Status: Confirmed
             Assigned to: None
             Open/Closed: Open
                 Release: 1.12
         Discussion Lock: Any



From: Lowell Johnson <address@hidden>
To: "Chad Mitchell" <address@hidden>
Cc: address@hidden, "Alex J. Dragt" <address@hidden>
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:

> 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
of order=29, 0 < q < 1000.  The array functions do indeed avoid the "blips"
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
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

   /* 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)



Bug-gsl mailing list


Reply to this item at:


  Message sent via/by Savannah

reply via email to

[Prev in Thread] Current Thread [Next in Thread]