Re: [Bug-gsl] bug report on about the Mathieu function

[Lowell Johnson]

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