# NIntegration and GlobalAdaptive method?

I integrate with different methods and obtain different results. How to understand which one is correct?

T := 20000;
I3s := 1./(2*3^2);
\[Sigma][w_]:=1/(2Pi) 1/(x3s[w]w) 2^9 Pi^3 3^2 (x3s[w])^7 (((x3s[w])^2+1)(7(x3s[w])^2+27)^2)/((x3s[w])^2+9)^6 E^(-4x3s[w]*ArcCot[x3s[w]/3])/(1-E^(-2\[Pi]*x3s[w]));
nB[w_]:=1/(Exp[w/(kB*T)]-1);
G3s := (1.8980259927018785*10^8)/(4.1341*10^16);
x3s[w_] := Sqrt[0.5/(w - I3s)];


and finally the integral

NIntegrate[w^2 \[Sigma][w]/w Exp[-(w/(kB*T))](w1*nB[w1](G3s/((w-w1)^2+G3s^2/4)+G3s/((w+w1)^2+G3s^2/4))),{w,I3s,Infinity},{w1,0,Infinity}]


• Please, post your code to understand what you are doing... – José Antonio Díaz Navas Mar 25 '18 at 11:50
• How to paste mathematica code here? – Timur Zalialiutdinov Mar 25 '18 at 11:54
• T := 20000 I3s := 1./(2*3^2) G3s := (1.8980259927018785*10^8)/(4.1341*10^16) – Timur Zalialiutdinov Mar 25 '18 at 12:53
• kB := 3.16681*10^-6 – Timur Zalialiutdinov Mar 25 '18 at 13:44
• Kindly edit your question to include all the additional information you gave in the comments, for completeness. – J. M.'s ennui Mar 25 '18 at 14:33

Have a look at the function, where the major contributions to the integral are. (Attention: This is a Log-Plot)

fs[w_, w1_] =
FullSimplify[
w^2 \[Sigma]3sw[w]/
w Exp[-(w/(kB*T))] (w1*
nB[w1] (G3s/((w - w1)^2 + G3s^2/4) +
G3s/((w + w1)^2 + G3s^2/4))), {w > I3s, w1 > 0}];

Plot3D[Log@fs[w, w1], {w, I3s, 2}, {w1, 0, 2}, PlotRange -> All,
PlotPoints -> 50, MaxRecursion -> 5]


The function peaks at w1 is close to w. (Spikes in the plot are due to low resolution)

Plot3D[Log@fs[w, w + r], {w, I3s, 2}, {r, - I3s, I3s},
PlotRange -> All, PlotPoints -> 50, MaxRecursion -> 5]


So you can limit the region of integration

NIntegrate[fs[w, w1], {w, I3s, 2}, {w1, w - I3s, w + I3s},
MaxRecursion -> 50]

(*   0.017838   *)


Integrate to a higher w, to get slightly better integral.

NIntegrate[fs[w, w1], {w, I3s, 4}, {w1, w - I3s, w + I3s},
MaxRecursion -> 50, WorkingPrecision -> 50, PrecisionGoal -> 10,
AccuracyGoal -> 10] // Quiet // N[#, 6] &

(*   0.0178421   *)

• Your last integral returns multiple errors of type: The integrand has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries ..., How do you get the numerical result then? What is your MMA version? BY using Method->"LocalAdaptive" I obtain for your integral the value 2.32069*10^-9, which is consistent with my result... – José Antonio Díaz Navas Mar 27 '18 at 12:23
• I am using MMA Version 8.0 . I rationalized all constants, so I can raise WorkingPrcision  as high as needed. Try higher WorkingPrecision. I had no problems with the code shown. I find it no good idea to use "LocalAdaptive" and integrate to infinty. It does not recognize the peakline you see in my pictures. Do re2 = Reap[NIntegrate[.., EvaluationMonitor :> Sow[{w, w1, Log@f[w, w1]}]]][[2, 1]]; with your integral to see, there are not more integration points at the peakline with : ListPointPlot3D[re2, PlotRange -> {{I3s, 1}, {0, 1}, {-80, 25}}]  but my integral resolves it exactly. – Akku14 Mar 27 '18 at 15:36
• No way, even doing what you say. I got the same errors. I have redone it by MonteCarlo and taking a mean of 50 runs, I got 2.26271*10^-6 for your integral. On the other hand, GlobalAdaptive is not indicated for peaked functions, more appropriate LocalAdaptive, so I cannot figure out why you find it "no good idea"... – José Antonio Díaz Navas Mar 27 '18 at 15:58
• related – José Antonio Díaz Navas Mar 27 '18 at 16:27
• @JoséAntonioDíazNavas Subdividing the interval, as I did in my answer to that related question, I get a result similar to Akku14's. – Michael E2 Mar 27 '18 at 20:11

