In [1]:
start_time = 0;
end_time = 100;

# How much time passes between each successive calculation
time_step = 1/4 # In years
end_step = int(((end_time - start_time)/ time_step));

# The number of zebras when the simulation is started
initial_zebras = 30000;

# The number of lions when the simulation is started
initial_lions = 15;

# The number of zebras that must be killed for a lion to be born. (Zebras / Lion)

zebras_per_lion = 1000;

# The chance a zebra will die when a zebra lion cross paths
fighting_chance = 0.50;

meeting_chance = 0.02;

zebra_growth_rate = 0.20;

lion_death_rate = 0.10;

# Initialization

zebras_over_time = fill(0.0, end_step+1);
lions_over_time = fill(0.0, end_step+1);
model_time = fill(0.0, end_step+1);

zebras = initial_zebras;
lions = initial_lions;

zebras_over_time[1]	= zebras;
lions_over_time[1] = lions;

In [2]:
# Run the model

for sim_step = 1:end_step

	sim_time = start_time + sim_step * time_step;
	model_time[sim_step] = sim_time;

	# First we must calculate our flows (our rates)
	zebra_births = zebras * zebra_growth_rate;
	zebra_deaths  = min(zebras, meeting_chance*fighting_chance*zebras*lions);

	lion_births = 1/zebras_per_lion * zebra_deaths;
	lion_deaths = lions * lion_death_rate;

	# Update the stock levels
	zebras = zebras + zebra_births - zebra_deaths;
	lions = lions + lion_births - lion_deaths

	# Stock values always update in the next time step
	zebras_over_time[sim_step+1] = zebras;
	lions_over_time[sim_step+1] = lions;
end

In [16]:
using Gadfly
using DataFrames

In [40]:
d = DataFrame(ID = 1:size(d, 1), T = model_time, L=lions_over_time, Z=zebras_over_time);
[head(d),tail(d)]


Out[40]:
IDTLZ
110.2515.030000.0
220.518.031500.0
330.7521.8732130.0
441.026.70983131529.169
551.2532.4602356556043929413.61504439561
661.538.76194084828707525748.60929503161
739699.00.020858287929950391.7362620572629801e-6
839799.250.0187724591373175062.0831523141764537e-6
939899.50.0168952132239768152.4993917180947974e-6
1039999.750.0152056919020014112.9988477841536824e-6
11400100.00.0136851227122572663.59816134542975e-6
124010.00.0123166104415239534.317301201720193e-6

In [53]:
d[:ZN] = d[:Z]/1000
d


Out[53]:
IDTLZZN
110.2515.030000.030.0
220.518.031500.031.5
330.7521.8732130.032.13
441.026.70983131529.16931.529169000000003
551.2532.4602356556043929413.6150443956129.41361504439561
661.538.76194084828707525748.6092950316125.74860929503161
771.7544.86640746765506520917.6704498412320.91767044984123
882.049.76477397765659615716.19728304243115.716197283042431
992.2552.6094266356796111038.30668386224711.038306683862247
10102.553.155673828779487438.7781639668697.4387781639668695
11112.7551.7942391035862364972.4011390755434.972401139075543
12123.049.190232528389853391.4640317284193.391464031728419
13133.2545.9394783188747842401.4877947501882.4014877947501883
14143.542.448761451786991778.55438890053871.7785543889005386
15153.7538.958859616442971379.2909568459691.379290956845969
16164.035.6003296823785861117.79312063524871.1177931206352487
17174.2532.43823475025383943.41370864920210.9434137086492022
18184.529.500438028706146826.06969690133870.8260696969013387
19194.7526.794088404843833747.58945727330640.7475894572733064
20205.024.314989344446552696.79756864086570.6967975686408657
21215.2522.05291666456929666.73082780164980.6667308278016497
22225.519.99465859194445653.04339952988890.6530433995298889
23235.7518.125766530943235653.07828124263760.6530782812426376
24246.016.43156532237125665.31849296882770.6653184929688277
25256.2514.897731032908112689.05994878860440.6890599487886044
26266.513.510612227443321724.21764072030460.7242176407203046
27276.7512.257397241819447771.21493174390650.7712149317439065
28287.011.12618839540958830.92704032061040.8309270403206105
29297.2510.106020063803093904.66194045026040.9046619404502604
30307.59.186843374634277994.16901132881880.9941690113288187
&vellip&vellip&vellip&vellip&vellip&vellip

In [54]:
ds = stack(d, [:ZN, :L], [:ID, :T]);
[head(ds),tail(ds)]


Out[54]:
variablevalueIDT
1ZN30.010.25
2ZN31.520.5
3ZN32.1330.75
4ZN31.52916900000000341.0
5ZN29.4136150443956151.25
6ZN25.7486092950316161.5
7L0.0208582879299503939699.0
8L0.01877245913731750639799.25
9L0.01689521322397681539899.5
10L0.01520569190200141139999.75
11L0.013685122712257266400100.0
12L0.0123166104415239534010.0

In [55]:
plot(ds, x = :T, y = :value, color =:variable, Geom.line )


Out[55]:
T -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 ZN L variable -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 -600 -580 -560 -540 -520 -500 -480 -460 -440 -420 -400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 820 840 860 880 900 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 1140 1160 1180 1200 -1000 0 1000 2000 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 value

In [ ]: