```
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from ipywidgets import interactive,IntSlider,FloatSlider,IntText,FloatText
from IPython.display import Javascript
def resize_colab_cell():
  display(Javascript('google.colab.output.setIframeHeight(0, true, {maxHeight: 5000})'))
get_ipython().events.register('pre_run_cell', resize_colab_cell)
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>


We denote with $X$ and $Y$ the number of preys and predators, respectively. We now make the following assumptions:
* The prey population *Y* have an unlimited food supply at all times.
*  In the absence of predators, the prey population $X$ would grow proportionally to its size with a constant of proportionality α>0 which Volterra called coefficient of auto-increase.
* In the absence of preys, the predator population $Y$ would decline proportionally to its size and the constant of proportionality is β>0.
* When both predator and prey are present, a decline in the prey population and a growth in the predator population will occur, each at a rate proportional to the frequency of encounters between individuals of the two species. These two constants of proportionality are indicate by δ>0 and γ>0.
    
We can translate these words into a system of two coupled ODEs, (remember $X=$ preys while $Y$= predators)
    \begin{cases}
      \frac{dX}{dt}=\alpha X- \delta XY\\
      \frac{dY}{dt}=-\beta Y + \gamma XY
    \end{cases}

The stationary points of this model can be found by rewriting the above equations as follows
    \begin{cases}
      \frac{dX}{dt}=(\alpha - \delta Y)X\\
      \frac{dY}{dt}=-(\beta - \gamma X)Y
    \end{cases}  
  and imposing the two equations to be equal to zero simultaneously
  \begin{cases}
      X(\alpha - \delta Y)=0\\
      Y(-\beta  +\gamma X)=0
 \end{cases}
  

ln this model there are two stationary points given by:

$P_{1}=(0,0)$
$P_{2}=(\beta/\gamma,\alpha/\delta)$

In order to determine the nature of the stationary point we perform a linear analysis of our system be defining the *Jacobian* matrix $J$ as 
\begin{equation}
J =
\begin{bmatrix}
  \frac{\partial F(X,Y)}{\partial X} & 
    \frac{\partial F(X,Y)}{\partial Y}\\[1ex] 
  \frac{\partial G(X,Y)}{\partial X} & 
    \frac{\partial G(X,Y)}{\partial Y}
\end{bmatrix}
\end{equation}
where we have defined 
\begin{equation}
    F(X,Y)= \alpha X- \delta XY,\quad\text{and}\quad  G(X,Y)=-\beta Y + \gamma XY.
\end{equation}


The Jacobian matrix becomes 
\begin{equation}
J(X,Y)=
\begin{bmatrix}
  \alpha-\delta Y& 
    -\delta X\\[1ex] 
  \gamma Y& -\beta +\gamma X
\end{bmatrix},
\end{equation}

\begin{equation}
J(0,0)=
\begin{bmatrix}
  \alpha& 
    0\\[1ex] 
  0& -\beta 
\end{bmatrix}.
\end{equation}


\begin{equation}
    \lambda_1=\alpha, \qquad \text{and}\qquad \lambda_2=-\beta.   \end{equation}
    
 The eigenvectors associated with the two eigenvalues are given by solving the following system of equations
\begin{equation}
\begin{bmatrix}
  \alpha & 
    0\\[1ex] 
  0& -\beta
\end{bmatrix}
\begin{bmatrix}
  x\\[1ex] 
  y\end{bmatrix}=\begin{bmatrix}
  \lambda x\\[1ex] 
  \lambda y\end{bmatrix},
\end{equation}
which means 
\begin{equation}
    \begin{cases}
      \alpha x = \lambda x\\
      -\beta y = \lambda y
    \end{cases}
    \qquad 
    \begin{cases}
     (\alpha -\lambda) x = 0\\
     (-\beta -\lambda)y = 0
    \end{cases}
\end{equation}

Therefore, if $\lambda=\alpha$ 
\begin{equation}
    \begin{cases}
      0=0\\
      (-\beta-\alpha)y =0
    \end{cases}
    \qquad 
    \begin{cases}
      x \quad \text{any value}\\
      y =0
    \end{cases}
\end{equation}
and if $\lambda=-\beta$ 
\begin{equation}
    \begin{cases}
      (\alpha + \beta)y =0 \\
      0=0
    \end{cases}
    \qquad 
    \begin{cases}
      x=0\\
      y \quad \text{any value}
    \end{cases}
\end{equation}
In conclusion, the eigenvectors associated with the two eigenvalues are given by 
\begin{equation}
    \textbf{e}_1=\left( {\begin{array}{c}
   1\\
   0\\
  \end{array} } \right)
\qquad \text{and} \qquad 
    \textbf{e}_2=\left( {\begin{array}{c}
   0\\
   1\\
  \end{array} } \right).
\end{equation}

Knowledge of eigenvalues and eigenvectors allows us to write down the general solution for our system
\begin{equation}
\left( {\begin{array}{c}
   X(t) \\
   Y(t) \\
  \end{array} } \right) =C_1e^{-\lambda_1 t}\mathbf{e_1} +C_2 e^{\lambda_2 t}\mathbf{e_2}=C_1e^{-\alpha t}\left( {\begin{array}{c}
   0 \\
   1 \\
  \end{array} } \right) +C_2 e^{\beta t}\left( {\begin{array}{c}
   1 \\
   0 \\
  \end{array} } \right).
\end{equation} 
We can see that if $t\to+\infty$, the second exponential diverges while the first one converges to zero.


```
def lotkaVolterra(N,t,a,d,b,g):
	dN0=a*N[0]-d*N[0]*N[1]
	dN1=g*N[0]*N[1]-b*N[1]
	return np.array([dN0,dN1])

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
a=2.0
d=0.06
b=1.0
g=0.05
T=np.linspace(1,40,10000)
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
fig, ax = plt.subplots()
ax.plot(0,0,"ko")
ax.plot(b/g,a/d,"bo")
ax.set_xlabel('X')
ax.set_ylabel('Y')
plt.show() 
```


```
fig, ax = plt.subplots()
ax.plot(0,0,"ko")
ax.plot(b/g,a/d,"bo")
 
