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

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);

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}];
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?

• Try adding the modes to the integration intervals: {t1, -Infinity, 0, t02, Infinity} and the same for t2. 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. Oct 31 '21 at 23:36

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"}]


• The documentation to "MultidimensionalRule" says only about the hypercube $\left[-\frac{1}{2},\frac{1}{2}\right]^d \text{, } d\in \mathbb{N},d>1$. Is your usage to an improper integral documented? Nov 1 '21 at 5:58
• Every method has its limitations: ListPlot[Thread[{Table[t02, {t02, 0, 160, 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, 0, t02, Infinity}, {t2, -Infinity, 0, t02, Infinity}, Method -> {"MultidimensionalRule", "Generators" -> 9}], {t02, 0, 160, 1}]}], PlotRange -> Full, Joined -> True, Frame -> True, FrameLabel -> {"t02-t01", "Overlap"}] ( 160 instead of 20) produces a bad plot Nov 1 '21 at 7:05
• @user64494 To follow the advice to "find its numeric support" when $t_{01}-t_{02}$ is too large for the standard recursive subdivision, usually one subdivides further at $\mu\pm a\sigma$, with $a$ around $5$ or $6$ for a precision goal of around 6. As you are a user who has complained at others' interest in ridiculously large or small numbers, I wonder why 160 is at all useful. f\[Tau]gaussian[1, 0, 160] is on the order of $10^{-5660}$. If you just want to stymie Mma, use 10^8. --- For the (centered) unit cube and infinite intervals, read the rest of the NIntegrate documentation. Nov 1 '21 at 12:58
• Thank you for your reply. You wrote " For the (centered) unit cube and infinite intervals, read the rest of the NIntegrate documentation". I repeat I don't find any example of the usage of "MultidimensionalRule" to improper integrals in the documentation to NIntegrate. Can you give a concrete reference? TIA. Nov 1 '21 at 18:13
• @user64494 This is too short to completely document the workings of NIntegrate, but if you read it carefully, you will discover how NIntegrate handles infinite intervals. You may have to experiment to fill in missing details, but such is the nature of the documentation that we have. The option IntegrationMonitor, described on this site, is a useful tool for investigating the operation of NIntegrate. Good luck with your researches! Nov 2 '21 at 18:55

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"}]