16
$\begingroup$

I have an ugly, six dimensional function that I need to integrate numerically. It works, but it currently take twelve hours to complete the calculation. Is there any good way to parallelize the calculation to run across multiple cores?

The best I've managed to come up with is:

Sum[
    ParallelTable[
        NIntegrate[calc[a,b,c,d,e,f],
                   {a,-1,1},
                   {b,-1,1},
                   {c,-1,1},
                   {d,-1,1},
                   {e,-1,1},
                   {f,-1+i/4,-1+(i+1)/4}],
    {i,0,7}]]

EDIT: The value for calc is given below:

calc[a,b,c,d,e,f] = (1.97531*10^15 (3. Q Cos[(3 Q)/20000] - 
     20000. Sin[(3 Q)/20000])^2)/(1.58025*10^24 + 6.32099*10^16 Q^2 + 
   Q^6 + (-1.58025*10^24 + 7.90123*10^15 Q^2 - 3.55556*10^8 Q^4) Cos[(
     3 Q)/10000] + (-4.74074*10^20 Q - 2.96296*10^12 Q^3) Sin[(3 Q)/
     10000])*
    UnitStep[5/2 - a]^2 UnitStep[5/2 + a]^2 UnitStep[5 - b]^2 UnitStep[
      5 + b]^2 UnitStep[5 - a - 1944 c] UnitStep[5 + a + 1944 c] UnitStep[
      5 - b - 1944 d] UnitStep[5 + b + 1944 d] UnitStep[
      40 + a - 1600 e] UnitStep[40 - a + 1600 e] UnitStep[
      45/2 + b - 1600 f] UnitStep[45/2 - b + 1600 f] /. Q-> (8000000 \[Pi] Sin[
          1/2 ArcCos[
            1/2 Sqrt[Cos[2 c] + Cos[2 d]] Sqrt[Cos[2 e] + Cos[2 f]] + 
             Sin[c] Sin[e] + Sin[d] Sin[f]]])

There are two immediate things to note. First, there's a removable singularity at (c==e && d==f). I've also tried using a piecewise function to plug this discontinuity, but it doesn't seem to have a significant effect on the speed.

One other thought that has come to mind is using the UnitSteps to directly find the range of integration. I hadn't done that before mostly out of laziness and the understanding that Mathematica automatically breaks up piecewise functions for integration, anyway.

|improve this question
$\endgroup$
  • 1
    $\begingroup$ That is a decent way to parallelize it, but I'd suggest also looking at various integration methods to speed it up---some are faster than others for multidimensional. Take a look at Monte Carlo methods in particular $\endgroup$ – Szabolcs Feb 10 '12 at 23:11
  • 1
    $\begingroup$ The more important question is, how fast is calc? Is it a complex function of simple expressions or does it use sophisticated Mathematica methods too? $\endgroup$ – halirutan♦ Feb 11 '12 at 0:09
  • $\begingroup$ You will also want to exploit any structure in your integral (e.g. symmetry), and be mindful of variable substitutions that might make your problem slightly more tractable. $\endgroup$ – J. M. will be back soon♦ Feb 11 '12 at 3:43
  • $\begingroup$ I don't think that that paralellization will speed it up. That would imply that the time taken to integrate is proportional to the size of the region. Perhaps if you fixed the target accuracy of the parts... $\endgroup$ – Rojo Feb 11 '12 at 17:32
  • 1
    $\begingroup$ I'd personally suggest an initial treatment of PiecewiseExpand[], followed by FullSimplify[]... $\endgroup$ – J. M. will be back soon♦ Feb 13 '12 at 17:16
12
$\begingroup$

The best you can do it is to speed up your function. Your calc is using a replace, but it's better if you use With:

calc[a_, b_, c_, d_, e_, f_] := 
 With[{Q = (8000000 \[Pi] Sin[
       1/2 ArcCos[
         1/2 Sqrt[Cos[2 c] + Cos[2 d]] Sqrt[Cos[2 e] + Cos[2 f]] + 
          Sin[c] Sin[e] + 
          Sin[d] Sin[f]]])}, (1.97531*10^15 (3. Q Cos[(3 Q)/20000] - 
         20000. Sin[(3 Q)/20000])^2)/(1.58025*10^24 + 
      6.32099*10^16 Q^2 + 
      Q^6 + (-1.58025*10^24 + 7.90123*10^15 Q^2 - 
         3.55556*10^8 Q^4) Cos[(3 Q)/10000] + (-4.74074*10^20 Q - 
         2.96296*10^12 Q^3) Sin[(3 Q)/10000])*
   UnitStep[5/2 - a]^2 UnitStep[5/2 + a]^2 UnitStep[5 - b]^2 UnitStep[
     5 + b]^2 UnitStep[5 - a - 1944 c] UnitStep[
    5 + a + 1944 c] UnitStep[5 - b - 1944 d] UnitStep[
    5 + b + 1944 d] UnitStep[40 + a - 1600 e] UnitStep[
    40 - a + 1600 e] UnitStep[45/2 + b - 1600 f] UnitStep[
    45/2 - b + 1600 f]]

For working with parallel, it's better if you distribute your definitions before calling parallel:

DistributeDefinitions[calc]

Then try:

Total[ParallelTable[
   NIntegrate[
    calc[a, b, c, d, e, f], {a, -1, 1}, {b, -1, 1}, {c, -1, 
     1}, {d, -1, 1}, {e, -1, 1}, {f, -1 + i/4, -1 + (i + 1)/4}], {i, 
    0, 7}]] // AbsoluteTiming
|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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