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.

My assignment and code

As part of an assignment for one of my classes, I'm trying to run a "massive" Monte Carlo simulation in Parallel on the follow model:

{H'[t] == r (1 - H[t]/k) - d H[t] P[t], 
 P'[t] == -s P[t] + e H[t] P[t], 
 H[0] == H0, P[0] == P0}

I've been given values for d, s and e but am required to use distributions for r, k, H0 and P0. I've implemented this as follows:

If[Kernels[] == {}, LaunchKernels[]];

kdis = RandomVariate[NormalDistribution[150, 20], 10]; 
rdis = RandomVariate[NormalDistribution[0.4, 0.003], 10];
H0dis = RandomVariate[UniformDistribution[{50, 200}], 10];
P0dis = RandomVariate[UniformDistribution[{20, 70}], 10];
d = 0.01;
s = 0.3;
e = 0.02;
SetSharedVariable[d, s, e];

distr = Tuples[{kdis, rdis, H0dis, P0dis}];

Eqs[{k_, r_, H0_, P0_}] := 
 NDSolve[{H'[t] == r (1 - H[t]/k) - d H[t] P[t], 
          P'[t] == -s P[t] + e H[t] P[t], 
          H[0] == H0, P[0] == P0}, {H, P}, {t, 0, 600}]
DistributeDefinitions[Eqs];

AbsoluteTiming[solset = ParallelMap[Eqs, distr];]

To display the results I create an envelope enclosing all runs for both predator and prey populations:

Envelope[n_] =
  Map[{Min@#, Max@#} &, 
   Transpose[Map[Flatten[{H[n], P[n]} /. #] &, solset]]];
DistributeDefinitions[Envelope];

AbsoluteTiming[plotset = ParallelMap[Envelope, Range[0, 600, 1]];]

plotsetH = plotset[[All, 1]];
plotsetP = plotset[[All, 2]];
Show[ListPlot[{plotsetH[[All, 1]], plotsetH[[All, 2]]}, 
               Joined -> True, PlotRange -> {{0, 600}, {0, 25}}],
     ListPlot[{plotsetP[[All, 1]], plotsetP[[All, 2]]}, 
               Joined -> True, PlotRange -> {{0, 600}, {0, 25}}]]

The result looks decent:

The result

The problem

Now my issue is that I'm currently hardly using proper distributions, but only 10 values per variable, which I doubt is what my instructor meant when he asked for a massive Monte Carlo simulation. The notebook takes about 2 minutes to evaluate and eats up almost all of the memory 8gb memory I have available in that time. Using 11 instead of 10 values per variable results in an error about kernels appearing dead, which I assume is caused by my memory running out. Are there any ways to optimize this kind of code, or to decrease memory usage?

Thanks in advance!

share|improve this question

1 Answer 1

up vote 7 down vote accepted

With

eqn[{k_, r_, H0_, P0_}] := {H'[t] == r (1 - H[t]/k) - d H[t] P[t], 
 P'[t] == -s P[t] + e H[t] P[t], H[0] == H0, P[0] == P0}
d = 0.01;
s = 0.3;
e = 0.02;

I would define one simulation as

sim := Module[
 {k = RandomVariate[NormalDistribution[150, 20]], 
  r = RandomVariate[NormalDistribution[0.4, 0.003]],
  H0 = RandomVariate[UniformDistribution[{50, 200}]], 
  P0 = RandomVariate[UniformDistribution[{20, 70}]]},
 Flatten@NDSolve[eqn[{k, r, H0, P0}], {H, P}, {t, 0, 600}]
 ]

and then perform the Monte Carlo simulation using

AbsoluteTiming[solset = ParallelTable[sim, {10000}];]

{13.021745, Null}

For the envelops parallelization gives no speed advantage, therefore

envelope[n_] = {Min@#, Max@#} & /@ (Transpose[{H[n], P[n]} /. solset]);
AbsoluteTiming[plotset = envelope[#] & /@ Range[0, 600, 1];]

{97.699588, Null}

and then your plot code

plotsetH = plotset[[All, 1]];
plotsetP = plotset[[All, 2]];
Show[ListPlot[{plotsetH[[All, 1]], plotsetH[[All, 2]]}, 
  Joined -> True, PlotRange -> {{0, 600}, {0, 25}}], 
 ListPlot[{plotsetP[[All, 1]], plotsetP[[All, 2]]}, Joined -> True, 
  PlotRange -> {{0, 600}, {0, 25}}]]

plot


Even Faster

If you don't need the InterpolatingFunctions for other calculations, you can get the plots even faster with

sim := Module[
 {k = RandomVariate[NormalDistribution[150, 20]], 
  r = RandomVariate[NormalDistribution[0.4, 0.003]],
  H0 = RandomVariate[UniformDistribution[{50, 200}]], 
  P0 = RandomVariate[UniformDistribution[{20, 70}]],
  sol},
 sol = Flatten@NDSolve[eqn[{k, r, H0, P0}], {H, P}, {t, 0, 600}];
 {H[#], P[#]} /. sol & /@ Range[0, 600, 1]
 ]

and then

AbsoluteTiming[
 solset = ParallelTable[sim, {10000}];
 plotset = Function[n, {Min@#, Max@#} & /@ Transpose[solset[[All, n]]]] /@ 
  Range[Length@solset[[1]]];
]
{39.557263, Null}      
share|improve this answer
    
Thank you! I wouldn't have thought of defining it like that, it certainly looks quite a bit cleaner as well :) –  Jason Houtekamer Aug 30 at 17:19

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.