How to improve precision of NIntegrate while integrating a simple 2D Gaussian distribution?

Question

I have the following expression to integrate

$$ Overlap=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}dt_{1}dt_{2} \left(\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{1}-t_{01})^2}}{\sqrt[4]{\pi }}\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{2}-t_{02})^2}}{\sqrt[4]{\pi }}+\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{2}-t_{01})^2}}{\sqrt[4]{\pi }}\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{1}-t_{02})^2}}{\sqrt[4]{\pi }}\right)^2 $$

This is just two Gaussian distributions and I think there is no problem for Mathematica, I have known that when $t_{01}=t_{02}$ the integral result is 4 and when $t_{02}-t_{01}\rightarrow\infty$ the integral approaches to 2. If we set $t_{01}=0$ and changes $t_{02}$ the ideal integral results would be as follows

However, If I directly integrate the expression using NIntegrate

f\[Tau]gaussian[\[CapitalGamma]_,t_,t0_:5]:=(E^(-(1/2) (t-t0)^2 \[CapitalGamma]^2) Sqrt[\[CapitalGamma]])/\[Pi]^(1/4);

ListPlot[Thread[{Table[t02, {t02, 0, 20, 1}], 
   Table[NIntegrate[(f\[Tau]gaussian[1, t1, 0]*
         f\[Tau]gaussian[1, t2, t02] + 
        f\[Tau]gaussian[1, t2, 0]*
         f\[Tau]gaussian[1, t1, t02])^2, {t1, -Infinity, Infinity}
        , {t2, -Infinity, Infinity}], {t02, 0, 20, 1}]}], 
 PlotRange -> Full, Joined -> True, Frame -> True,
 FrameLabel -> {"t02-t01", "Overlap"}]

when $t_{02}-t_{01}>9$ the results turn to be incorrect. The expression to be integrated could expand into 4 terms and If I integrate each of them independently the results would be much better

nor1 = Table[
   NIntegrate[(f\[Tau]gaussian[1, t1, 0]*
       f\[Tau]gaussian[1, t2, t02])^2, {t1, -Infinity, Infinity}
       , {t2, -Infinity, Infinity}], {t02, 0, 20, 0.5}];
nor2 = Table[
   NIntegrate[(f\[Tau]gaussian[1, t2, 0]*
       f\[Tau]gaussian[1, t1, t02])^2, {t1, -Infinity, Infinity}
       , {t2, -Infinity, Infinity}], {t02, 0, 20, 0.5}];
nor3 = Table[
   NIntegrate[(f\[Tau]gaussian[1, t1, 0]*f\[Tau]gaussian[1, t2, t02]*
      f\[Tau]gaussian[1, t2, 0]*
      f\[Tau]gaussian[1, t1, t02]), {t1, -Infinity, Infinity}
       , {t2, -Infinity, Infinity}], {t02, 0, 20, 0.5}];
nor4 = Table[
   NIntegrate[(f\[Tau]gaussian[1, t2, 0]*f\[Tau]gaussian[1, t1, t02]*
      f\[Tau]gaussian[1, t1, 0]*
      f\[Tau]gaussian[1, t2, t02]), {t1, -Infinity, Infinity}
       , {t2, -Infinity, Infinity}], {t02, 0, 20, 0.5}];
t0mat = Table[t02, {t02, 0, 20, 0.5}];
ListPlot[{Thread[{t0mat, nor1}], Thread[{t0mat, nor2}], 
  Thread[{t0mat, nor3}], Thread[{t0mat, nor4}],
  Thread[{t0mat, nor1 + nor2 + nor3 + nor4}]}, PlotRange -> Full, 
 Joined -> True, Frame -> True,
 PlotStyle -> {Dashed, Dashed, Dashed, Dashed, Thick},
 FrameLabel -> {"t02-t01", "Overlap"}, 
 PlotLegends -> {"nor1", "nor1", "nor3", "nor4", 
   "nor1+nor2+nor3+nor4"}]

from the integration result for each term, we can see that the integral for the expression

$$ Overlap=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}dt_{1}dt_{2} \left(\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{1}-t_{01})^2}}{\sqrt[4]{\pi }}\frac{\sqrt{\Gamma } (e^{-\frac{1}{2} \Gamma ^2 (t_{2}-t_{02})^2}}{\sqrt[4]{\pi }}\right)^2 $$

will be incorrect for certain $t_{02}-t_{01}$ and then will affect the final result. I think this is just a very simple 2D Gaussian distribution and should be integrated easily. I don't know how to use Mathematica to do such kinds of integral and obtain relatively correct results. Can anyone help me with this?

Answer 1

Try adding the modes of the Gaussians to the integration intervals. The problem looks like undersampling. The integrand probably underflows over a large portion of the domain, so you’ve got to help NIntegrate find its numeric support. It should be an otherwise easy integral.

fτgaussian[Γ_, t_, 
   t0_ : 
    5] := (E^(-(1/2) (t - 
           t0)^2 Γ^2) Sqrt[Γ])/π^(1/4);

ListPlot[Thread[{Table[t02, {t02, 0, 20, 1}], Table[
    NIntegrate[(fτgaussian[1, t1, 0]*
         fτgaussian[1, t2, t02] + 
        fτgaussian[1, t2, 0]*fτgaussian[1, t1, t02])^2,
     {t1, -Infinity, 0, t02, Infinity},
     {t2, -Infinity, 0, t02, Infinity},
     Method -> {"MultidimensionalRule", "Generators" -> 9}],
    {t02, 0, 20, 1}]}],
 PlotRange -> Full, Joined -> True, Frame -> True, 
 FrameLabel -> {"t02-t01", "Overlap"}]

Answer 2

I don't know good numeric methods for improper multiple integrals. My best is

f\[Tau]gaussian[\[CapitalGamma]_, t_,t0_] := (E^(-(1/2) (t - 
       t0)^2 \[CapitalGamma]^2) Sqrt[\[CapitalGamma]])/\[Pi]^(1/
  4);

ListPlot[Thread[{Table[t02, {t02, 0, 20, 1}], Table[NIntegrate[(f\[Tau]gaussian[1, t1, 0]*
     f\[Tau]gaussian[1, t2, t02] + 
    f\[Tau]gaussian[1, t2, 0]*
     f\[Tau]gaussian[1, t1, t02])^2, {t1, -Infinity, 
  Infinity}, {t2, -Infinity, Infinity}, 
 Method -> {"AdaptiveQuasiMonteCarlo", "BisectionDithering" -> 0},
  PrecisionGoal -> 2, AccuracyGoal -> 3, MaxRecursion -> 100, 
 WorkingPrecision -> 20], {t02, 0, 20, 1}]}], PlotRange -> Full, 
 Joined -> True, Frame -> True, FrameLabel -> {"t02-t01", "Overlap"}]

, not so bad when comparing with the exact result

ListPlot[Thread[{Table[t02, {t02, 0, 20, 1}], 
   Table[Integrate[(f\[Tau]gaussian[1, t1, 0]*f\[Tau]gaussian[1, t2, t02] +   f\[Tau]gaussian[1, t2, 0]*f\[Tau]gaussian[1, t1, t02])^2, {t1, -Infinity, Infinity}, {t2, -Infinity, 
  Infinity}], {t02, 0, 20, 1}]}], PlotRange -> Full, 
 Joined -> True, Frame -> True, FrameLabel -> {"t02-t01", "Overlap"}]