# NIntegrate a Triple Integral... Different Integration Strategies

I'm trying to evaluate a triple integral using NIntegrate. It takes the form:

\begin{align} \int_{\omega_0-5\delta}^{\omega_0+5\delta} \mathrm{d}\omega \int_{\omega_0-5\delta}^{\omega_0+5\delta}\mathrm{d}\omega' \int_{-1+\epsilon}^{1+\epsilon}\mathrm{d}s\: g(\omega)g(\omega')F(s) \end{align} where \begin{align} F(s) \propto \frac{1}{(s+2n+1)(s+2n-1)}\left(\frac{1+s}{1-s}\right)^{i \omega}\left( \frac{s+2n+1}{s+2n-1} \right)^{-i\omega'} \end{align} The $$\epsilon$$ is there as a cut-off for divergences near $$\pm1$$. Anyway, I've tried integrating this with the automatic NIntegrate with the following settings, which take a long time:

(e.g. Method -> {"AdaptiveMonteCarlo", "SymbolicProcessing" -> 0} and AccuracyGoal -> 6, PrecisionGoal -> 6)

and don't end up being accurate anyway (I have a physical intuition as to the maximum value of such integrals) (even though I don't get any warnings).

Method -> {"AdaptiveMonteCarlo",
"SymbolicProcessing" ->
0}, AccuracyGoal -> MinRecursion -> 10^8, MaxRecursion -> Infinity


which is faster, but gives convergence errors:

NIntegrate obtained 1.41829 -0.141122\ I and 0.005438535119730085 \
for the integral and error estimates."


and still yields a slightly wrong answer. Does anyone have any suggestions about the integration strategy I could use in this case, or options on the MonteCarlo integration that would help?

(I have attached the example script here since it didn't format correctly by copy-pasting. The quantity Ineq should be less than 1 but as shown in the .nb, it is 1.10927)

• This is the link to the Mathematica script I am using: 1drv.ms/u/s!Ak_7kWMZ0x15g80EhS2Xnk852B9_rg?e=ucx3m6 Commented Nov 15, 2019 at 0:07
• Is this ok \[Omega]0, \[Delta], Int1^2, Int2^2 , Ineq ={0.001,1.49957287164036,0.395907960255947,0.973870207974076}.? Commented Nov 15, 2019 at 12:05
• Hi Alex, that seems to be an okay answer, however strangely, the integral doesn't seem to converge (to a single value) when increasing the upper limit of integration on the w, w' integrals (which is expected since g(w) have exponentially decaying tails). I've set the upper limits to infinity and gone back to 'LocalAdaptive'. Commented Nov 15, 2019 at 22:58
• Then you incorrectly formulated the question. Do you want to study the convergence of the integral for $\delta \rightarrow Infinity$? Commented Nov 16, 2019 at 11:03
• Please add the Mathematica code to your question, not just LaTeX formulas. Commented Dec 15, 2019 at 15:21

It is necessary to take all parameters in a rational form and change the integration strategy, then the integrals converge

ClearAll["Global*"]
\[Omega]0 = 1/100;
\[Delta] = 1/10;
a = 1;
\[Epsilon] = $$MachineEpsilon; n = 10^(-3); w0M =$$MachineEpsilon;
w0P = \[Omega]0 + 5*\[Delta];
Acc1 = 10;
Acc2 = 5;
MP = 10000000;
gw = (Sqrt[\[Omega]]*Exp[-((\[Omega] - \[Omega]0)^2/(4*\[Delta]^2))])/(Exp[-(\[Omega]0^2/(2*\[Delta]^2))]*\[Delta]^2 + Sqrt[Pi/2]*\[Omega]0*\[Delta]*(1 + Erf[\[Omega]0/(Sqrt[2]*\[Delta])]))^2^(-1);
gwd = (Sqrt[\[Omega]d]*Exp[-((\[Omega]d - \[Omega]0)^2/(4*\[Delta]^2))])/(Exp[-(\[Omega]0^2/(2*\[Delta]^2))]*\[Delta]^2 + Sqrt[Pi/2]*\[Omega]0*\[Delta]*(1 + Erf[\[Omega]0/(Sqrt[2]*\[Delta])]))^2^(-1);
Aww = (1/(Pi*a))*Sqrt[\[Omega]d/(a*(\[Omega]/a))]*(1/((s + 2*n - 1)*(s + 2*n + 1)))*((1 + s)/(1 - s))^(I*\[Omega])*((s + 2*n + 1)/(s + 2*n - 1))^(I*\[Omega]d);
Bww1 = ((-(Pi*a)^(-1))*Sqrt[\[Omega]d/(a*(\[Omega]/a))]*(1/((s + 2*n - 1)*(s + 2*n + 1)))*((1 + s)/(1 - s))^(I*\[Omega]))/
((s + 2*n + 1)/(s + 2*n - 1))^(I*\[Omega]d);
Bww2 = ((-(Pi*a)^(-1))*Sqrt[\[Omega]/(a*(\[Omega]d/a))]*(1/((s + 2*n - 1)*(s + 2*n + 1)))*((s + 2*n + 1)/(s + 2*n - 1))^(I*\[Omega]))/
((1 + s)/(1 - s))^(I*\[Omega]d);
B10 = NIntegrate[gw*gwd*(Bww1/(2*Sinh[Pi*(\[Omega]d/a)])), {\[Omega], w0M, w0P}, {\[Omega]d, w0M, w0P}, {s, -1 + \[Epsilon], 1 - \[Epsilon]},
Method -> {"AdaptiveQuasiMonteCarlo", "SymbolicProcessing" -> 0}, WorkingPrecision -> 15, MaxRecursion -> 20];
B01 = NIntegrate[gw*gwd*(Bww2/(2*Sinh[Pi*(\[Omega]/a)])), {\[Omega], w0M, w0P}, {\[Omega]d, w0M, w0P}, {s, -1 + \[Epsilon], 1 - \[Epsilon]},
Method -> {"AdaptiveQuasiMonteCarlo", "SymbolicProcessing" -> 0}, WorkingPrecision -> 15, MaxRecursion -> 20];
Ineq = 1/2 + (1/4)*Abs[B10]^2 + (1/4)*Abs[B01]^2;
Print[StringForm["\[Omega]0, \[Delta], \!$$\*SuperscriptBox[\(Int1$$, $$2$$]\), \!$$\*SuperscriptBox[\(Int2$$, $$\(2$$$$\\$$\)]\), Ineq =.",
{N[n], Abs[B10]^2, Abs[B01]^2, Ineq}]];
(*\[Omega]0, \[Delta], Int1^2, Int2^2 , Ineq ={0.001,1.49957287164036,0.395907960255947,0.973870207974076}.*)

• Thanks Alex, sorry for my naivety but what is the importance of writing parameters in rational form? Commented Nov 15, 2019 at 23:00
• To use option WorkingPrecision ->p with an arbitrary p, all numbers must be taken with the same accuracy, integers or rationals. Commented Nov 16, 2019 at 10:52