# How do I get my simulation to run?

I am simulating a birth-death population process. My code is

RESULTS = {};
For[iterator = 1, iterator < 21, iterator++,
n = 25; t = 0;
results = {{t, n}};
mu = 1; lam = 1;
While[t < 200,
death = lam*n;
birth = mu*n;
rate = death + birth;
deltaT[r_] := -1/r*Log[RandomReal[]];
t1 = deltaT[rate];
t = t + t1;
rand = RandomReal[]*rate;
Which[
rand <= death, n = n - 1,
True, n = n + 1
];
results = Append[results, {t, n}];
];
AppendTo[RESULTS, results];
]


I am getting an error message that says 'Infinite expression 1/0 encountered. >>'.

I'm not quite sure why this is. I think it might be because the population dies out so the rate would be zero, hence deltaT[r_] would be trying to work out -1/0.

Is this correct?

How do I get around this?

• I understand that's the reason why. I therefore want the simulation to run up until the point to which the rate becomes 0, then stop. How do I adapt my code so that it does this? @HenrikSchumacher – Mlo27 Apr 14 '18 at 15:49
• @Mlo27 For rate to be zero, you have to have n = 0, right? So you need to change your Which statement so that n does not become zero. – C. E. Apr 14 '18 at 15:50
• You could put a If[Abs[rate]<100. $MachineEpsilon,Break[]] before the point of division. – Henrik Schumacher Apr 14 '18 at 16:01 • Have you seen this? mathematica.stackexchange.com/questions/158212/… – OkkesDulgerci Apr 14 '18 at 17:32 • @OkkesDulgerci I have. I don’t want to change code massively though. I just want to add a command so that the simulation stops when I get the error. – Mlo27 Apr 14 '18 at 17:49 ## 2 Answers As @Henrik mentioned in comment, you can break simulation if rate<0. RESULTS = {}; For[iterator = 1, iterator < 21, iterator++, n = 25; t = 0; results = {{t, n}}; mu = 1; lam = 1; While[t < 200, death = lam*n; birth = mu*n; rate = death + birth; If[Abs[rate] < 100.$MachineEpsilon, Break[]];
deltaT[r_] := -1/r*Log[RandomReal[]];
t1 = deltaT[rate];
t = t + t1;
rand = RandomReal[]*rate;
Which[rand <= death, n = n - 1, True, n = n + 1];
results = Append[results, {t, n}];];
AppendTo[RESULTS, results];]

ListLinePlot[RESULTS, InterpolationOrder -> 0] I have tried to make a functional version of the original code that incorporates the observation made in the comments by multiple discussants that n=0 breaks the original code.

This implementation uses NestWhileList to substitute for the innermost While and NestList replaces the outermost For.

Block[{test, compose},

With[{no = 25, to = 0, te = 200, mu = 1, lam = 1, iterator = 20, seedo = 123456987},

(* input {t, n, r}; returns 'True' is t < te (=200) and r \[NotEqual] 0 *)
test = #1 < te && #3 != 0 &;

(* input {t, n}; returns {t, n, r} *)
compose = {##, (lam + mu) #2} &;

(* get rid of the first entry in the output (= NestList's second argument (check documentation) *)
data = Rest@NestList[

(* create the seed to be used in BlockRandom (used for reproducibility) *)
With[{seed = #[] + 1},

(* output = {{seed, {{t, n, r},.. }}.. } *)
{seed, BlockRandom[

NestWhileList[

With[{t = #[], n = #[], r = #[[-1]], rands = RandomReal[{0, 1}, 2]},

(* perhaps redundant With - helps readability *)
With[{death = lam n},

(* produce next {t, n, r} from {t, n} *)
Apply[compose][{
t - 1/r Log[rands[]],
Which[
r rands[[-1]] <= death, n - 1,
True, n + 1
]
}]

]
] &,

(* init = {t, n, r} *)
Apply[compose][{to, no}],

(* check if t < 200 and r\[NotEqual]0 *)
Apply[test][#] &], RandomSeeding -> #[]]}

(* init NestList *)
] &, {seedo, {}}, iterator];

]

] The plots were produced with:

ListLinePlot[
Part[#, -1, All, {1, 2}],
Frame -> True,
PlotLegends -> Placed[Row[{HoldForm[seed] -> Part[#, 1]}], Bottom],
InterpolationOrder -> 0,
PlotRange -> {{0, 200}, Automatic}
] & /@ data // Partition[#, 4, 4, {1, 1}, {}] & // Grid