ax.set_xlim(-1,70)
ax.set_ylim(-1,70)
my_xlim=ax.get_xlim()[1]
my_ylim=ax.get_ylim()[1]
dx=my_xlim/10
dy=my_ylim/10
[ax.arrow(0,6+i*dy,0,-dy,width=.3,head_width=1.5) for i in range(10)]
[ax.arrow(.1+i*dx,0,dx,0,width=.3,head_width=1.5) for i in range(10)]

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.grid()
plt.show() 
```

\begin{cases}
      \frac{dX}{dt}=(\alpha - \delta Y)X\\
      \frac{dY}{dt}=-(\beta + \gamma X)Y
    \end{cases}  

$\frac{dX}{dt}>0$ if $Y<\frac{\alpha}{\delta}$

and $\frac{dY}{dt}>0$ if $X<\frac{\beta}{\gamma}$



We now focus on the second stationary point $P_2=(\alpha/\delta , \beta/\gamma)$ and we evaluate the Jacobian 
\begin{equation}
J(\beta /\gamma,\alpha /\delta)=
\begin{bmatrix}
  \alpha-\delta\cdot \frac{\alpha}{\delta}& 
    -\delta \cdot \frac{\beta}{\gamma}\\[1ex] 
  \gamma \cdot \frac{\alpha}{\delta}& -\beta +\gamma \cdot\frac{\beta}{\gamma}
\end{bmatrix}=\begin{bmatrix}
  0& -\delta \cdot \frac{\beta}{\gamma}
    \\[1ex] 
  \gamma\cdot \frac{\alpha}{\delta}& 0 
\end{bmatrix}.
\end{equation}
The eigenvalues are obtained by evaluating the determinant of the following matrix
\begin{equation}
    M=\begin{bmatrix}
  -\lambda& -\delta \cdot \frac{\beta}{\gamma}
    \\[1ex] 
  \gamma\cdot \frac{\alpha}{\delta}& -\lambda 
\end{bmatrix},
\end{equation}
which is given by 
\begin{equation}
    \text{det}M = \lambda^2 + \gamma\frac{\alpha}{\delta}\cdot \delta\frac{\beta}{\gamma}=\lambda^2 +\alpha\beta,
\end{equation}
whose solutions are given by 
\begin{equation}
    \lambda_1=+i\sqrt{\alpha\beta}\qquad \text{and} \qquad \lambda_2=-i\sqrt{\alpha\beta}.
\end{equation}
The eigenvalues are complex and the real part is zero. We conclude that the point is a *center*.


```
fig, ax = plt.subplots()
ax.plot(0,0,"ko")
ax.plot(b/g,a/d,"bo")

