# 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 – j.foobles Nov 15 '19 at 0:07
• Is this ok \[Omega]0, \[Delta], Int1^2, Int2^2 , Ineq ={0.001,1.49957287164036,0.395907960255947,0.973870207974076}.? – Alex Trounev Nov 15 '19 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'. – j.foobles Nov 15 '19 at 22:58
• Then you incorrectly formulated the question. Do you want to study the convergence of the integral for $\delta \rightarrow Infinity$? – Alex Trounev Nov 16 '19 at 11:03
• Please add the Mathematica code to your question, not just LaTeX formulas. – Anton Antonov Dec 15 '19 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*\[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*\[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? – j.foobles Nov 15 '19 at 23:00
• To use option WorkingPrecision ->p with an arbitrary p, all numbers must be taken with the same accuracy, integers or rationals. – Alex Trounev Nov 16 '19 at 10:52