[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Hankel transform and Bessel functions zeros

From: Michael Carley
Subject: [Help-gsl] Hankel transform and Bessel functions zeros
Date: Fri, 26 Aug 2011 11:28:14 +0100 (BST)
User-agent: Alpine 2.00 (LNX 1167 2008-08-23)

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':

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:

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 ;


`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:
Home page:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]