ax.set_xlim(-1,70)
ax.set_ylim(-1,70)
my_xlim=ax.get_xlim()[1]
my_ylim=ax.get_ylim()[1]
dx=my_xlim/10
dy=my_ylim/10
[ax.arrow(0,6+i*dy,0,-dy,width=.3,head_width=1.5) for i in range(10)]
[ax.arrow(.1+i*dx,0,dx,0,width=.3,head_width=1.5) for i in range(10)]

plt.vlines(b/g,0,70,color='b',linestyle='dashed')
plt.hlines(a/d,0,70,color='b',linestyle='dashed')
plt.text(b/g,20,r'$\frac{\beta}{\gamma}$',fontsize=20)
plt.text(2,4+a/d,r'$\frac{\alpha}{\delta}$',fontsize=20)
ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$')
plt.ylabel('$predators\quad (Y)$');
```

\begin{cases}
      \frac{dX}{dt}=(\alpha - \delta Y)X\\
      \frac{dY}{dt}=-(\beta - \gamma X)Y
    \end{cases} 

$\frac{dX}{dt}>0$  if $Y\lt\frac{\alpha}{\delta}$

$\frac{dY}{dt}>0$ if $X\gt\frac{\beta}{\gamma}$





```
circ=lambda theta,r,hx,hy: [hx+r*np.cos(theta),hy+r*np.sin(theta)]
fig, ax = plt.subplots()
ax.plot(0,0,"ko")
ax.plot(b/g,a/d,"bo")
print(b/g,a/d)
ax.set_xlim(-1,70)
ax.set_ylim(-1,70)
my_xlim=ax.get_xlim()[1]
my_ylim=ax.get_ylim()[1]
dx=my_xlim/10
dy=my_ylim/10
[ax.arrow(0,6+i*dy,0,-dy,width=.3,head_width=1.5) for i in range(10)]
[ax.arrow(.1+i*dx,0,dx,0,width=.3,head_width=1.5) for i in range(10)]

plt.vlines(b/g,0,70,color='b',linestyle='dashed')
plt.hlines(a/d,0,70,color='b',linestyle='dashed')
plt.text(b/g,20,r'$\frac{\beta}{\gamma}$',fontsize=20)
plt.text(60,3+a/d,r'$\frac{\alpha}{\delta}$',fontsize=20)
plt.arrow(*circ(0.75*np.pi,15,b/g,a/d),0,-5,width=.3,head_width=1.5)
plt.arrow(*circ(-0.75*np.pi,15,b/g,a/d),0,-5,width=.3,head_width=1.5)
plt.arrow(*circ(0.25*np.pi,15,b/g,a/d),0,5,width=.3,head_width=1.5)
plt.arrow(*circ(-0.25*np.pi,15,b/g,a/d),0,5,width=.3,head_width=1.5)
#
plt.arrow(*circ(0.75*np.pi,15,b/g,a/d),-5,0,width=.3,head_width=1.5)
plt.arrow(*circ(-0.75*np.pi,15,b/g,a/d),5,0,width=.3,head_width=1.5)
plt.arrow(*circ(0.25*np.pi,15,b/g,a/d),-5,0,width=.3,head_width=1.5)
plt.arrow(*circ(-0.25*np.pi,15,b/g,a/d),5,0,width=.3,head_width=1.5)
ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$')
plt.ylabel('$predators\quad (Y)$')
```


```
N=np.array([50,1])
T=np.linspace(0,16,1000)
Nt=odeint(lotkaVolterra,N,T,args=(a,d,b,g))
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
ax=plt.subplot()
plt.plot(Nt[:,0],Nt[:,1])
plt.plot(Nt[0,0],Nt[0,1],'ro',markerSize=10)
plt.vlines(b/g,0,300,color='b',linestyle='dashed')
plt.hlines(a/d,0,300,color='b',linestyle='dashed')
ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$')
plt.ylabel('$predators\quad (Y)$')
```


```
N_ini=np.array([[16,26.666667],[44,73.3333333],[72,120],[100,166.66666667],[140,233.33333333]])
Nts=[]
for No in N_ini:
  Nt=odeint(lotkaVolterra,No,T,args=(a,d,b,g))
  Nts.append(Nt)
  
  
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
fig, ax = plt.subplots()
for Nt in Nts:
  ax.plot(Nt[:,0],Nt[:,1])
  ax.plot(Nt[0,0],Nt[0,1],'ro',markerSize=10)

