The Lotka–Volterra equations, also known as the predator–prey equations, are a pair of first-order, non-linear, differential equations. They are frequently used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey. You can read more about this from Wikipedia http://en.wikipedia.org/wiki/Lotka-Volterra_equation.
Their populations change with time according to the following pair of equations:
$ dx/dt = x ( \alpha - \beta y ) $
$ dy/dt = -y( \gamma - \delta x ) $
Where, $x$ is the number of prey (say rabbits), $y$ is the number of predators (say foxes). $dx/dt, dy/dt$ gives the rate of change of their respective populations over time $t$. $\alpha, \beta, \gamma, \delta $ are the parameters describing the interaction between the two species. This models makes some assumptions about the predator and prey. You can read about them from the Wikipedia page mentioned above.
The above equations can be written in a slightly different form to interpret the physical meaning of the four parameters used.
1.Equation for prey
$dx/dt = \alpha x - \beta x y $
The prey are supposed to have unlimited supply of food and $\alpha x$ represents the rate of population growth of prey. Rate of decrease of population of prey is assumed to be proportional to the rate at which predator and prey meet and is given by $ \beta y x$
2.Equation for predator
$ dy/dt = \delta x y -\gamma y $
For the predators, $\delta x y $ gives the rate of growth of predator population. Note that this is similar to the rate of decrease of population of prey. The second term $\gamma y $ gives the rate of population decrease for predators due to natural death or emigration.
A simple python code for solving these equations is shown below.
In [2]:
# importrequired libraries
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
In [3]:
# set the initial parameters
alpha = 1.
beta = 1.2
gamma = 4.
delta = 1.
In [4]:
#define the time stepping scheme - euler forward, as used in earlier lessons
def euler_step(u, f, dt):
"""Returns the solution at the next time-step using Euler's method.
Parameters
----------
u : array of float
solution at the previous time-step.
f : function
function to compute the right hand-side of the system of equation.
dt : float
time-increment.
Returns
-------
u_n_plus_1 : array of float
approximate solution at the next time step.
"""
return u + dt * f(u)
In [5]:
# define the function that represents the Lotka-Volterra equations
def f(u):
"""Returns the rate of change of species numbers.
Parameters
----------
u : array of float
array containing the solution at time n.
Returns
-------
dudt : array of float
array containing the RHS given u.
"""
x = u[0]
y = u[1]
return numpy.array([x*(alpha - beta*y), -y*(gamma - delta*x)])
In [6]:
# set time-increment and discretize the time
T = 15.0 # final time
dt = 0.01 # set time-increment
N = int(T/dt) + 1 # number of time-steps
x0 = 10.
y0 = 2.
t0 = 0.
# set initial conditions
u_euler = numpy.empty((N, 2))
# initialize the array containing the solution for each time-step
u_euler[0] = numpy.array([x0, y0])
# use a for loop to call the function rk2_step()
for n in range(N-1):
u_euler[n+1] = euler_step(u_euler[n], f, dt)
In [7]:
time = numpy.linspace(0.0, T,N)
x_euler = u_euler[:,0]
y_euler = u_euler[:,1]
We will now plot the variation of population for each species with time.
In [8]:
plt.plot(time, x_euler, label = 'prey ')
plt.plot(time, y_euler, label = 'predator')
plt.legend(loc='upper right')
#labels
plt.xlabel("time")
plt.ylabel("number of each species")
#title
plt.title("predator prey model")
A better understanding of the system behaviour can be obtained by a phase plot of the population of predators vs. the population of prey. It will tell us if the system sustains or collapses over time. For the choice of parameters $ \alpha, \beta, \gamma $ and $ \delta $ made above, we see that the maximum population of each species keeps increasing each cycle. You can read more about that in the Wikipedia link mentioned above.
In [9]:
plt.plot(x_euler, y_euler, '-->', markevery=5, label = 'phase plot')
plt.legend(loc='upper right')
#labels
plt.xlabel("number of prey")
plt.ylabel("number of predators")
#title
plt.title("predator prey model")
Vary the interactions between species by changing $\alpha, \beta, \gamma, \delta$ and see what happens to the population evolution as well as phase plots.
Also, try the same exercise with a fourth order time stepping method called "Runge-Kutta 4" whose algorithm is given below.
In [10]:
def RK4(u,f,dt):
# Runge Kutta 4th order method
"""Returns the solution at the next time-step using Runge Kutta fourth order (RK4) method.
Parameters
----------
u : array of float
solution at the previous time-step.
f : function
function to compute the right hand-side of the system of equation.
dt : float
time-increment.
Returns
-------
u_n_plus_1 : array of float
approximate solution at the next time step.
"""
#calculate slopes
k1 = f(u)
u1 = u + (dt/2.)*k1
k2 = f(u1)
u2 = u + (dt/2.)*k2
k3 = f(u2)
u3 = u + dt*k3
k4 = f(u3)
return u + (dt/6.)*(k1 + 2.*k2 + 2.*k3 + k4)
In [2]:
from IPython.core.display import HTML
css_file = 'styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[2]: