Introduction

This example will illustrate using system dynamics with Julia by using the classic predator/prey problem.

The model simulates the interaction between foxes (preditors) and rabbits (their prey), over the course of 100 years.

Download the code

Background

System dynamics is one of the more popular forms of modeling, allowing us to better understand all sorts of systems in fields that vary from economics to environmental science.

One classic problem attempts to model the relationship between a predator and their prey. This can be described by a system of two first order non-linear differential equations, known as the Lotka-Volterra equation.

Where is the population size of the prey and is the population of the predator. In this example, the predator is the fox, who preys on the rabbit. As the population of foxes grows, the population of rabbits dwindles, and as a result, eventually rabbits die off. As rabbits die, foxes run out of food. Foxes starve, and their numbers drop. As their predators disappear, rabbit populations rebound, and so on, in cyclical fashion.

The rate of change in the population of foxes is inversely proportional to the population of foxes, and directly proportional to the population of rabbits. The effect is the opposite for the rabbits.

The solution to this problem is difficult to solve analytically, due to its non-linearity, which is why system dynamics is an ideal way to study the behavior of the system.

Setup

Let's get started. Open Julia Studio and create a new directory, called RabbitsFoxes, and create a file called main.jl. Open main.jl if it does not do so automatically.

The first thing we will do is initialize the parameters of the simulation:

# Time parameters
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 rabbits when the simulation is started. (Rabbits)
initial_rabbits = 30000

# The number of foxes when the simulation is started (Foxes)
initial_foxes = 15

# The number of rabbits that must be killed for a fox to be born. (Rabbits/Fox)
rabbits_killed_per_fox_birth = 1000

# The chance a rabbit will die when a rabbit fox cross paths. (dmnl)
chance_a_rabbit_will_die_during_a_meeting = 0.50

# The chance that a rabbit and fox will cross paths. (dmnl)
chance_of_rabbit_and_fox_meeting = 0.02

# The percent of the rabbit population that will be born this time step. (1/Year)
rabbit_growth_rate = 0.20

# The percent of the fox population that will die this time step from old age. (1/Year)
fox_death_rate = 0.10

# Initialization
rabbits_over_time = fill(0.0, end_step+1)
foxes_over_time = fill(0.0, end_step+1)
model_time = fill(0.0, end_step+1)

rabbits = initial_rabbits
foxes = initial_foxes

rabbits_over_time[1] = rabbits
foxes_over_time[1] = foxes

With these values set, the system is set up to oscillate. Feel free to experiment with the values of the initial parameters. Doing so will allow you to see how their values affect the outcomes. That's the beauty of simulation.

Simulation

Now let's write the equations that govern the simulation. First create a loop, starting at 1 and running until end_step, and print some output when the simulation is finished.

# Run the model
for sim_step = 1:end_step

end

println("Time,Rabbits (Thousands),Foxes")
for i = 1:end_step
    print(model_time[i])
    print(",")
    print(rabbits_over_time[i]/1000)
    print(",")
    println(foxes_over_time[i])
end

Now what equations do we need to write? First, for both rabbits and foxes, we need to figure out how much is added to the system, and how much leaves the system. In system dynamics, these are known as flows. If inflow exceeds outflow, then accumulation will occur. If outflow exceeds inflow, our stocks will decrease. If inflow and outflow are balanced, we are at steady state and the size of our stocks will remain constant. Within the loop, we will add the following calculations:

# Get the time from the step
sim_time = start_time + sim_step * time_step
model_time[sim_step] = sim_time

# First we must calculate our flows (our rates)
rabbit_births = rabbits * rabbit_growth_rate
rabbits_eaten = min(rabbits, chance_a_rabbit_will_die_during_a_meeting * chance_of_rabbit_and_fox_meeting * foxes * rabbits)

fox_births = 1/rabbits_killed_per_fox_birth * rabbits_eaten
fox_deaths = foxes * fox_death_rate

The next step is to update our stocks by adding our flows:

# Then we update our stocks
foxes = foxes + fox_births - fox_deaths
rabbits = rabbits + rabbit_births - rabbits_eaten

# Stock values always update in the next time step
rabbits_over_time[sim_step+1] = rabbits
foxes_over_time[sim_step+1] = foxes

main.jl should look like this:

# main.jl

#Time parameters
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 rabbits when the simulation is started. (Rabbits)
initial_rabbits = 30000

# The number of foxes when the simulation is started (Foxes)
initial_foxes = 15

# The number of rabbits that must be killed for a fox to be born. (Rabbits/Fox)
rabbits_killed_per_fox_birth = 1000

# The chance a rabbit will die when a rabbit fox cross paths. (dmnl)
chance_a_rabbit_will_die_during_a_meeting = 0.50

# The chance that a rabbit and fox will cross paths. (dmnl)
chance_of_rabbit_and_fox_meeting = 0.02

# The percent of the rabbit population that will be born this time step. (1/Year)
rabbit_growth_rate = 0.20

# The percent of the fox population that will die this time step from old age. (1/Year)
fox_death_rate = 0.10

rabbits_over_time = fill(0.0, end_step+1)
foxes_over_time = fill(0.0, end_step+1)
model_time = fill(0.0, end_step+1)

rabbits = initial_rabbits
foxes = initial_foxes

rabbits_over_time[1] = rabbits
foxes_over_time[1] = foxes

# Run the model
for sim_step = 1:end_step
    # Get the time from the step
    sim_time = start_time + sim_step * time_step
    model_time[sim_step] = sim_time

    # First we must calculate our flows (our rates)
    rabbit_births = rabbits * rabbit_growth_rate
    rabbits_eaten = min(rabbits, chance_a_rabbit_will_die_during_a_meeting * chance_of_rabbit_and_fox_meeting * foxes * rabbits)

    fox_births = 1/rabbits_killed_per_fox_birth * rabbits_eaten
    fox_deaths = foxes * fox_death_rate

    # Then we update our stocks
    foxes = foxes + fox_births - fox_deaths
    rabbits = rabbits + rabbit_births - rabbits_eaten

    # Stock values always update in the next time step
    rabbits_over_time[sim_step+1] = rabbits
    foxes_over_time[sim_step+1] = foxes
end

println("Time,Rabbits (Thousands),Foxes")
for i = 1:end_step
    print(model_time[i])
    print(",")
    print(rabbits_over_time[i]/1000)
    print(",")
    println(foxes_over_time[i])
end

When you run the code, you should see the output in the console. You should see oscillations. Experiment with changing parameters to see how you can affect the outcome.

Conclusion

Although it's a simplified example, this model describes the behavior of many natural systems. Julia computes the results quickly, efficiently, and simply.