[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] bug report on about the Mathieu function
From: |
Lowell Johnson |
Subject: |
Re: [Bug-gsl] bug report on about the Mathieu function |
Date: |
Wed, 23 Jan 2013 21:27:05 -0600 |
User-agent: |
KMail/4.9.5 (Linux/3.2.0-36-generic; KDE/4.9.5; x86_64; ; ) |
On Wednesday, January 23, 2013 3:50:26 PM Takeshi Morita wrote:
> Dear Sirs/Madams,
>
> My name is Takeshi Morita, a Japanese physics researcher.
>
> I found a bug in GSL about the Mathieu characteristic function and Mathieu
> function.
>
> I needed to calculate the Mathieu characteristic function;
> gsl_sf_mathieu_a(r, q, &result)
> gsl_sf_mathieu_b(r, q, &result)
> at large q.
>
> Then I found that
> gsl_sf_mathieu_a(r, q, &result) - gsl_sf_mathieu_b(r, q, &result)
> and/or
> gsl_sf_mathieu_b(r+1, q, &result) - gsl_sf_mathieu_a(r, q, &result)
> becomes NEGATIVE, although they have to be always POSITIVE.
> (See p.724 (a) in the Handbook of Mathematical Functions by Abramowitz and
> Stegun.)
>
> The attached file is a sample code to see this bug.
>
> For example, when q=2000.0,
> gsl_sf_mathieu_b(70, q, &result) - gsl_sf_mathieu_a(69, q, &result)
> becomes -355.55....
>
> I tried several q and r, and this error happenes when q is roughly larger
> than 1024 and r is around 70.
>
> Related to this issue, gsl_sf_mathieu_ce and gsl_sf_mathieu_se
> do not work properly.
>
> I compiled the code by using gcc (version 4.5.3) on Cygwin on Win 7.
> I asked my friend and he also confirmed the same bug on his following PC
> environment; Linux 3.2.0-4-amd64 #1 SMP Debian 3.2.32-1 x86_64
> gcc version 4.7.2 (Debian 4.7.2-5)
> gsl version 1.15
>
> So I think this bug is independent of the PC environment.
>
> I would appreciate it if you could confirm the bug and fix the problem.
>
> Sincerely yours,
>
>
> Takeshi Morita
Hi Takeshi,
My name is Lowell Johnsonn, and I'm the contributor of the GSL Mathieu
function routines. As you may already be aware, computing the characteristic
values from the recurrence relations and continued fractions is quite
unstable, computationally. Root-finding methods can easily converge to the
"wrong" solution (usually a neighbor solution off by one order from the one
desired). The process is dependent on an accurate initial guess for the
solution as well as techniques to determine whether the "correct" solution has
been found. When the desired solution isn't found, the initial guess must be
slightly modified and retried.
An alternative method, which is much more stable, but requires additional
memory, is to find the characterstic values as an array of eigenvalues of the
tridiagonal matrix that represents the truncated series of recurrence
relations. This method finds all characteristic values from order 0 to n for a
given q.
I'm not familiar with all of the applications for Mathieu functions, but from
my experience in elliptical wave equations, I consider the eigenvalue approach
to be preferred. Since I usually need the series of solutions from 0 to n
anyway, it's more efficient to compute and store them all at once. And there's
no concern about computing the wrong solution for a given order. The primary
problem (which hasn't applied to my situation) is if you need orders large
enough to hit memory constraints (i.e., orders >> 1000).
I've attached a modified version of your test routine, adding in the
gsl_sf_mathieu_a_array() (and _b_ version also).
Hopefully, the array solutions will work for your application. If not, we may
be able to add additional tests into the continued fraction functions. But it
is a difficult problem.
If someone else has suggestions and is willing to modify the functions in
mathieu_charv.c, that would also be much appreciated. My available time for
improving the Mathieu functions is very limited at this time, but I do want to
see them continue to be useful and as accurate as possible.
Regards,
--
Lowell D. Johnson
char_val_array.c
Description: Text Data