ax.vlines(b/g,0,340,color='b',linestyle='dashed')
ax.hlines(a/d,0,340,color='b',linestyle='dashed')


ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$')
plt.ylabel('$predators\quad (Y)$')
```


```
ymax=100
xmax=100
n_points=25
x=np.linspace(0,xmax,n_points)
y=np.linspace(0,ymax,n_points)
X,Y=np.meshgrid(x,y)
dX,dY=lotkaVolterra([X,Y],0,a,d,b,g)

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
hypotenusa =lambda x,y : np.sqrt(x**2+y**2)
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
h=hypotenusa(dX,dY)
h[h==0]=1
dX=dX/h
dY=dY/h
fig, ax = plt.subplots()
Q = ax.quiver(X,Y,dX, dY, h, pivot='mid', cmap=plt.cm.jet)
plt.vlines(g/b,0,100,color='b',linestyle='dashed')
plt.hlines(a/d,0,100,color='b',linestyle='dashed')


ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$');
plt.ylabel('$predators\quad (Y)$');

```


```
N_ini=np.column_stack((X.diagonal()[4:20:4],Y.diagonal()[4:20:4]))
Nts=[]
for No in N_ini:
  Nt=odeint(lotkaVolterra,No,T,args=(a,d,b,g))
  Nts.append(Nt)
fig, ax = plt.subplots()
for Nt in Nts:
  ax.plot(Nt[:,0],Nt[:,1],color='r',alpha=0.8)

Q = ax.quiver(X,Y,dX, dY, h, pivot='mid', cmap=plt.cm.jet)

ax.set_aspect('equal')
plt.xlabel('$preys\quad (X)$')
plt.ylabel('$predators\quad (Y)$');
```


```
N=np.array([50,1])
Nt=odeint(lotkaVolterra,N,T,args=(a,d,b,g))
plt.plot(T,Nt[:,0])
plt.plot(T,Nt[:,1])
plt.xlabel('time')
plt.ylabel('N')
plt.legend(['Prey','Predator'])


```

**Species average number**

\begin{equation}
    \frac{dX}{dt}=(\alpha -\delta Y(t))X,
\end{equation}

\begin{equation}
    \frac{1}{X}\frac{dX}{dt}=\alpha -\delta Y(t),
\end{equation}

\begin{equation}
    \int_{X(0)}^{X(T)}\frac{dX}{X}=\int_0^T(\alpha -\delta Y(t))dt.
\end{equation}
The left-hand side can be immediately integrated and it gives
\begin{equation}
    \log X(T) - \log X(0)=0,
\end{equation}
which gives zero because the system is periodic and therefore the number of predator at time $t=0$ is equal to the number of predators at time $t=T$.
The integral on the right-hand side gives 
\begin{equation}
    \int_0^T(\alpha -\delta Y(t))dt=\alpha T - \delta\int_0^TY(t)dt,
\end{equation}
which cannot be evaluated exactly since we do not know the expression of $Y(t)$. By combining these two results we have
\begin{equation}
    0=\alpha T -\delta\int_0^TY(t)dt, 
\end{equation}
which can be written as 
\begin{equation}
    \frac{1}{T}\int_0^TY(t)dt=\frac{\alpha}{\delta}.
\end{equation}
The left-hand side is the average number of preys over a period of length $T$. If we consider the second equation, by doing the same calculations we obtain
\begin{equation}
    \frac{1}{T}\int_0^TX(t)dt=\frac{\beta}{\gamma}.
\end{equation}


```
def update(N0X,N0Y,T,alpha,beta,delta,gamma):
    N0=np.array([N0X,N0Y])
    T=np.linspace(0,T,1000)
    sol=odeint(foo,N0,T,args=(alpha,delta,beta,gamma))
    fig,axes=plt.subplots(nrows=1,ncols=2,figsize=(25,10))
    [axes[0].plot(T,sol[:,i]) for i in (0,1)]
    axes[0].grid()
    axes[0].set_title(r"$\overline{X}\quad (prey) = %.3f\quad\overline{Y}\quad(predator) = %.3f$"%(beta/gamma,alpha/delta))
    axes[0].set_xlabel("time",fontsize=20)
    axes[0].set_ylabel("species",fontsize=20)
    axes[0].legend(['prey','predator'],fontsize=20)
    axes[1].plot(sol[:,0],sol[:,1])
    axes[1].set_xlabel("predator",fontsize=20)
    axes[1].set_ylabel("prey",fontsize=20)
    axes[1].grid()
    plt.show()