Following @Anton Antonov's setup of rationalizing the parameters:

Clear[kB, T, I3s, σ, nB, G3s, x3s]
kB = 316681*10^-11;
T = 20000;
I3s = 1/(2*3^2);
σ[w_] := 1/(2 Pi) 1/(x3s[w] w) 2^9 Pi^3 3^2 (x3s[w])^7 (((x3s[w])^2 + 1) (7 (x3s[w])^2 + 27)^2)/((x3s[w])^2 + 9)^6 E^(-4 x3s[w]*ArcCot[x3s[w]/3])/(1 - E^(-2 π*x3s[w]));
nB[w_] := 1/(Exp[w/(kB*T)] - 1);
G3s = (18980259927018785*10^-8)/(41341*10^12);
x3s[w_] := Sqrt[(1/2)/(w - I3s)];


### Addressing the near-singularity

It's not hard to see from a term in the integrand that w1 == w has a delta-like, near-singularity:

G3s/((w - w1)^2 + G3s^2/4) // N
(*  4.59115*10^-9/(5.26966*10^-18 + (w - 1. w1)^2)  *)


If we add w to the iterator for w1, NIntegrate will subdivide the interval at w, which will address a substantial part of the problem.

NIntegrate[w^2 σ[w]/w *
Exp[-(w/(kB*T))] (w1*nB[w1] (G3s/((w - w1)^2 + G3s^2/4) + G3s/((w + w1)^2 + G3s^2/4))),
{w, I3s, Infinity}, {w1, 0, w, Infinity},
MaxRecursion -> 20,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None}
] // AbsoluteTiming


[No messages]

{0.130505, 0.0177845}


### Assessing convergence

To test the accuracy, standard tricks include increasing PrecisionGoal, WorkingPrecision, and sometimes MinRecursion. The default for PrecisionGoal on a 2D integral is 6. Let's try 8:

NIntegrate[w^2 σ[w]/w *
Exp[-(w/(kB*T))] (w1*nB[w1] (G3s/((w - w1)^2 + G3s^2/4) + G3s/((w + w1)^2 + G3s^2/4))),
{w, I3s, Infinity}, {w1, 0, w, Infinity},
MaxRecursion -> 20, PrecisionGoal -> 8,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None}
] // AbsoluteTiming


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

{0.936032, 0.0178419}


The message NIntegrate::slwcon is just a warning and does not indicate an error. In both cases, there was no message indicating that the error goals were not met. However, there is a large difference between the two results, a difference in the 3rd significant figure. This suggests the error estimate is not accurate and casts some doubt on the convergence of numerical integration. Considering the advice in the message, the integrand does not appear to be oscillating and so it could be that the WorkingPrecision is too small. I've also found in such cases that increasing MinRecursion can remove the warning. Increasing it causes denser sampling throughout the domain, which can help the error estimate. (It has the downside of increasing sampling in places where it's not needed. For each increment of 1, it usually doubles the sampling in each dimension, so the sampling grows exponentially with the level of recursion.)

Usually I do these things step-by-step by hand and think about the results at each step. But here's a way to increase MinRecursion until it does no more good. It compares one level with the previous level until the relative change is less than 10^-6. Hence it overshoots the maximum useful MinRecursion, and we have to step it back at the end.

ClearAll[nonconvQ];
nonconvQ[{___, r1_, r2_}] := Abs[(r1 - r2)/(10^-16 + 10^-6 r1)] > 1;
nonconvQ[{___}] := True;

PrintTemporary@Dynamic@{mr, Clock[Infinity]};  (* monitor *)
mr = 0;    (* initial value for MinRecursion *)
pg = 8;    (* PrecisionGoal *)
mrresults = {};
Print[PrecisionGoal -> pg];  (* to remind me/check that it's set as desired *)
While[nonconvQ@mrresults && mr <= 10,  (* upper limit of 10 *)
AppendTo[
mrresults,
Check[
res = NIntegrate[
w^2 σ[w]/w Exp[-(w/(kB*T))] *
(w1*nB[w1] (G3s/((w - w1)^2 + G3s^2/4) + G3s/((w + w1)^2 + G3s^2/4))),
{w, I3s, Infinity}, {w1, 0, w, Infinity},
MinRecursion -> mr, MaxRecursion -> 20, PrecisionGoal -> pg,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None}],
Print[MinRecursion -> mr]; res]
];
mr += 2;
] // AbsoluteTiming
mr = mr - 4

PrecisionGoal->8


NIntegrate::slwcon: ....

MinRecursion->0


NIntegrate::slwcon: ....

MinRecursion->2


NIntegrate::slwcon: ....

General::stop: ....

MinRecursion->4

{2.97894, Null}

2  (* <-- mr *)


We found that MinRecursion -> 2 is the best, although we still get NIntegrate::slwcon warnings. (What works may also depend on the PrecisionGoal and WorkingPrecision, but this gives us a good starting point.) We can examine the results and their differences:

mrresults

{0.017841865511943975,
0.01784207781124873,
0.01784207380554739} // Differences

(*  {2.12299*10^-7, -4.0057*10^-9}  *)


So let's try raising WorkingPrecision. I usually double machine precision as a first try, to 30 or 32. To see whether WorkingPrecision -> 30 is working, we can increasingly raise PrecisionGoal. Instead of checking for convergence as we did for MinRecursion, we'll examine the results stored in hpg (High PrecisionGoal) for stability. The differences of the results show that the last two results for PrecisionGoal of 10 and 12 agree to 9 significant digits; probably, therefore, a precision goal of at least 9 has been reached in the final result.

PrintTemporary@Dynamic@{pg2, Clock[Infinity]};  (* monitor *)
Print[MinRecursion -> mr];  (* check it's set as desired (mr determined in prev. code) *)
hpg = Table[
Check[
res = NIntegrate[
w^2 σ[w]/w Exp[-(w/(kB*T))] *
(w1*nB[w1] (G3s/((w - w1)^2 + G3s^2/4) + G3s/((w + w1)^2 + G3s^2/4))),
{w, I3s, Infinity}, {w1, 0, w, Infinity},
MinRecursion -> mr, MaxRecursion -> 20,
PrecisionGoal -> pg2, WorkingPrecision -> 30,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None}],
Print[PrecisionGoal -> pg2]; res],
{pg2, pg, 12, 2}  (* tries 8, 10, 12 *)
]; // AbsoluteTiming
hpg

MinRecursion->2

{279.558, Null}

{0.0178419717276098694849061712810,  (* these are the results in hpg *)
0.0178420778761599013231241349024,
0.0178420779116383260908787283014} // Differences

{1.061485500318382179636214*10^-7, 3.54784247677545933990*10^-11} (* Differences@hpg *)


This shows that @Akku14's result of 0.0178421 seems accurate to the digits shown.

• However, I still would like an explanation of why by using LocalAdaptive, or MonteCarlo, one obtain so different result...If I do the same as you, I obtain your result, but the OP question was related about the reliability of the result of NIntegrate by using different methods, i.e., what is the correct result? – José Antonio Díaz Navas Mar 28 '18 at 11:17
• @JoséAntonioDíazNavas I think the erroneous results in both cases are explained in this tutorial. I think "LocalAdaptive" does work, esp. if you use {w1, 0, w, Infinity}, but not as well as "GlobalAdaptive", as explained here. – Michael E2 Mar 28 '18 at 13:35
• I'm not sure why one would think "MonteCarlo" would work after seeing Akku's plot. When the support is more spread out, one expects convergence to be roughly $\sigma/\sqrt{n}$, where $\sigma$ is a bound on the standard deviation of the $n$ random function samples. I don't know how having such a skewed distribution affects the error, but I'd guess it makes it worse. It also makes it more likely that the sampling will missed the support along the line w == w1. – Michael E2 Mar 28 '18 at 13:35
• Note AccuracyGoal -> 5 would make it stop when the absolute error estimate is below 10^-5, which is likely to happen with a small sample size. You might need a sample size over 10^10, maybe 10^14 to get a smallish error. This works somewhat NIntegrate[integrand, {w, I3s, 0.31}, {w1, 0, w, 0.31}, Method -> {"MonteCarlo", "MaxPoints" -> 10^7}], although it varies wildly with each evaluation. – Michael E2 Mar 28 '18 at 13:35
• Thanks for your efforts. However, I know and read that tutorial, and there is no place whether it is clear that GlobalAdaptive is better than LocalAdaptive. Further "In general the global adaptive strategy has better performance than the local adaptive one. In some cases though, the local adaptive strategy is more robust and/or gives better performance.". Thus, why is not the analysed case one of those? Why does not it overestimate the integral by subdividing the intervals and doing GlobalAdaptive with GaussKronrodRule? – José Antonio Díaz Navas Mar 28 '18 at 17:05

