This is an example of an environmental system where the state variables are two species populations, one of which eats the other one. Say we have a population of rabbits $R$ and foxes $F$. These are our "stocks". The ODEs describing this system can be written as:
$\frac{dR}{dt} = aR - bRF$
$\frac{dF}{dt} = cRF - dF$
The coefficient $a$ is the number of births per rabbit per unit time. $b$ describes the death rate of rabbits per encounter with foxes (essentially, rabbit deaths per time per fox). $c$ is the efficiency with which new foxes are created by eating rabbits (fox births per rabbit eaten per time). And finally $d$ is the natural death rate of foxes in the absence of a food source.
Notice some important assumptions--rabbits only die by being eaten, and foxes are only born when rabbits are eaten. If the foxes have another food source, or the rabbits have another method of dying off, you might need some new equations.
With the preamble out of the way, let's get started...
In [2]:
import sys
sys.path.append('./../')
from stockflow import simulation
import matplotlib.pyplot as plt
import numpy as np
tmax = 20
dt = 0.1
t = np.arange(0,tmax,dt)
(These path hacks are just because I don't have stockflow installed globally...). Now initialize states and choose some parameter values.
In [3]:
s = simulation(t)
s.stocks({'rabbits': 10, 'foxes': 5})
a = 0.5
b = 0.3
c = 0.3
d = 0.9
In [4]:
s.flow('rabbit_births', start=None, end='rabbits', f=lambda t: a*s.rabbits)
s.flow('rabbit_deaths', start='rabbits', end=None, f=lambda t: b*s.rabbits*s.foxes)
s.flow('fox_births', start=None, end='foxes', f=lambda t: c*s.rabbits*s.foxes)
s.flow('fox_deaths', start='foxes', end=None, f=lambda t: d*s.foxes)
s.run() # not using discrete=True here because the flow functions are continuous
If you're familiar with this example, this might seem like a verbose way of setting it up, and it is. But if you have a lot of stocks/flows this could be an intuitive way of setting up a system and accessing variables, kind of like combining STELLA with Python.
Let's check out the dynamics.
In [5]:
plt.plot(t, s.rabbits, t, s.foxes)
plt.xlabel('Time (units)')
plt.ylabel('Population')
plt.legend(["Rabbits","Foxes"])
plt.show()
These equations are interesting because they usually show oscillating behavior. Here the populations both come precariously close to extinction--of course, our continuous equations allow populations less than 1, which doesn't make much sense.
What happens if rabbits become better at escaping and/or foxes become less efficient hunters (or digesters)?
In [6]:
b = 0.1
c = 0.1
s.run()
In [7]:
plt.plot(t, s.rabbits, t, s.foxes)
plt.xlabel('Time (units)')
plt.ylabel('Population')
plt.legend(["Rabbits","Foxes"])
plt.show()
With these coefficients the populations are more stable, but still oscillating. Let's end with a phase plot of the two state variables. This is an interesting application of quiver from StackOverflow:
In [9]:
plt.quiver(s.rabbits[:-1], s.foxes[:-1], s.rabbits[1:]-s.rabbits[:-1], s.foxes[1:]-s.foxes[:-1], scale_units='xy', angles='xy', scale=1)
plt.xlabel('Rabbit Population')
plt.ylabel('Fox Population')
plt.show()
Just going around in circles...