Re: [Help-gsl] Underflow in gsl_sf_bessel_Jn_array

[Jonny Taylor]

> > 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
> address@hidden 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.