Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

migrated from stackoverflow.com Feb 11 '12 at 2:59

This question came from our site for professional and enthusiast programmers.

1  
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 –  Szabolcs Feb 10 '12 at 23:11
1  
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? –  halirutan Feb 11 '12 at 0:09
    
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. –  J. M. Feb 11 '12 at 3:43
    
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... –  Rojo Feb 11 '12 at 17:32
1  
I'd personally suggest an initial treatment of PiecewiseExpand[], followed by FullSimplify[]... –  J. M. Feb 13 '12 at 17:16

1 Answer 1

up vote 11 down vote accepted

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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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