URL:
<http://savannah.gnu.org/bugs/?25075>
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
_______________________________________________________
Details:
From: Lowell Johnson <[email protected]>
To: "Chad Mitchell" <[email protected]>
Cc: [email protected], "Alex J. Dragt" <[email protected]>
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 protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?25075>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
