# How to use parallel computation in mathematica for numerical Integration

I have a multi dimensional integral and I need to integrate it numerically and than plot it over an extra variable for zeta = 0.5.

g[kx_, ky_, k1x_, k1y_, qx_, qy_, zeta_, zz_] := (1/((1 + Exp[-Sqrt[k1x^2 + k1y^2] + zeta*k1x*Cos[zz] +
zeta*k1y*Sin[zz]])*(1 +
Exp[-Sqrt[kx^2 + ky^2] + zeta*kx*Cos[zz] +
zeta*ky*Sin[zz]])*(1 +
Exp[Sqrt[(kx + qx)^2 + (ky + qy)^2] + zeta*(kx + qx)*Cos[zz] +
zeta*(ky + qy)*Sin[zz]]*(1 +
Exp[Sqrt[(k1x - qx)^2 + (k1y - qy)^2] +
zeta*(k1x - qx)*Cos[zz] + zeta*(k1y - qy)*Sin[zz]]))));
PolarPlot[NIntegrate[g[kx, ky, k1x, k1y, qx, qy, 0.5, zz], {kx, -1, 1}, {ky, -1,
1}, {k1x, -1, 1}, {k1y, -1, 1}, {qx, -1, 1}, {qy, -1, 1}], {zz,
0, 2 Pi}]


my problem is that, it take many hours (more than 10 hours) without answer. my solution was to calculate the integral for some discrete set of angles, and than plot it, but I am searching for more better way.
based on searching in site,

Parallelizing Numerical Integration in Mathematica

it may be possible to use parallel computation to plot it. my knowledge in Mathematica is limited and I couldn't understand what the role of index "i" is in that link. Would someone please tell me how can I use all cores to decrease the time of calculation.

thanks

• The parameter i is the ParallelTable index in that code.It's computing 7 integral values and by the use of ParallelTable, these seven are being distributed (automatically) to the available cores. So if I wanted to create a table of 100 say with 4 cores, using ParallelTable would (approximately) run a total of 25 on each core and (approximately) reduce the overall execution time to 1/4 of time using the non-parallel Table command (if all calculation instances were approximately the same) – Dominic Apr 8 at 11:25
• @Dominic would you please write the code for parallelizing, I copy that code for the above integral, but it seems that index "i" is outside of the inner bracket of integral and doesn't seems to be seventh integral interval. – Arian Apr 8 at 12:31

## 1 Answer

Some comments:

(1) When you pass a function to a plotting routine, in general the routine will use a significant number of points for the plot, perhaps many more than you want so I just set up a table of 12 zz points below,

(2) Next, your integration as it stands is a bit troublesome so I included "Method->LocalAdaptive" which significantly increases the amount of time just to integrate one zz point but produces no errors as in the case of not doing so,

(3) When you supply even a single machine-precision number to Mathematica such as 0.5, every subsequent calculation involving that number will be in machine precision so instead I specified zeta=1/2 below which is arbitrary precision as well as a table of arbitrary precision values of zz. That way you can if you wish, integrate subsequently with arbitrary precision arithmetic.

Next, even one integration takes a bit of time:

AbsoluteTiming[
NIntegrate[
g[kx, ky, k1x, k1y, qx, qy, 1/2, 3/2], {kx, -1, 1}, {ky, -1,
1}, {k1x, -1, 1}, {k1y, -1, 1}, {qx, -1, 1}, {qy, -1, 1},
Method -> "LocalAdaptive"];
]
{109.281, Null}


So about 109 seconds but tolerable I think so I set up just 12 points running in parallel:

(*
set up shared kernel index to keep track of table index
*)
SetSharedVariable[theIndex]
Dynamic@theIndex
theIndex = 0;
(*
create list of 12 zz data points to integrate over in [0,2 Pi]
*)
myZData = Array[# &, 12, {0, 2 Pi}];
(*
print status
*)
PrintTemporary["Processing: ", Dynamic@theIndex];
(*
Pass a (parallel) table to the main kernel to distribute to other \
cores
*)
timing = AbsoluteTiming[
results = ParallelTable[
theIndex = myZ;
theVal =
NIntegrate[
g[kx, ky, k1x, k1y, qx, qy, 1/2, myZData[[myZ]]], {kx, -1,
1}, {ky, -1, 1}, {k1x, -1, 1}, {k1y, -1, 1}, {qx, -1,
1}, {qy, -1, 1}, Method -> "LocalAdaptive"];
{myZData[[myZ]], theVal},
{myZ, 1, Length@myZData}
];
];
timing[]/60.

7.41614


So about 7 1/2 minutes running on a 4.5 GHz quad-core machine running in parallel.

ListPolarPlot[results] • thank you so much, I expected that plot to be a circle, but I used "MonteCarloRule" method which had some error and were different for different runs. – Arian Apr 8 at 18:23