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

    
import numpy
import matplotlib.pyplot as plt



In [49]:

    
# set the initial parameters
alpha = 1.
beta = 1.2
gamma = 4
delta = 1



In [50]:

    
#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 [51]:

    
# 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.
    """
    t = t0
    x = u[1]
    y = u[2]
    return numpy.array([1., x*(alpha - beta*y), -y*(gamma - delta*x)])



In [52]:

    
# 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, 3))

# initialize the array containing the solution for each time-step
u_euler[0] = numpy.array([t0, 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 [53]:

    
time = u_euler[:,0]
x_euler = u_euler[:,1]
y_euler = u_euler[:,2]



In [54]:

    
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")
plt.show()



In [55]:

    
plt.plot(x_euler, y_euler, 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")
plt.show()

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

    
def RK4(u,f,dt):
    # Runge Kutta 4th order method
    #u1 = numpy.zeros((len(u)))
    #u2 = numpy.zeros((len(u)))
    #u3 = numpy.zeros((len(u)))
    #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)