0
$\begingroup$

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

I then tried 'AdaptiveMonteCarlo': 

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)

Improve this question
$\endgroup$
5
  • $\begingroup$ This is the link to the Mathematica script I am using: 1drv.ms/u/s!Ak_7kWMZ0x15g80EhS2Xnk852B9_rg?e=ucx3m6 $\endgroup$
    – j.foobles
    Commented Nov 15, 2019 at 0:07
  • $\begingroup$ Is this ok \[Omega]0, \[Delta], Int1^2, Int2^2 , Ineq ={0.001,1.49957287164036,0.395907960255947,0.973870207974076}.? $\endgroup$
    – Alex Trounev
    Commented Nov 15, 2019 at 12:05
  • $\begingroup$ 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'. $\endgroup$
    – j.foobles
    Commented Nov 15, 2019 at 22:58
  • $\begingroup$ Then you incorrectly formulated the question. Do you want to study the convergence of the integral for $\delta \rightarrow Infinity$? $\endgroup$
    – Alex Trounev
    Commented Nov 16, 2019 at 11:03
  • $\begingroup$ Please add the Mathematica code to your question, not just LaTeX formulas. $\endgroup$
    – Anton Antonov
    Commented Dec 15, 2019 at 15:21

1 Answer 1

Reset to default
2
$\begingroup$

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}.*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks Alex, sorry for my naivety but what is the importance of writing parameters in rational form? $\endgroup$
    – j.foobles
    Commented Nov 15, 2019 at 23:00
  • 1
    $\begingroup$ To use option WorkingPrecision ->p with an arbitrary p, all numbers must be taken with the same accuracy, integers or rationals. $\endgroup$
    – Alex Trounev
    Commented Nov 16, 2019 at 10:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.