# Radial integration of expensive function with Bessel weights

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.

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

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 this methods (types of methods):

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:

For direction 2, I would advise trying to develop some form of a 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):

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: