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.
