I've tried the answers in similar posts but they don't seem to work. As per title, I need to double integrate a complicated quickly oscillatory function. I've checked and there are no poles, the function is well behaved and falls to 0 quickly. The function has three variables: r, td and k. I want to find a plot in terms of k and integrate r and td. I need to integrate td from 0 to infinity and r from td to infinity (thus the Boole[r > td] and both limits set from 0 to infinity). I've tried different integration methods like QuasiMonteCarlo, which yields some result but with a ton of error, specially for big k, or LevinRule which is the most natural but yields an error and an absurd result (a super big number to the power of a super big number). The error is:
XXX is a Levin function of differential order 72 which exceeds value of \option "MaxOrder" -> 50. Treating XXX as a non-Levin function
Where XXX is a long expression related to my input (but weirdly changed in some places)
The code I'm using (for the time being for a given k) is
F0[td_, r_] := 2 (r^2 - td^2)^2 (r^2 + 6 r + 12);
F1[td_, r_] :=
2 (r^2 - td^2) (-r^2 (r^3 + 4 r^2 + 12 r + 24) +
td^2 (r^3 + 12 r^2 + 60 r + 120));
F2[td_, r_] :=
1/2 (r^4 (r^4 + 4 r^3 + 20 r^2 + 72 r + 144) -
2 td^2 r^2 (r^4 + 12 r^3 + 84 r^2 + 360 r + 720) +
td^4 (r^4 + 20 r^3 + 180 r^2 + 840 r + 1680));
Itdr[td_, r_] :=
Exp[td/2] + Exp[-td/2] + (td^2 - r^2 - 4 r)/(4 r) Exp[-r/2];
Integrand[td_, r_] :=
k^3/(12 \[Pi]) (Exp[-r/2] Cos[k td])/(
r^3 Itdr[td, r]) (SphericalBesselJ[0, k r] F0[td, r] +
SphericalBesselJ[1, k r]/(k r) F1[td, r] +
SphericalBesselJ[2, k r]/(k r)^2 F2[td, r]);
k = 0.1;
limit = Infinity;
NIntegrate[
Boole[rd > td] Integrand[td, rd], {rd, 0, limit}, {td, 0, limit}, Method -> "LevinRule"]
I've also tried changing the limits to be finite (after all, the function drops quickly) but this doesn't work specially well. Any idea on what I should try next? Any help is greatly appreciated.
Edit: some plots of the function. For 'small' k=1/10 the function has this form (for different values of td)
You can see that it converges nice and easy. For higher values of k the oscillation is super fast. For k=100:
It still converges, but the fast oscillation makes it hard to find a reasonable integration