Optimizing Monte Carlo simulation of a Pred-Prey model

Question

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 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!

Answer 1

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}}]]

---

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}