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 [39]:

    

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


    

    
Out[39]:




T
L
Z
id
1
0.25
15.0
30000.0
1
2
0.5
18.0
31500.0
2
3
0.75
21.87
32130.0
3
4
1.0
26.709831
31529.169
4
5
1.25
32.46023565560439
29413.61504439561
5
6
1.5
38.761940848287075
25748.60929503161
6
7
99.0
0.02085828792995039
1.7362620572629801e-6
396
8
99.25
0.018772459137317506
2.0831523141764537e-6
397
9
99.5
0.016895213223976815
2.4993917180947974e-6
398
10
99.75
0.015205691902001411
2.9988477841536824e-6
399
11
100.0
0.013685122712257266
3.59816134542975e-6
400
12
0.0
0.012316610441523953
4.317301201720193e-6
401

In [38]:

    

ds = stack(d, [:Z, :L], :T);
[head(ds),tail(ds)]


    

    
Out[38]:




variable
value
T
1
Z
30000.0
0.25
2
Z
31500.0
0.5
3
Z
32130.0
0.75
4
Z
31529.169
1.0
5
Z
29413.61504439561
1.25
6
Z
25748.60929503161
1.5
7
L
0.02085828792995039
99.0
8
L
0.018772459137317506
99.25
9
L
0.016895213223976815
99.5
10
L
0.015205691902001411
99.75
11
L
0.013685122712257266
100.0
12
L
0.012316610441523953
0.0

In [ ]: