Radial integration of expensive function with Bessel weights

Question

I need to calculate the integral $$I = \int_0^R f(r)J_n\left(\frac{z_{nm}r}{R}\right)rdr$$

where $J_n$ is the $n^{\mathrm{th}}$ order Bessel functions of the first kind, $z_{nm}$ is its $m^{\mathrm{th}}$ zero and $f(r)$ is a real function which is somewhat similar to $J_n$ (but not the same, it is quite complicated and usually involves terms with $J_n^2$ and sometimes $\exp(J_n)$).

As $f(r)$ is extremely expensive and this integral must be evaluated very many times, I am looking for the best (very fast, but still reasonably accurate) numerical method to solve it. Currently, I am using the trapezoidal rule with 11 points. But I am investigating other methods, such as Clenshaw–Curtis and Gauss–Kronrod (with low order).

But I am wondering if there is a method particularly suited to such integrals, especially given that it is similar to those required to compute Hankel transforms.

Answer 1

For the Hankel transform, one can classify the methods into four major groups:

1. Numerical quadrature-based.
2. Fourier-based ones.
3. Asymptotic expansion of Bessel into sines and cosines.
4. Projection-slice methods.

The following paper gives a nice overview of these methods (types of methods):

• M. J. Cree and P. J. Bones, "Algorithms to numerically evaluate the Hankel transform," Comp. Math. Appl., vol. 26, no. 1, pp. 1–12, Jul. 1993.

According to this paper (p. 3, Sec. 4):

The main disadvantage of the trapezoidal rule is its poor computational efficiency. ... the trapezoidal method was found to be almost always as reliable as any other method tested. ...we found it useful as a benchmark to test more efficient algorithms.

So, two directions are possible:

1. Either find a more efficient numerical quadrature rule.
2. Follow along with the Hankel transform-like direction.

In fact, this integral is very similar to the $n$th order Hankel transform. However, I see as a major complication here that $f(r)$ is not only extremely expensive, but oscillatory, and, potentially, highly oscillatory. In the given formulation, $f(r)$ is almost a black-box with some known properties which are actually bad. A lot of the standard techniques for evaluation Hankel transforms use the fact that $f(r)$ is eventually monotonic.

Because of that, I would look into the direction of Filon-type quadrature, that is used for highly-oscillatory integrals. But I think you would have to use some knowledge about the behaviour of $f(r)$ to make it work. That is regarding the direction 1. The following reference (as well as the one for the next part) will be helpful for general understanding:

• R. Barakat and E. Parshall, "Numerical evaluation of the zero-order Hankel transform using Filon quadrature philosophy," Appl. Math. Lett., vol. 9, no. 5, pp. 21–26, Sep. 1996.

For direction 2, I would advise trying to develop some form of projection-slice methods. I am personally not aware of a ready method of this type developed for your integral.

The following references might be useful (for the Filon quadrature, as well):

• J. Markham and J.-A. Conchello, "Numerical evaluation of Hankel transforms for oscillating functions," J. Opt. Soc. America A, vol. 20, no. 4, pp. 621–630, 2003. Here, the authors discuss 6 methods for numerical calculation of the $0$th order Hankel transform of oscillating functions.

Another interesting approach worth trying (me stumbling upon this paper was the main motivation to write this answer) is to apply Mellin transform. However, it might not be worth it if it required too many evaluations of $f(r)$ in the first place:

• L. Knockaert, "Fast Hankel transform by fast sine and cosine transforms: the Mellin connection," IEEE Trans. Signal Processing, vol. 48, no. 6, pp. 1695–1701, Jun. 2000.