What am I missing in this highly oscillatory integral?

I want to numerically integrate this equation (in python without calling Mathematica):

$\int_0^\infty {\rm d}k f(k) J_v(r k) J_n(s k)$

where $f(k)$ is a non-smooth function, $J_v$ are the Bessel function of the fist kind, and $r$ and $s$ are just constants. In practice I have to integrate

$\int_a^b {\rm d}k f(k) J_v(r k) J_n(s k)$

with usual values for $a \sim 10^{-8}$ and $b \sim 10^{16}$, where $f(k)→0$ for $b≥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 $∼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∼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"}}]

Clearly 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!

• It has been a while since I last looked at Levin's paper, but in the meantime: have you seen this? – J. M. will be back soon Jul 29 '15 at 15:51
• @J. M. Apropos of Levin, I was against him, and on the losing team, in a doubles ping pong game last month. – Daniel Lichtblau Jul 29 '15 at 16:12
• @Daniel, As in David Levin of Tel Aviv himself? I didn't know he played table tennis… ;) – J. M. will be back soon Jul 29 '15 at 16:22
• @J. M. Yes, that Levin. Basement here – Daniel Lichtblau Jul 29 '15 at 17:00
• Hmm, in that case, apart from Levin's papers, you might consider looking for papers by Arieh Iserles on oscillatory integration. If all else fails, you can always fall back on Longman's method. – J. M. will be back soon Jul 29 '15 at 18:23