# Optimizing Monte Carlo simulation of a Pred-Prey model

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 == H0, P == 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 == H0, P == 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?

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 == H0, P == 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}

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[]];
]

{39.557263, Null}

• 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 '14 at 17:19