[Help-gsl] Hankel transform and Bessel functions zeros

[Michael Carley]

I had a query about the Discrete Hankel Transform a while ago which turned 
out to be an error on my side rather than in GSL but investigating it has 
turned up an issue.

There is a loss of accuracy in computing the zeros of Bessel functions in 
gsl_sf_bessel_zero_Jnu. Running the code at the end of this email to 
estimate the error in the eighth order Bessel function at the computed 
zeros gives:

i  log10(fabs(Jn))

1 -14.69
2 -14.75
3 -14.81
4 -15.42
5 -8.16
6 -14.67
7 -14.89
8 -15.04
9 -15.61
10 -15.54
11 -14.20
12 -14.36
13 -14.60
14 -14.67
15 -14.70

and you can see that the Bessel function at the fifth zero is about 10^-8, 
rather than 10^-15. There are similar results for different orders of 
Bessel function. Looking at the code for computing the zeros, it seems 
that this is the point where the method for estimating the zeros changes. 
There is a similar (though not so large) variation in the error at the 
higher roots.

If I use the method of `A Fast Algorithm for the Calculation of the Roots 
of Special Functions':

http://dx.doi.org/10.1137/06067016X

I get much better results with a roughly constant error and if I use this 
method to compute the Bessel function zeros for the Hankel transform, that 
is also improved. I had already implemented the root finding function in 
my Gaussian Quadrature library:

http://www.paraffinalia.co.uk/Software/gqr.shtml

It seems to me that it might be useful to include an implementation of the 
root finding algorithm in GSL since it could be used for a wide range of 
special functions, not just Bessel functions, and would also improve the 
accuracy of the Hankel transform, and I would be happy to recode my 
implementation to GSL standards or pass it on to someone else who could do 
so.

/*test code for zeros of Bessel functions*/

#include <math.h>
#include <stdio.h>

#include <gsl/gsl_sf_bessel.h>

int main(int argc, char **argv)

{
  double x0, Jn ;
  int n = 8, i ;

  for ( i = 1 ; i < 16 ; i ++ ) {
    x0 = gsl_sf_bessel_zero_Jnu(n, i) ;
    Jn = gsl_sf_bessel_Jn(n, x0) ;
    fprintf(stdout, "%d %1.2f\n", i, log10(fabs(Jn))) ;
  }

  return 0 ;
}

gives

--
`To tell the truth, let us be honest at least, it is some considerable
          time since I last knew what I was talking about.'
No MS attachments: http://www.gnu.org/philosophy/no-word-attachments.html
Home page: http://people.bath.ac.uk/ensmjc/