# 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. Commented Oct 31, 2021 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? Commented Nov 1, 2021 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 Commented Nov 1, 2021 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. Commented Nov 1, 2021 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. Commented Nov 1, 2021 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! Commented Nov 2, 2021 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"}]