I want to numerically integrate this equation (in python):
$\int_{0}^{\infty}{\rm d}k f(k) J_v(r k)J_v(s k) $,
where f(k) is a non-smooth function, and $J_v$ are the Bessel function of the fist kind. In practice I have to integrate
$\int_{a}^{b}{\rm d}k f(k) J_v(r k)J_v(s k) $,
with usual values for $a \sim 10^{-8}$ and $b \sim 10^{16}$, where $f(k) \rightarrow 0$ for $b \ge 10^{16}$. I wrote a code based on Levin's paper that partially solves my problem, with this method I can integrate up to 50 without much computational power. A exponential change of variables $k = 10^u$, helps me to integrate up to $\sim 10^2$.
Now, I can compute the integral for larger upper limits with more computational power but the error of the integral dramatically increase when I reach $b \sim 10^3$, pretty far for what I want to do.
I have tested this method in the case where $f(k) = k^{-p}$ and using the analytical solution of this paper and using Mathematica with the following instruction:
NIntegrate[f[k]BesselJ[v,r k]BesselJ[v, s k], {k, a, b}, Method -> {"OscillatorySelection", Method -> "LevinRule", "FourierFiniteRangeMethod" -> {"GlobalAdaptive"}}]����
Somehow Mathematica deals with the integral nicely, much better than my code.
So, my question is: what am I missing? What does Mathetica do that I don't?
Any help or ideas are really appreciated. Thanks in advance!