As @JoséAntonioDíazNavas mentioned it seems that your integral is zero, or very close to zero.

Using these definitions with rationalized numbers:

Clear[kB, T, I3s, \[Sigma], nB, G3s, x3s]
kB = 316681*10^-11;
T = 20000;
I3s = 1/(2*3^2);
\[Sigma][w_] :=
1/(2 Pi) 1/(x3s[w] w) 2^9 Pi^3 3^2 (x3s[
w])^7 (((x3s[w])^2 + 1) (7 (x3s[w])^2 + 27)^2)/((x3s[w])^2 +
9)^6 E^(-4 x3s[w]*ArcCot[x3s[w]/3])/(1 - E^(-2 \[Pi]*x3s[w]));
nB[w_] := 1/(Exp[w/(kB*T)] - 1);
G3s = (18980259927018785*10^-8)/(41341*10^12);
x3s[w_] := Sqrt[(1/2)/(w - I3s)];


I get the following estimates.

NIntegrate[
w^2 \[Sigma][w]/
w Exp[-(w/(kB*T))] (w1*
nB[w1] (G3s/((w - w1)^2 + G3s^2/4) +
G3s/((w + w1)^2 + G3s^2/4))), {w, I3s, Infinity}, {w1, 0,
Infinity},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 10000,
"SingularityHandler" -> None}, MinRecursion -> 4,
MaxRecursion -> 20, WorkingPrecision -> 30, PrecisionGoal -> 4]

(* NIntegrate::slwcon ...*)
(* NIntegrate::eincr...
NIntegrate obtained 0.000036654314889954921292480476695230. and \
6.1891071264343163397372600942630.*^-6 for the integral and error \
estimates.
*)
(* 0.0000366543148899549212924804766952 *)

NIntegrate[
w^2 \[Sigma][w]/
w Exp[-(w/(kB*T))] (w1*
nB[w1] (G3s/((w - w1)^2 + G3s^2/4) +
G3s/((w + w1)^2 + G3s^2/4))), {w, I3s, Infinity}, {w1, 0,
Infinity},
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule",
"MaxErrorIncreases" -> 10000}, MinRecursion -> 4,
MaxRecursion -> 20]

(* NIntegrate::slwcon...
NIntegrate::pincer...
NIntegrate obtained 0.0002643264047745091 and    0.0000943240383048598 \
for the integral and error estimates.*)
(* 0.0002643264047745091 *)


You can experiment with using higher precision and "MaxErrorIncreases" and see do you consistently get closer to 0. (Or any other value...)

By using your definitions in the notebook provided, I have done:

NIntegrate[
w^2 \[Sigma]3sw[w]/
w^1 Exp[-(w/(kB*T))] (w1^1*
nB[w1] (G3s/((w - w1)^2 + G3s^2/4) +
G3s/((w + w1)^2 + G3s^2/4))), {w, I3s, Infinity}, {w1, 0, Infinity},
MinRecursion -> 40, MaxRecursion -> 50,
Method -> {"LocalAdaptive", Method -> "MultidimensionalRule"},
PrecisionGoal -> \[Infinity], AccuracyGoal -> 5]

(* 1.50775*10^-8 *)
`

The result was obtained with no errors. I used MMA 11.3 and macOS 10.13.13