# 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.

# Details

• $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:

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.

• 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. – David Z Dec 7 '15 at 6:17
• 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)? – Brian Borchers Dec 7 '15 at 19:21