Short version

I'm computing the zero-order Hankel transform of a function $f(r)$, $$F(k) = \int_0^\infty f(r)J_0(kr)r\,\mathrm{d}r$$ I know $f(r_n)$ at selected $r_n$ but it is impractical to compute $f$ at additional values of $r$. I need the transformed function $F(k)$ to have the best possible accuracy at large values of $k$, especially up to $k\sim 10^4$ but the higher the better, and I'm looking for algorithms or techniques that might aid this process.


I know a few things about this function on theoretical grounds:

  • $f(0) = 1$ and $f'(0) = 0$
  • $0 < f(r) \leq 1$ for all $r$
  • $f(r)\to 0$ as $r\to\infty$, roughly like $r^{-4}$
  • the function is smooth and has no singularities

Beyond that, what I have to work with are a set of numerical abcissas $r_n$ and corresponding function values $f_n$. In theory the abcissas are supposed to be logarithmically spaced, $r_n = r_0 \delta^n$, but I think some values got mangled along the way and they're a little off from exact logarithmic spacing. Anyway, the point is, it's nearly impossible to evaluate the function at arbitrary values of $r$. For purposes of this question, I have a list of points $(r_n, f_n)$ and that's it. Effectively, the abcissas $r_n$ are $0$ (because I know $f(0) = 1$), $10^{-8}$, and then a set of 5000 values spanning the range from roughly $4\times 10^{-5}$ to $900$, roughly logarithmically spaced. A log-log plot of a representative data set looks like this:

plot of very smoothly varying function

Note that the function is very well-behaved: smooth, slowly varying, high sample rate in the region where it deviates from linearity.

Up to now I've been working with a set of transformed values which I think were generated using something like linear interpolation and Gaussian quadrature with a constant weight function. The original code to do the transformation is missing, which is why I'm not 100% sure if that was the method used. Anyway, some of the computed values for $F(k)$ at $k\gtrsim 1000$ are negative; for example:

k^2      F(k)
1032.56  7.69276e-08
1179.5   -6.37344e-08
1347.35  -2.88371e-07
1539.08  -2.47056e-07
1758.09  -4.05768e-07
2008.28  -2.65087e-07
2294.06  -2.60826e-07
2620.51  -2.93459e-07
2993.42  -2.86159e-08
3419.4   -3.24329e-07
3905.99  5.79954e-08
4461.83  -2.02765e-07

I know, on theoretical grounds, that this cannot be true. Theoretically, $F(k) > 0$ for all $k$. Hence the need for high accuracy.

Are there any specialized algorithms or techniques that can be used to improve the accuracy of $F(k)$, especially at large $k$? Speed is not an issue; the amount of computing power I can throw at this problem is far more than I can imagine it needing.


For large values of $k$, the Hankel transform depends primarily on $f(r)$ at small values of $r$. Your smallest values of $r$ aren't small enough to get good coverage.

For example at $k=1 \times 10^{5}$, and your smallest value of $r=4 \times 10^{-5}$, $kr=4$, so you've basically skipped over the first cycle of the Bessel function.

Since you know that $f(r)$ is smooth, you might try smoothly interpolating values of $f(r)$ and then applying a library routine for computing the Hankel transform from the interpolated values. The accuracy of this approach will depend on how well your interpolated values of $f(r)$ match the actual values of $f(r)$ for small $r$. If your values of $f(r)$ aren't very precise (as you indicate), then this probably won't work well.

You might get a more helpful answer if you gave us more background on your problem.

  • $\begingroup$ Thanks, I suppose I should have been more explicit that the function is nearly constant at small $r$, so I can get a very good approximation using linear interpolation within $r\in[0,4\times 10^{-5}]$. Let me come up with some more information to edit into the question. $\endgroup$
    – David Z
    Dec 7 '15 at 6:17
  • $\begingroup$ Are the values of $f(r)$ measured data (and if so how precise are they) or are they precisely computed values (say accurate to 15 digits)? $\endgroup$ Dec 7 '15 at 19:21
  • $\begingroup$ They're computed values. There is a bit of uncertainty inherent in the computation procedure, but I think you can take them to be accurate to at least 6 digits, maybe more. $\endgroup$
    – David Z
    Dec 7 '15 at 19:40
  • $\begingroup$ Try fitting a cubic spline through the points. If it comes out smooth, then run a library routine for the Hankel transform on the interpolated result. If it's not smooth, then try a smoothed spline. $\endgroup$ Dec 7 '15 at 19:43
  • $\begingroup$ See also marineemlab.ucsd.edu/~kkey/pubs/Geophysics%202012%20Key.pdf You'll need to be able to evaluate f at arbitrary points for either of the approaches discussed in this paper. $\endgroup$ Dec 7 '15 at 22:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.