Lotka - Volterra Model

Model for the prey-predator population:

\begin{align} x' & = x (r-ay) \\ y' & = y(-m+bx) \end{align}

where $r$ is the per capita growth rate i.e., when there is no predator the prey dynamics is governed by $x'=rx$. In absence of prey the predator dies via $y'=-my$, where $m$ is the per capita mortality rate. The interaction between the prey and predator is proportional $xy$, change in the population of prey due to this interaction is given by $-axy$ and change in the population of predator due to this interaction is given by $bxy$. The parameter $a$ depends upon the number of interaction and success of those intearctions. $b=\epsilon a$, where $\epsilon$ is the conversion efficiency of the predator population.

Analysis of the Lotka-Voltera Model

Equilibrium Points: $(0,0)$ and $(m/b,r/a)$

Nullclines:
x-nullclines is where $x'=0$, $x=0$ and $y=r/a$
y-nullclines is where $y'=0$, $y=0$ and $x=m/b$

Purely vertical flow when $x'=0$ and $y'=r/a(-m+bx)$. This flow is upward for $x > \frac{m}{b}$ and downward for $x < \frac{m}{b}$.

Purely horizontal flow when $y'=0$ and $x'=m/b(r-ay)$. This flow is leftward for $y > \frac{r}{a}$ and rightward for $y < \frac{r}{a}$.

In [1]:
# import all the packages required 
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# lets assign some values to r, a, b and m
# you can assign any other values 
r=1; a = 0.1; m = 1.5; b = 0.2;


# x and y values where you want arrows
x=np.array([m/b,m/b,m/b-3,m/b+3])
y=np.array([r/a-3.0,r/a+3.0,r/a,r/a])


# velocity vector
u=x*(r-a*y); v=y*(-m+b*x) 

# normalize the lengths of the vectors 
u=u/np.sqrt(u**2 + v**2); v=v/np.sqrt(u**2 + v**2)

plt.quiver(x,y,u,v,color='b',scale=10,headlength=5)

plt.axhline(r/a,color='k',ls='--')
plt.axvline(m/b,color='k',ls='--')

plt.ylim(5,15.0);plt.xlim(0,16.0)
plt.text(7.1,3.5,'m/b',fontsize=18)
plt.text(6.1,14.0,'y-nullcline',fontsize=18)
plt.text(11.1,9.0,'x-nullcline',fontsize=18)

plt.text(-2.5,9.8,'r/a',fontsize=18)
plt.show;

Stability Analysis:

(i) Jacobian at $(0,0)$ is $$ A_{(0,0)}= \begin{bmatrix} r & 0 \\ 0 & -m \end{bmatrix} $$

The eigenvalues will be $\lambda_{1,2}=r,-m$. As $r$ and $m$ are positive numbers therefore, $(0,0)$ is a SADDLE. Orbits near $(0,0)$ will veer away.

(ii) Jacobian at $(m/b,r/a)$ is $$ A_{(m/b,r/a)}= \begin{bmatrix} 0 & \frac{-am}{b} \\ \frac{br}{a} & 0 \\\end{bmatrix} $$ The eigenvalues will be $\lambda_{1,2}=\pm i\sqrt{rm} $. Pure imaginary roots, we are not guaranteed to have a center in case of nonlinear systems.

In [3]:
# Solve the differential equations  
from scipy.integrate import odeint


def fun(X,t):
    # lets assign some values to r, a, b and m
    r=1; a = 0.1; m = 1.5; b = 0.075; 
    x=X[0];y=X[1];
    dx=x*(r-a*y); dy=y*(-1.0*m+b*x) 
    
    return [dx,dy]


plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 10, 500)

no_of_orbits=5;x0=5;y0=5;  
for i in np.arange(0,5):
    Xinit=np.array([0.5+x0,0.5+y0])
    Xsolu = odeint(fun,Xinit, t)


    plt.plot(Xsolu[:,0], Xsolu[:,1],lw=2)
    x0=x0+2.5;y0=y0+2.5;


# stuff required for plotting arrows/vectors on the plot    
xmax = plt.xlim(xmin=0)[1];ymax = plt.ylim(ymin=0)[1]       
xl = np.linspace(0.0,xmax,20)
yl = np.linspace(0.0,ymax,20)
x1 , y1  = np.meshgrid(xl, yl) 
[u,v]= fun([x1,y1],0); 


# normalize the lengths of the vectors 
uvn=np.sqrt(u**2 + v**2); uvn[uvn==0]=1.
u=u/uvn; v=v/uvn;

plt.quiver(x1,y1,u,v,uvn,scale=40,cmap=plt.cm.cool)
 
plt.xlim(-1,60); plt.ylim(-1,35)   
plt.xlabel('Prey',fontsize=20); plt.ylabel('Predator',fontsize=20)
plt.show;