timeSlider=FloatText(value=7.72,min=1,max=100,step=1,description="T")
N0XSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0X")
N0YSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0Y")
alphaSlider=FloatSlider(value=2.5, min=0.1,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&alpha;</MATH>')
betaSlider=FloatSlider(value=1, min=0.01,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&beta;</MATH>')
deltaSlider=FloatSlider(value=0.07, min=0.01,max=0.1,step=0.01,readout_format='.3f',description=r'<MATH>&delta;</MATH>')
gammaSlider=FloatSlider(value=0.14, min=0.,max=0.3,step=0.01,description=r'<MATH>&gamma;</MATH>')

foo=lotkaVolterra
interactive_plot = interactive(update, N0X=N0XSlider,N0Y=N0YSlider, T=timeSlider,
                               alpha=alphaSlider, beta=betaSlider,gamma=gammaSlider,delta=deltaSlider)
output = interactive_plot.children[-1]
output.layout.height = '650px'
interactive_plot
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    interactive(children=(FloatSlider(value=1.0, description='N0X', max=10.0, step=1.0), FloatSlider(value=1.0, de…


  \begin{cases}
      \frac{dX}{dt}=(\alpha-\rho) X- \delta XY\\
      \frac{dY}{dt}=-(\beta+\psi) Y + \gamma XY
    \end{cases}


```
def lotkaVolterra2(N,t,a,d,b,g,r,p):
	dN0=(a-r)*N[0]-d*N[0]*N[1]
	dN1=g*N[0]*N[1]-(b+p)*N[1]
	return np.array([dN0,dN1])
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
def update2(N0X,N0Y,T,alpha,beta,delta,gamma,rho,psi):
    N0=np.array([N0X,N0Y])
    T=np.linspace(0,T,1000)
    sol=odeint(foo,N0,T,args=(alpha,delta,beta,gamma,rho,psi))
    fig,axes=plt.subplots(nrows=1,ncols=2,figsize=(25,10))
    [axes[0].plot(T,sol[:,i]) for i in (0,1)]
    axes[0].grid()
    axes[0].set_title(r"$\overline{X}\quad (prey) = %.3f\quad\overline{Y}\quad(predator) = %.3f$"%((beta+psi)/gamma,(alpha-rho)/delta))
    axes[0].set_xlabel("time",fontsize=20)
    axes[0].set_ylabel("species",fontsize=20)
    axes[0].legend(['prey','predator'],fontsize=20)
    axes[1].plot(sol[:,0],sol[:,1])
    axes[1].set_xlabel("predator",fontsize=20)
    axes[1].set_ylabel("prey",fontsize=20)
    axes[1].grid()
    plt.show()
timeSlider=FloatText(value=7.72,min=1,max=100,step=1,description="T")
N0XSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0X")
N0YSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0Y")
alphaSlider=FloatSlider(value=2.5, min=0.1,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&alpha;</MATH>')
betaSlider=FloatSlider(value=1, min=0.01,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&beta;</MATH>')
deltaSlider=FloatSlider(value=0.07, min=0.01,max=0.1,step=0.01,readout_format='.3f',description=r'<MATH>&delta;</MATH>')
gammaSlider=FloatSlider(value=0.14, min=0.,max=0.3,step=0.01,description=r'<MATH>&gamma;</MATH>')
rhoSlider=FloatSlider(value=0., min=0.,max=10,step=0.01,description=r'<MATH>&rho;</MATH>')
psiSlider=FloatSlider(value=0., min=0.,max=10,step=0.01,description=r'<MATH>&psi;</MATH>')
foo=lotkaVolterra2
interactive_plot = interactive(update2, N0X=N0XSlider,N0Y=N0YSlider, T=timeSlider,
                               alpha=alphaSlider, beta=betaSlider,gamma=gammaSlider,delta=deltaSlider,rho=rhoSlider,psi=psiSlider)
output = interactive_plot.children[-1]
output.layout.height = '650px'
interactive_plot
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    interactive(children=(FloatSlider(value=1.0, description='N0X', max=10.0, step=1.0), FloatSlider(value=1.0, de…


# Predator Prey with logistic growth

$
\frac{dX}{dt}=\alpha X(1-\frac{X}{\omega})-\beta XY\\
$

$
\frac{dY}{dt}=\gamma XY - \delta Y\\
$


```
def lotkaVolterraMod(N,t,a,w,b,d,g):
	dN0=a*N[0]*(1-N[0]/w)-b*N[0]*N[1]
	dN1=g*N[0]*N[1]-d*N[1]
	return np.array([dN0,dN1])

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
a=2
w=4
d=1
b=2
g=1
No=np.array([3,1])
T=np.linspace(1,100,10000)
print(d<g*w,a*(g*w-d)/(b*g*w),d/g)

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>


    True 0.75 1.0



```
Nt=odeint(lotkaVolterraMod,No,T,args=(a,w,b,d,g))
Nt2=odeint(lotkaVolterraMod,No,T,args=(a-0.5,w,b,d,g))
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
fig, ax = plt.subplots()
ax.plot(Nt[:,0],Nt[:,1],color='r',alpha=0.8)
#ax.plot(Nt2[:,0],Nt2[:,1],color='b',alpha=0.8)
X=np.linspace(np.min(Nt[:,0]),np.max(Nt[:,0]))
plt.plot(X,a*(w - X)/(b*w))
plt.vlines(d/g,0,np.max(Nt[:,1]))
ax.set_xlabel('preys');
ax.set_ylabel('predators');
ax.grid()

```


```
print(Nt[-1])
fig, ax = plt.subplots()
ax.plot(T,Nt[:,0],color='r')
ax.plot(T,Nt[:,1],color='k')
ax.set_xlabel('time');
ax.set_ylabel('N');
ax.grid()
```


```
ymax=6
xmax=3
n_points=25
x=np.linspace(0,xmax,n_points)
y=np.linspace(0,ymax,n_points)
X,Y=np.meshgrid(x,y)
dX,dY=lotkaVolterraMod([X,Y],0,a,w,d,b,g)

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```

h=hypotenusa(dX,dY)
h[h==0]=1
dX=dX/h
dY=dY/h
fig, ax = plt.subplots()
Q = ax.quiver(X,Y,dX, dY, h, pivot='mid', cmap=plt.cm.jet)
ax.set_xlabel('preys');
ax.set_ylabel('predators');
```


```
Nts=[]
N_ini=np.array([[3,1],[3,2],[3,3],[3,4]])
for No in N_ini:
  Nt=odeint(lotkaVolterraMod,No,T,args=(a,w,d,b,g))
  Nts.append(Nt)
```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



```
fig, ax = plt.subplots()
for Nt in Nts:
  ax.plot(Nt[:,0],Nt[:,1],color='r',alpha=0.8)

Q = ax.quiver(X,Y,dX, dY, h, pivot='mid', cmap=plt.cm.jet)
ax.set_xlabel('preys');
ax.set_ylabel('predators');
```


```
a=2
w=4
d=1
b=2
g=1
No=np.array([3,1])
T=np.linspace(1,100,10000)
print(d<g*w,a*(g*w-d)/(b*g*w),d/g)

```


    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>



    <IPython.core.display.Javascript object>


    True 0.75 1.0



```
def update3(N0X,N0Y,T,alpha,omega,beta,delta,gamma):
    N0=np.array([N0X,N0Y])
    T=np.linspace(0,T,1000)
    sol=odeint(foo,N0,T,args=(alpha,omega,delta,beta,gamma))
    fig,axes=plt.subplots(nrows=1,ncols=2,figsize=(25,10))
    SSy=(alpha/delta)*(1-beta/(omega*gamma))
    SSx=beta/gamma
    axes[0].hlines(SSx,0,T[-1],color='k',lw=2,ls=':')
    axes[0].hlines(SSy,0,T[-1],color='r',lw=2,ls=':')

    [axes[0].plot(T,sol[:,i]) for i in (0,1)]
    axes[0].set_xlabel("time",fontsize=20)
    axes[0].set_ylabel("species",fontsize=20)
    axes[0].legend([r'prey',r'predator',r'$SS_{prey}$',r'$SS_{predator}$'],fontsize=20)
    axes[0].grid()

    axes[1].plot(sol[:,0],sol[:,1])
    xlim=axes[1].get_xlim()
    ylim=axes[1].get_ylim()
    axes[1].vlines(SSx,ylim[0],ylim[1],color='k',lw=2,ls=':',label=r'$SS_{prey}$')
    axes[1].hlines(SSy,xlim[0],xlim[1],color='r',lw=2,ls=':',label=r'$SS_{predator}$')
    axes[1].legend()
    axes[1].set_xlabel("predator",fontsize=20)
    axes[1].set_ylabel("prey",fontsize=20)
    axes[1].grid()

    plt.show()

timeSlider=FloatText(value=20,min=1,step=1,description="T")
N0XSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0X")
N0YSlider=FloatSlider(value=1, min=0.,max=10,step=1,description="N0Y")
alphaSlider=FloatSlider(value=2, min=0.,max=4,step=0.1,readout_format='.3f',description=r'<MATH>&alpha;</MATH>')
omegaSlider=FloatSlider(value=4, min=0.,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&omega;</MATH>')

betaSlider=FloatSlider(value=1, min=0.01,max=4,step=0.1,readout_format='.3f',description=r'<MATH>&beta;</MATH>')
deltaSlider=FloatSlider(value=1, min=0.01,max=2,step=0.01,readout_format='.3f',description=r'<MATH>&delta;</MATH>')
gammaSlider=FloatSlider(value=1, min=0.,max=2,step=0.01,description=r'<MATH>&gamma;</MATH>')

#N0XSlider=FloatText(value=1, min=0.,max=10,step=1,description="N0X")
#N0YSlider=FloatText(value=1, min=0.,max=10,step=1,description="N0Y")
##alphaSlider=FloatText(value=2, min=0.,max=4,step=0.1,readout_format='.3f',description=r'<MATH>&alpha;</MATH>')
#omegaSlider=FloatText(value=4, min=0.,max=10,step=0.1,readout_format='.3f',description=r'<MATH>&omega;</MATH>')

#betaSlider=FloatText(value=1, min=0.01,max=4,step=0.1,readout_format='.3f',description=r'<MATH>&beta;</MATH>')
#deltaSlider=FloatText(value=1, min=0.01,max=2,step=0.01,readout_format='.3f',description=r'<MATH>&delta;</MATH>')
#gammaSlider=FloatText(value=1, min=0.,max=2,step=0.01,description=r'<MATH>&gamma;</MATH>')



foo=lotkaVolterraMod
interactive_plot = interactive(update3, N0X=N0XSlider,N0Y=N0YSlider, T=timeSlider,
                               alpha=alphaSlider,omega=omegaSlider, beta=betaSlider,gamma=gammaSlider,delta=deltaSlider)
output = interactive_plot.children[-1]
output.layout.height = '650px'


interactive_plot
```

491483b4c3ee97b64eb56a2946ca683ffffcfa79
ipynb
Jupyter Notebook
GrowthModels/predator_prey.ipynb
cohmathonc/biosci670coh
514a11088a4feb87e668a8ed6a5da11283676f8c
CC-BY-4.0
2021-02-23T21:12:57.000Z
2021-02-23T21:12:57.000Z
GrowthModels/predator_prey.ipynb
cohmathonc/biosci670coh
514a11088a4feb87e668a8ed6a5da11283676f8c
CC-BY-4.0
GrowthModels/predator_prey.ipynb
cohmathonc/biosci670coh
514a11088a4feb87e668a8ed6a5da11283676f8c
CC-BY-4.0
2021-03-09T19:24:38.000Z
2021-03-09T19:24:38.000Z
Qwen/Qwen-72B
1. YES
2. YES