### Tutorial 2 – Rabbits and Foxes

#### 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.

#### 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.