How to make the NIntegate of two variables as efficient as possible?

Question

I asked my question here (Multivariable NIntegrate gives different value if I integrate it separately), I think the best way to get help is by putting the original code here. So, my question is how can I speed up the running of this code? Since I am going to plot the ContourPlot the first and the third variable, it is almost impossible for my computer to run it.

gEnergy = Interpolation[{{977.893, 108.902}, {606.976, 108.902}, {296.819, 108.902}, {125.051,108.902}, {73.127, 106.281}, {36.842, 102.467}, {17.4871, 92.951}, {7.81995,87.4579}, {3.60275, 86.3988}, {1.92659, 82.2894}, {1.06143, 75.5626}, {0.60247, 70.2362}, {0.246369, 59.9484}, {0.199973, 57.0971}, {0.1884, 26.827}, {0.162314, 22.8975}, {0.127878, 19.5436}, {0.0793737, 16.681}, {0.0555056, 15.1319}, {0.0365685, 13.5605}, {0.0220314, 11.8597}, {0.0121379, 11.0237}, {0.00668718, 11.0237}, {0.00317408, 11.0237}, {0.00159913, 10.8902}, {0.000805655, 10.7584}, {0.000405896, 10.1226},{0.000292428, 9.2951}, {0.000210681, 8.1293}, {0.000161109, 7.10971}, {0.000119583, 5.77969}, {0.0000887606, 4.64159}, {0.000074225, 4.21055}, {0.0000620698, 3.91375},{0.0000503808, 3.77328}, {0.000035231, 3.77328}, {0.0000246369, 3.77328}, {0.0000157547, 3.77328}, {9.213*10^-6, 3.77328}}, InterpolationOrder -> 1];

g = 2;
mPlank = 1.22*10^(19);
vlight = 3*10^8;
x0[m_] = 4.26*10^(12)*m;
c[m_, mp_, d_] = (6*d*m)/(\[Pi]^2 mp);
a[d_, v_] = (v/vlight)/(2 d);
s[m_, mp_, d_, v_] = (\[Pi] d/m)^2 *(Sqrt[1 - (mp/m)^2]/a[d, v]) *((Sinh[2 \[Pi] a[d, v]*c[m, mp, d]])/(Cosh[2 \[Pi] a[d, v]*c[m, mp,d]] - Cos[2 \[Pi] Sqrt[c[m, mp, d] - (a[d, v]*c[m, mp, d])^2]]));
f[m_?NumericQ, mp_?NumericQ, d_?NumericQ] := Module[{st, xf, J},
  st[x_?NumericQ] := (x^(3/2))/(2*vlight^3*Sqrt[\[Pi]]) NIntegrate[s[m, mp, d, v] v^2 E^((-x v^2)/(4*vlight^2)), {v, 0, \[Infinity]}, PrecisionGoal -> 4];
  xf = y /. FindRoot[Log[(63 Sqrt[5])/(32 Sqrt[y] \[Pi]^3)*(g/Sqrt[gEnergy[m/y]])* mPlank*m*st[y]] == y, {y, 20, 15}];
  J = NIntegrate[(Sqrt[gEnergy[m/x]]*st[x])/x^2, {x, xf, x0[m]}, PrecisionGoal -> 4];
  (8.77*10^-11)/J]
 

In the end, I want to do ContourPlot:

ContourPlot[f[m, 0.01, d] == 0.1199, {m, 0.1, 10000}, {d, 0.000001, 1}, ScalingFunctions -> {"Log", "Log"}]

One single point take more than one second to run with two errors (NIntegrate::ncvb, General::stop):

f[1, 0.01, 0.001] // Timing

(*{1.32858, 0.000216154}*)

At this point, I think the errors are not too much important, because the final result is reasonable number (or maybe I am wrong).

I am so sorry that I put my whole code here. I tried a lot to make this code as efficient as I can, however, the best I did was this code. Please note that as I explained here , mixing the two integration in one single integration does not give a correct result.

Answer 1

There are things in Mathematica worth trying, but I doubt it will result in massive speedup.

First, as the integral is quickly decaying, setting the manual cutoff can result in increase of efficiency of adaptive methods. I would integrate it from 0 to something like 100*vlight^2/x, the contributions from up there should be meaningless. Accuracy goal may be worth to be set to infinity, and only precision goal to be tweaked. You can consider using LaplaceTransform instead of writing the integral manually to encourage engine to tweak the method of integration to something suitable for this weight ($e^{-xv^2}$).

Another issue I see here is in the cosine (Cos[2 \[Pi] Sqrt[c[m, mp, d] - (a[d, v]*c[m, mp, d])^2]]), but my idea for optimization is rather not mathematica related. To understand this thing, it's better to write it as Cos[2 \[Pi] Sqrt[ - (a[d, v]*c[m, mp, d])^2 - c[m, mp, d]]]. You may then write the denominator of your integrand as

    (Cosh[2 \[Pi] a[d, v]*c[m, mp, d]] - 
    Cosh[2 \[Pi] Sqrt[(a[d, v]*c[m, mp, d])^2-c[m, mp, d]]])

which allows you to either do some perturbation theory if that's an option, or use hyperbolic identities to simplify it further

Now instead of nasty sum you have everything in terms of product of Sinhes (with tiny bit horrible arguments.

So, the thing in your integrand (skipping weight) is now $$ \frac{\sinh(R)}{2\sinh(R/2+\sqrt{R^2 - d}/2)\sinh(R/2-\sqrt{R^2 - d}/2)} $$ with $R=$ a[d, v]*c[m, mp, d] But top part is $$ \sinh(R) = \sinh((R/2-\sqrt{R^2 + d}/2)+(R/2-\sqrt{R^2 - d}/2)) =\\ =\sinh(R/2-\sqrt{R^2 - d}/2)\cosh((R/2+\sqrt{R^2 - d}/2))+ \\ + \sinh(R/2+\sqrt{R^2 - d}/2)\cosh((R/2-\sqrt{R^2 - d}/2)) $$ That simplifies with denominator giving you two hyperbolic cotangents.

I'm not sure what problem is that but I guess you see it's physical interpretation. Mathematica should be able to do that more effectively, but you can go further and now comes the really cool stuff.

Your integral is even, which means you can write it as half of the same integrand integrated from -inf to inf. Integrand is also quickly decaying (due to exponential prefactor). Therefore you can close the contour in the complex plane. Now you are able to find simple poles of the integrand (not terribly hard, these are just coshes with a little mashup) and via Cauchy residue theorem find the answer to your integral as a simple (infinite) sum.

Each term of the sum would be weighted by your exponential prefactor. It's very likely you will be able to evaluate huge part of that sum analytically (or maybe not) but computing infinite sum, especially with the handy exponential-decay prefactor is much faster than adaptive numerical integration.