Please can somebody advise on using gsl_sf_bessel_Jn_array for
small z?
I need to generate Bessel functions up to some n_max (of order 50)
for
arbitrary z. I run into problems for small z because the function
uses
a downward recurrence to populate the array, and suffers underflow
for
large n.
Could you send a small example program showing the problem to
bug-gsl@gnu.org for reference. Thanks.
The following is sufficient to demonstrate the issue.
double besselArray[51];
double rho = 1e-6;
gsl_sf_bessel_Jn_array(0, 50, rho, besselArray);
The problem is that the array is populated by downward recurrence, but
J_51 and J_50 are way smaller than DBL_MIN, resulting in an error
within the GSL library:
gsl: gamma.c:1454: ERROR: underflow
Default GSL error handler invoked.
Because of the downward recurrence, it is _not_ sufficient just to
ignore the error in my own code: as written at the moment, _none_ of
the array entries will be correctly populated if this error condition
occurs.
My current, very rough-and-ready solution looks something like this:
memset(besselArray, 0, sizeof(besselArray));
if (fabs(rho) < 1e-30)
besselArray[0] = 1;
else if (fabs(rho) < 1e-7)
gsl_sf_bessel_Jn_array(0, MIN(n_max, 3), rho, besselArray);
else
gsl_sf_bessel_Jn_array(0, MIN(n_max, 20), rho, besselArray);
... but it's rather heuristic and I'm sure there are more elegant
solutions that would correctly populate all elements of the array that
are >= DBL_MIN, rather than using an arbitrary cutoff.
