This is an assignment in Python, I contributed to a numerical Python MOOC from George Washington University. The link to this assignment on github is here.
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.
Equations describing the system
Their populations change with time according to the following pair of equations:
Where, is the number of prey (say rabbits), is the number of predators (say foxes). gives the rate of change of their respective populations over time . 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
The prey are supposed to have unlimited supply of food and 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
2.Equation for predator
For the predators, 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 gives the rate of population decrease for predators due to natural death or emigration.
Numerical solution using Python
A simple python code for solving these equations is shown below.
# importrequired libraries import numpy import matplotlib.pyplot as plt %matplotlib inline
# set the initial parameters alpha = 1. beta = 1.2 gamma = 4. delta = 1.
#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)
# 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)])
# 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)
time = numpy.linspace(0.0, T,N) x_euler = u_euler[:,0] y_euler = u_euler[:,1]
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")