Integration of a complicated oscillatory function

Question

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

Answer 1

I tried the following:

1. Change the definition of your integrand function to use Set rather than SetDelayed, i.e. Integrand[td_, r_] = ...; there is no need for that static expression to be continuously re-evaluated.
2. Rather than with Boole, write your limits of integration explicitly: NIntegrate[Integrand[td, rd], {td, 0, limit}, {rd, td, limit}].
3. Try the default automatically selected methods before LevinRule which may be slower.

With those in mind, and using your definitions for the rest of the function:

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 = 1/10;

Table[
  {limit, NIntegrate[Integrand[td, rd], {td, 0, limit}, {rd, td, limit}]}, 
  {limit, {1, 5, 100, 200, 300, 600, 1000}}
]

As you mentioned, the integrand appears pretty well-behaved and converges quickly, but higher values of limit make the calculation painfully slow. Perhaps you don't need them though:

ListLogPlot[
  results,
  PlotRange -> All, Joined -> True,
  Mesh -> All, MeshStyle -> Directive[PointSize[0.02], Red]
]

---

It is typically preferable to avoid names starting with capital letters for user-defined quantities (e.g. Integrand) to clearly distinguish them from the built-ins, which are all capitalized.