# Parallelization of Monte Carlo numerical integration

I'm trying to parallelize some integral. What I do: I define integrand:

Integrand105[p0_, p1_, t1_] =  2 (2*Pi)^4 (2/Pi)^4*(p0*
p1 (ArcTan[p0/2]/(4*Pi)) (ArcTan[
p1/2]/(4*Pi))/((p0^2 - 2*p0*p1*t1 + p1^2 + 1)^4*(p0^2 +
1)*(p1^2 + 1)));


Give a definition to all the kernels:

LaunchKernels[];
DistributeDefinitions[Integrand105];


The integral itself:

res = Total[
ParallelTable[
NIntegrate[
Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0,
Infinity}, {t1, -1 + i, i}, Method -> "MonteCarlo",
MaxPoints -> 1000000, PrecisionGoal -> 5], {i, 0, 1}]] //
AbsoluteTiming


That means that all parallelization consists in splitting the domain of integration. Ok. What I get: - from first kernel:

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 0.0024072378420660233 and 5.0413357069178085*^-6 for the integral and error estimates.


-from second one:

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 0.05468151281718856 and 0.0007526057830408755 for the integral and error estimates.


And in the end:

{3.569,0.05708}


So, my first question is: How can I get common estimation for error? Of course, I can sum these two error bars, but it will be the highest border. Any ideas? But even with this way, unfortunately, I have troubles:

I found hear: NIntegrate error bound how I can obtain errors:

Unprotect[Message];
Message[NIntegrate::maxp, its_, int_, err_] := Sow[err]


If I apply it just to integral:

NIntegrate[
Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0,
Infinity}, {t1, -1, 1}, Method -> "MonteCarlo",
MaxPoints -> 2000000, PrecisionGoal -> 10] // Reap


I get:

{0.0583558, {{0.0007278}}}


Excellent, I thought! I did all kinds of things to apply it for my case:

Total[ParallelTable[
NIntegrate[
Integrand105[p0, p1, t1], {p0, 0, Infinity}, {p1, 0,
Infinity}, {t1, -1 + i, i}, Method -> "MonteCarlo",
MaxPoints -> 1000000, PrecisionGoal -> 5] // Reap, {i, 0,
1}]] // AbsoluteTiming


And... {5.217, {0.05672, {}}}

Ok, may be without Total - no - {5.69, {{0.002407, {}}, {0.05574, {}}}}.

Also I tried other ways, but I didn't find smth that would help me.

I believe in the wisdom of the community. Thank you very much in advance.

p.s. In common, I'm looking for a way to calculate a integral(not like in the example above, integral as in the example one can calculate by means of, for instance:

Method -> {"GlobalAdaptive", Method -> "CartesianRule"},PrecisionGoal -> 15


or with:

"SingularityHandler" -> "DuffyCoordinates"}, PrecisionGoal -> 15


Within 5 minutes on a laptop. But applying this method, say, to 7-dimensional integral you 100% fail (If I'm wrong, please, make me the happiest man and correct me). I have an access to a cluster and I try to find an approach to calculate this(7-dim) integral with highest as far as possible precision.

• I would think that part of the problem is that the denominator of your integrand is zero for t1 == (1 + p0^2 + p1^2)/(2 p0 p1). Oct 29, 2017 at 4:44
• Rather than splitting the domain, could you just let the different parallel executions all try to evaluate the same integral and then treat them as independent measurements whose errors you know how to sum because you know the sample size of each? Oct 29, 2017 at 8:40
• Using a precision goal larger than 2 with a Monte Carlo method is not a good idea. Oct 29, 2017 at 9:19