Systems of Ordinary Differential Equations

Objectives:

visualization of systems of ODE's, trajectiory, vector fields...
Classified Stabilities at equilibrium points.

Requirement:

IPython
NumPy/SciPy
Sympy
Bokeh: observation for Equilibrium points (?)



In [2]:

    
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D


from matplotlib import animation
from JSAnimation import IPython_display

from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets



In [29]:

    
#from directionarrow import *

2D differential Equations

\begin{eqnarray} \frac{ d x}{ d t} &=& 2 x (1-x/2 ) - x y\\ \frac{ d y}{ d t} &=& 3 y (1 -y/3) - 2 x y \end{eqnarray}



In [1]:

    
def f(x,y): return 2.*x*(1-x/2.) - x*y
def g(x,y): return 3.*y*(1.-y/3.) - 2.*x*y



In [30]:

    
fig1 = plt.figure(figsize=(9,6))
ax = fig1.add_subplot(1,1,1)
xmin=-0.5;  xmax=3.5;  ymin=-0.5; ymax=4.5
xmesh=linspace(xmin,xmax,25); ymesh=linspace(ymin,ymax,21)

plotwidth=5.0 # inches, horizontal distance between xmin and xmax on display
plotheight=3.9 # inches, vertical distance between ymin and ymax on display
alpha= plotwidth/(xmax-xmin);  beta = plotheight/(ymax-ymin) # scaling factors

for xp in xmesh:
  for yp in ymesh:
    vx=f(xp,yp); vy=g(xp,yp) # direction vector v=(vx,vy)
    ax.quiver(xp, yp, vx, vy, pivot='middle', headwidth=3, headlength=4,color="red")

plt.scatter((0,1,2,0),(3,1,0,0),c='b',marker='o',alpha=0.6)

ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_aspect(beta/alpha)
ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
plt.grid()                             # turn on a grid of light, dotted lines
plt.autoscale(enable=True,tight=True)  # remove extra whitespace around slopefield
plt.title("Direction Field for\nx'= 2x(1-x/2)-xy, y'= 3y(1-y/3)-2xy")









    Out[30]:





<matplotlib.text.Text at 0x11170f650>

2D differential Equations



In [20]:

    
fig1 = plt.figure(figsize=(9,6))
ax = fig1.add_subplot(1,1,1)
xmin=-0.5;  xmax=3.5;  ymin=-0.5; ymax=4.5
xmesh=linspace(xmin,xmax,25); ymesh=linspace(ymin,ymax,21)

plotwidth=5.0 # inches, horizontal distance between xmin and xmax on display
plotheight=3.9 # inches, vertical distance between ymin and ymax on display
alpha= plotwidth/(xmax-xmin);  beta = plotheight/(ymax-ymin) # scaling factors

for xp in xmesh:
  for yp in ymesh:
    vx=f(xp,yp); vy=g(xp,yp) # direction vector v=(vx,vy)
    #xda,yda=directionarrow([xp,yp],[vx,vy],[alpha,beta],0.1)
    #ax.plot(xda,yda,'b')
    h=0.1/sqrt(1.+vx**2+vy**2)
    
    ax.arrow(xp,yp,1.5*h*vx,1.5*h*vy, head_width=0.05, head_length=0.05, fc='r', ec='r')    


ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_aspect(beta/alpha)
ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
plt.grid()                             # turn on a grid of light, dotted lines
#plt.autoscale(enable=True,tight=True)  # remove extra whitespace around slopefield
plt.title("306ch4a Direction Field for\nx'= 2x(1-x/2)-xy, y'= 3y(1-y/3)-2xy")

plt.scatter((0,1,2,0),(3,1,0,0),c='b',marker='o')
#plt.scatter(0,3,1,1,2,0,0,0)









    Out[20]:





<matplotlib.collections.PathCollection at 0x108e1e5d0>



In [ ]:



In [31]:

    
from bokeh.plotting import output_notebook
from bokeh.plotting import line,hold,show,figure
output_notebook()









    




    
        
        
    
        
        BokehJS successfully loaded.



In [16]:

    
xmin=-0.5;  xmax=3.5;  ymin=-0.5; ymax=4.5
xmesh=linspace(xmin,xmax,25); ymesh=linspace(ymin,ymax,21)

plotwidth=5.0 # inches, horizontal distance between xmin and xmax on display
plotheight=3.9 # inches, vertical distance between ymin and ymax on display
alpha= plotwidth/(xmax-xmin);  beta = plotheight/(ymax-ymin) # scaling factors
hold()
for xp in xmesh:
  for yp in ymesh:
    vx=f(xp,yp); vy=g(xp,yp) # direction vector v=(vx,vy)
    #xda,yda=directionarrow([xp,yp],[vx,vy],[alpha,beta],0.1)
    #ax.plot(xda,yda,'b')

    h=0.1/sqrt(1.+vx**2+vy**2)
    
    #ax.arrow(xp,yp,2*h*vx,2*h*vy, head_width=0.05, head_length=0.05, fc='r', ec='r')  
    line((xp,xp+2*h*vx),(yp,yp+2*h*vy))
    
#ax.set_xlabel('x'); ax.set_ylabel('y')
#ax.set_aspect(beta/alpha)
#ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
#plt.grid()                             # turn on a grid of light, dotted lines
#plt.autoscale(enable=True,tight=True)  # remove extra whitespace around slopefield
#plt.title("306ch4a Direction Field for\nx'= 2x(1-x/2)-xy, y'= 3y(1-y/3)-2xy")
show()



In [ ]:

Equilibrium

$f(x,y)=g(x,y)=0$ implies $(x,y)=(0, 0), (0, 3), (1, 1), (2, 0)$ are equilibrium points.

What kind of the equilibrium points are, stable or unstable?

Condider the Jacobian: $$ J=\frac{\partial(f,g)}{\partial(x,y)}=\left(\begin{array}{cc} f_{x} & f_{y}\\ g_{x} & g_{y} \end{array}\right) $$ and $|J|$ repesents the determinant of $J$ with eigenvalues, $\lambda_1$ and $\lambda_2$:

Stability

$2\times2$ Jacobian Matrix, $\lambda_1,\lambda_2$ eigenvalues at equilibrium point $P=(x_0,y_0)$:

$\lambda_1<\lambda_2<0$: $P$ stable (sink);
$0<\lambda_1<\lambda_2$: $P$ unstable (source);
$\lambda_1<0<\lambda_2$: $P$ unstable (saddle point);
$\lambda=\pm b i$: $P$ stable (center of orbit);
$\lambda=a \pm b i$:
- $a<0$: asymptotically stable;
- $a>0$: unstable;
$\lambda_1=\lambda_2<0$:
- two independent eigenvectors: stable proper node;
- only one eigenvector: stable improper node;
$\lambda_1=\lambda_2>0$:
- two independent eigenvectors: unstable proper node;
- only one eigenvector: unstable improper node;



In [2]:

    
import numpy as np
from sympy import Symbol , diff, dsolve , Function , Derivative , Eq, lambdify, symbols,solve
from sympy.matrices import Matrix



In [3]:

    
x,y = symbols("x y")
eq1=Eq(2*x*(1-x/2)-x*y)
eq2=Eq(3*y*(1-y/3)-2*x*y)



In [4]:

    
eq_point=solve([eq1,eq2],[x,y]);print(eq_point)









    



[(0, 0), (0, 3), (1, 1), (2, 0)]



In [153]:



In [5]:

    
def f(x,y): return 2.*x*(1-x/2.) - x*y
def g(x,y): return 3.*y*(1.-y/3.) - 2.*x*y



In [7]:

    
def Jacobian_np(funcs,vars,eq_pt):
    """
    define the Jacobian (numeric verion)
    
    Input:
    func: vector field, e.g. [f,g] where x'=f(x,y),y'=g(x,y)
    vars: variables, e.g. [x,y]
    eq_pt: equilibrium point, e.g. [x0,y0]
    
    Output:
    out Jacobian matrix: d(f,g)/d(x,y)
    """
    J11=diff(funcs[0],vars[0]).subs({x:eq_pt[0],y:eq_pt[1]})
    J12=diff(funcs[0],vars[1]).subs({x:eq_pt[0],y:eq_pt[1]})
    J21=diff(funcs[1],vars[0]).subs({x:eq_pt[0],y:eq_pt[1]})
    J22=diff(funcs[1],vars[1]).subs({x:eq_pt[0],y:eq_pt[1]})
    A=np.zeros([2,2])
    A[0,0]=J11; A[0,1]=J12 # partial f partial x, partial f partial y
    A[1,0]=J21;  A[1,1]=J22   # partial g partial x, partial g partial y
    return A



In [5]:



In [8]:

    
A=Jacobian_np([f(x,y),g(x,y)],[x,y],[0,0])



In [9]:

    
from scipy import linalg as LA



In [10]:

    
eigvals,eigvecs=LA.eig(A);print(eigvals)









    



[ 2.+0.j  3.+0.j]



In [122]:



In [11]:

    
def Jacobian_smp(funcs,vars):
    """
    define the Jacobian (symbolic verion)
    
    Input:
    func: vector field, e.g. [f,g] where x'=f(x,y),y'=g(x,y)
    vars: variables, e.g. [x,y]
    
    Output:
    out Jacobian matrix: d(f,g)/d(x,y)
    """
    J11=diff(funcs[0],vars[0])
    J12=diff(funcs[0],vars[1])
    J21=diff(funcs[1],vars[0])
    J22=diff(funcs[1],vars[1])
    A=Matrix(([J11,J12],[J21,J22]))
    return A



In [13]:

    
def Jacobian_smp2(funcs,vars):
    """
    define the Jacobian (symbolic verion)
    
    Input:
    func: vector field, e.g. [f,g] where x'=f(x,y),y'=g(x,y)
    vars: variables, e.g. [x,y]
    
    Output:
    out Jacobian matrix: d(f,g)/d(x,y)
    """   
    dim=len(funcs)
    A=Matrix(dim,dim,lambda i,j:diff(funcs[i],vars[j]))
    return A



In [16]:



In [14]:

    
B=Jacobian_smp([f(x,y),g(x,y)],[x,y]);print(B)









    



Matrix([[-2.0*x - y + 2.0, -x], [-2.0*y, -2.0*x - 2.0*y + 3.0]])



In [16]:

    
B2=Jacobian_smp2([f(x,y),g(x,y)],[x,y]);print(B2)









    



Matrix([[-2.0*x - y + 2.0, -x], [-2.0*y, -2.0*x - 2.0*y + 3.0]])



In [134]:

    
MM = Matrix(3, 2, lambda i,j:Symbol('M_%d%d' % (i+1,j+1)));print MM









    



Matrix([[M_11, M_12], [M_21, M_22], [M_31, M_32]])



In [17]:

    
lam=B.eigenvals()
tt=lam.keys();print(tt)









    



[-2*x - 3*y/2 + sqrt(8*x*y + y**2 - 2*y + 1)/2 + 5/2, -2*x - 3*y/2 - sqrt(8*x*y + y**2 - 2*y + 1)/2 + 5/2]



In [18]:

    
lam=B2.eigenvals()
tt=lam.keys();print(tt)









    



[-2*x - 3*y/2 + sqrt(8*x*y + y**2 - 2*y + 1)/2 + 5/2, -2*x - 3*y/2 - sqrt(8*x*y + y**2 - 2*y + 1)/2 + 5/2]



In [19]:

    
C=B.subs({x:0,y:0});print(C)









    



Matrix([[2.00000000000000, 0], [0, 3.00000000000000]])



In [20]:

    
lam=C.eigenvals()
lam.keys()









    Out[20]:





[2, 3]



In [21]:

    
M =Matrix( [ [0, -1], [1, 0] ])
M.eigenvals()









    Out[21]:





{-I: 1, I: 1}



In [22]:

    
def eig_arr(funcs,vars,eq_pts):
    A=Jacobian_smp(funcs,vars)
    Anum=A.subs({vars[0]:eq_pts[0],vars[1]:eq_pts[1]})
    lam=Anum.eigenvals()
    return lam.keys()



In [27]:

    
def eig_arr2(funcs,vars,eq_pts):
    dim=len(funcs)
    A=Jacobian_smp2(funcs,vars)
    Anum=A.subs([(vars[i],eq_pts[i]) for i in np.arange(dim)])
    lam=Anum.eigenvals()
    return lam.keys()



In [24]:

    
eig_arr([f(x,y),g(x,y)],[x,y],[2,0])









    Out[24]:





[-1, -2]



In [28]:

    
eig_arr2([f(x,y),g(x,y)],[x,y],[2,0])









    Out[28]:





[-1, -2]



In [120]:

    
for eq_pt in eq_point:
    result=eig_arr([f(x,y),g(x,y)],[x,y],eq_pt)

    print " Eigenvalues of J at Point",eq_pt, ": ",result









    



 Eigenvalues of J at Point (0, 0) :  [2, 3]
 Eigenvalues of J at Point (0, 3) :  [-1, -3]
 Eigenvalues of J at Point (1, 1) :  [-1 + sqrt(2), -sqrt(2) - 1]
 Eigenvalues of J at Point (2, 0) :  [-1, -2]



In [29]:

    
for eq_pt in eq_point:
    result=eig_arr2([f(x,y),g(x,y)],[x,y],eq_pt)

    print " Eigenvalues of J at Point",eq_pt, ": ",result









    



 Eigenvalues of J at Point (0, 0) :  [2, 3]
 Eigenvalues of J at Point (0, 3) :  [-1, -3]
 Eigenvalues of J at Point (1, 1) :  [-1 + sqrt(2), -sqrt(2) - 1]
 Eigenvalues of J at Point (2, 0) :  [-1, -2]



In [ ]:

Exercise:

Consider the following ODE:

\begin{eqnarray} y_1' &=& a y_1\left(1-\frac{y_1}{K}\right)-by_1y_2 \\ y_2' &=& -r y_2+cy_1y_2-ey_2y_3 \\ y_3' &=& -f y_3+gy_2y_3 \end{eqnarray}

where the parameters are $a = 2, K = 1, b = 1, r = 1/2,c = 4, e = 2, f = 1,$ and $g = 2$.



In [ ]:



In [43]:

    
y1,y2,y3 = symbols("y1 y2 y3")
a=2;K=1;b=1;r=1/2.;c=4;e=2;f=1;g=2;

eq1=Eq(a*y1*(1-y1/K)-b*y1*y2)
eq2=Eq(-r*y2+c*y1*y2-e*y2*y3)
eq3=Eq(-f*y3+g*y2*y3)



In [31]:

    
eq3_point=solve([eq1,eq2,eq3],[y1,y2,y3]);print(eq3_point)









    



[(0.0, 0.0, 0.0), (0.0, 0.500000000000000, -0.250000000000000), (0.125000000000000, 1.75000000000000, 0.0), (0.750000000000000, 0.500000000000000, 1.25000000000000), (1.00000000000000, 0.0, 0.0)]



In [32]:

    
for eq_pt in eq3_point:
    result=eig_arr2([eq1,eq2,eq3],[y1,y2,y3],eq_pt)

    print " Eigenvalues of J at Point",eq_pt, ": ",result









    



 Eigenvalues of J at Point (0.0, 0.0, 0.0) :  [1, -1, 0]
 Eigenvalues of J at Point (0.0, 0.500000000000000, -0.250000000000000) :  [0]
 Eigenvalues of J at Point (0.125000000000000, 1.75000000000000, 0.0) :  [0]
 Eigenvalues of J at Point (0.750000000000000, 0.500000000000000, 1.25000000000000) :  [0]
 Eigenvalues of J at Point (1.00000000000000, 0.0, 0.0) :  [1/2 + sqrt(3)*I/2, -1, 1/2 - sqrt(3)*I/2]



In [ ]:



In [1]:

    
#from numpy import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline



In [2]:

    
from ivisual import *



In [3]:

    
a=2;K=1;b=1;r=1/2.;c=4;e=2;f=1;g=2;

def F(x,y,z): return a*x*(1-x/K)-b*x*y
def G(x,y,z): return -r*y+c*x*y-e*y*x
def H(x,y,z): return -f*z+g*y*z



In [4]:

    
scene = idisplay('3D Vector Field')

scene.autoscale = 0  # disable function of auto-zooming
scene.title = 'Vector Field'
scene.range = (4,4,4)

scene.background=color.white

xmin=-2;  xmax=2.;  ymin=-2.; ymax=2.;  zmin=-2.; zmax=2.
ngrid=21
xmesh=np.linspace(xmin,xmax,ngrid); ymesh=np.linspace(ymin,ymax,ngrid);zmesh=np.linspace(zmin,zmax,ngrid)
dh=(xmax-xmin)/(ngrid-1)

for xp in xmesh:
  for yp in ymesh:
    for zp in zmesh:    
        vx=F(xp,yp,zp); vy=G(xp,yp,zp); vz=H(xp,yp,zp) # direction vector v=(vx,vy)
        v=1./np.sqrt(1+vx*vx+vy*vy+vz*vz)*dh
        arrow(pos=(xp,yp,zp),axis=(v*vx,v*vy,v*vz),color=color.red,opacity=0.3)#,shaftwidth=0.1, fixedwidth=1)



In [7]:









    Out[7]:





<ivisual.visual.vector at 0x106a92390>



In [ ]:

Kermack-McKendrick Model: Strian SIR Model

It is used to explain the rapid rise and fall in the number of infective patients observed in epidemics

\begin{eqnarray*} \frac {d S}{d t} &=& -\beta S I\\ \frac {d I}{d t} &=& \beta S I -\gamma I\\ \frac {d R}{d t} &=& \gamma I \end{eqnarray*}

where $S, I$ and $R$ are the number of susceptible, infectious and recovered/immunized individ- uals respectively. $\beta$ is the transmission rate, $\gamma$ is the recovery rate.

Let $N$ denote the population size. Then $N$ is fixed since $$ \frac {d N}{d t} = \frac {d (S+I+R)}{d t}=0\to N\text{ is constant.} $$

\begin{eqnarray*} \frac {d S}{d t} &=& -\beta S I\\ \frac {d I}{d t} &=& (\beta S -\gamma) I\\ \Rightarrow ( S,I ) & = & \left\{\begin{array}{ll} ( 0,0 ) : & S- \text{nullclines}\\ \left( \frac{\gamma}{\beta} ,0 \right) : & I- \text{nullclines} \end{array}\right. \end{eqnarray*}

If $S(0) = S_0 < \frac{\gamma}{\beta}$, both $S(t)$ and $I(t)$ decreases and converges to a point on the $S$-axis.
f $S(0) = S_0 > \frac{\gamma}{\beta}$, $I(t)$ first increases in the region, $\left(\frac{\gamma}{\beta},1\right)$ and then decreases to 0, In this case, an outbreak occurs.

Define the basic reproduction number (epidemiological threshold) of this model: $$ R_0=\frac{N\gamma}{\beta}\sim\frac{S_0\gamma}{\beta}$$

Then \begin{eqnarray*} S_0 < \frac{\gamma}{\beta} &=& R_0<1\\ S_0 > \frac{\gamma}{\beta} &=& R_0>1\\ \end{eqnarray*}

$R_0<1$: each person who contracts to the disease will infect less than one person before dying or recovering;
$R_0>1$: the opposite occurs and there will be a outbreak of disease.



In [2]:

    
r=1.
beta=2.
S10,I10=0.2,0.8
S20,I20=0.8,0.2



In [3]:

    
def SIRmodel(gamma,beta,S0,I0,T=20,nt=1000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I=np.zeros(nt+1),np.zeros(nt+1)
    S[0]=S0;I[0]=I0
    for i in range(len(t)-1):
        S[i+1] = S[i]-dt*beta*S[i]*I[i]
        I[i+1] = I[i]+dt*(beta*S[i]*I[i]-gamma*I[i]) 
    return S,I



In [6]:

    
S1,I1=SIRmodel(r,beta,S10,I10)
S2,I2=SIRmodel(r,beta,S20,I20)



In [7]:

    
xy0=np.linspace(0,1,20)



In [24]:

    
plt.figure(figsize=[6,6])
plt.xlim(-0.1,1.1)
plt.ylim(-0.1,1.1)
plt.plot(S1,I1,'r',S2,I2,'b',xy0,1-xy0,'k')









    Out[24]:





[<matplotlib.lines.Line2D at 0x1066e8490>,
 <matplotlib.lines.Line2D at 0x1067fb090>,
 <matplotlib.lines.Line2D at 0x1067fb750>]



In [8]:

    
def SIRVecField(X,Y,xmin=-0.1,xmax=1.1,ymin=-0.1,ymax=1.1):
    xmesh=np.linspace(xmin,xmax,16); ymesh=np.linspace(ymin,ymax,16)
    for xp in xmesh:
        for yp in ymesh:
            vx=X(xp,yp); vy=Y(xp,yp) # direction vector v=(vx,vy)
            plt.quiver(xp, yp, vx, vy, pivot='middle', headwidth=4, headlength=4,color="red")



In [9]:

    
def f(x,y): return -beta*x*y
def g(x,y): return beta*x*y-r*y



In [12]:

    
plt.figure(figsize=[8,8])
plt.xlim(-0.1,1.1)
plt.ylim(-0.1,1.1)
plt.plot(S1,I1,'b',S2,I2,'b',xy0,1-xy0,'k')

plt.plot(r/beta*np.ones(len(xy0)),xy0,'k')
plt.plot(xy0,0*xy0,'k',0*xy0,xy0,'k',xy0,np.ones(len(xy0)),'k',np.ones(len(xy0)),xy0,'k')

plt.scatter((0.8,0.2),(0.2,0.8),c='b',marker='o',alpha=0.6)
plt.xlabel('S')
plt.ylabel('I')
plt.title('Strain S-I-R Epidemic-model',size=16)
SIRVecField(f,g)



In [ ]:

S-I-R Model with Birth

$N=1$
$\mu$: birth rate
Governed syetem equations: \begin{eqnarray*} \frac {d S}{d t} &=& \mu -\mu S -\beta S I\\ \frac {d I}{d t} &=& \beta S I -(\mu+\gamma) I\\ \frac {d R}{d t} &=& \gamma I+\mu R \end{eqnarray*}

\begin{eqnarray*} \frac {d S}{d t} &=& \mu-\mu S-\beta S I\\ \frac {d I}{d t} &=& (\beta S -\mu -\gamma) I\\ \Rightarrow ( S,I ) & = & \left\{\begin{array}{ll} \left( S, \frac{\mu}{\beta}\left (\frac{1}{S}-1 \right) \right) : & S- \text{nullclines}\\ \left( \frac{\mu+\gamma}{\beta} ,0 \right) : & I- \text{nullclines} \end{array}\right. \end{eqnarray*}

This results a triangle with vertices (0, 0), (1, 0) and (0, 1) on the $S I$-plane.

$P_1=(S_1^*,I_1^*)=(1,0)$ is always a equilibrium point.
- $S$- nullclines: $I=\left( S, \frac{\mu}{\beta}\left (\frac{1}{S}-1 \right) \right)$
- $I$- nullclines: $I=0$
another equilibrium point, $P_2=(S_2^*,I_2^*)$, is solution of \begin{eqnarray*} S\text{- nullclines: }I &=& \frac{\mu}{\beta}\left (\frac{1}{S}-1 \right)\\ I\text{- nullclines: }S &=& \frac{\mu+\gamma}{\beta} \\ \Rightarrow (S_2^*,I_2^*) &=& (\frac{\mu+\gamma}{\beta} , \frac{\mu}{\mu+\gamma}-\frac{\mu}{\beta}) \end{eqnarray*} exists if $\beta>\mu+\gamma$.

starting from $(S_1(0),I_1(0))$, $S(t)$ increases to $I(t)$ and decreases to equilibrium point, (no outbreak).
starting from $(S_2(0),I_2(0))$, I(t) increses to maximum where $\frac{\mu+\gamma}{\beta} $ and decreases to $I_2^*$, outbreak occurs.



In [32]:



In [1]:

    
from sympy import Symbol , dsolve , Function , Derivative , Eq, lambdify, symbols,Matrix,diff,solve,pprint
from sympy.abc import mu,gamma,beta
import sympy as sym



In [2]:

    
S,I,R = symbols("S I R")
DS=mu-mu*S-beta*S*I
DI=beta*S*I-(mu+gamma)*I
#FS=Function("FS")
#FI=Function("FI")
#FR=Function("FR")



In [11]:

    
EquiPnts=solve([DS,DI],[S,I]);
pprint(EquiPnts)









    



⎡        ⎛ ⎛     μ⋅(β - γ - μ)⎞                ⎞⎤
⎢        ⎜-⎜-μ + ─────────────⎟                ⎟⎥
⎢        ⎜ ⎝           β      ⎠   μ⋅(β - γ - μ)⎟⎥
⎢(1, 0), ⎜──────────────────────, ─────────────⎟⎥
⎣        ⎝          μ               β⋅(γ + μ)  ⎠⎦



In [15]:

    
EqMatrix=sym.Matrix([DS,DI]);pprint(EqMatrix)
JacMatrix=EqMatrix.jacobian([S,I]);pprint(JacMatrix)









    



⎡-I⋅S⋅β - S⋅μ + μ ⎤
⎢                 ⎥
⎣I⋅S⋅β - I⋅(γ + μ)⎦
⎡-I⋅β - μ     -S⋅β    ⎤
⎢                     ⎥
⎣  I⋅β     S⋅β - γ - μ⎦



In [15]:



In [23]:

    
def JacMatPnt(JacobianMatrix,Parameter,EquilibiumPnt):
    J=JacobianMatrix.subs([ (Parameter[0],EquilibiumPnt[0]), (Parameter[1],EquilibiumPnt[1] )])
    return J.eigenvals()



In [32]:

    
Eig1=JacMatPnt(JacMatrix,[S,I],EquiPnts[0]);pprint(Eig1)









    



⎧               __________                    __________   ⎫
⎪              ╱        2                    ╱        2    ⎪
⎨β   γ       ╲╱  (β - γ)       β   γ       ╲╱  (β - γ)     ⎬
⎪─ - ─ - μ - ─────────────: 1, ─ - ─ - μ + ─────────────: 1⎪
⎩2   2             2           2   2             2         ⎭



In [ ]:



In [48]:

    
Eig2=JacMatPnt(JacMatrix,[S,I],EquiPnts[1]);pprint(Eig2)









    



⎧                 ____________________________________________________________
⎪                ╱   ⎛ 2          2                  2      3       2         
⎨     β⋅μ      ╲╱  μ⋅⎝β ⋅μ - 4⋅β⋅γ  - 8⋅β⋅γ⋅μ - 4⋅β⋅μ  + 4⋅γ  + 12⋅γ ⋅μ + 12⋅γ
⎪- ───────── - ───────────────────────────────────────────────────────────────
⎩  2⋅(γ + μ)                                    2⋅γ + 2⋅μ                     

____________                      ____________________________________________
  2      3⎞                      ╱   ⎛ 2          2                  2      3 
⋅μ  + 4⋅μ ⎠           β⋅μ      ╲╱  μ⋅⎝β ⋅μ - 4⋅β⋅γ  - 8⋅β⋅γ⋅μ - 4⋅β⋅μ  + 4⋅γ  
────────────: 1, - ───────── + ───────────────────────────────────────────────
                   2⋅(γ + μ)                                    2⋅γ + 2⋅μ     

____________________________   ⎫
      2           2      3⎞    ⎪
+ 12⋅γ ⋅μ + 12⋅γ⋅μ  + 4⋅μ ⎠    ⎬
────────────────────────────: 1⎪
                               ⎭



In [47]:

    
sym.simplify(Eig2.keys()[0]*Eig2.keys()[1])









    Out[47]:





mu*(beta - gamma - mu)

Conclusion

If $\beta>\mu+\gamma$, the eigenvalues of Jacobian of systems at $P_2=(S_2^*,I_2^*)$ are all negative from above:

$\lambda_1<0$
$\lambda_1+\lambda_2=\mu(\beta-\gamma-\mu)>0$

Thus, if $P_2$ exists, it must be stable.



In [ ]:



In [ ]:



In [10]:

    
gamma=1.
beta=3.
mu=1.

S10,I10=0.1,0.9
S20,I20=0.9,0.1
S30,I30=1.,0.
xy0=np.linspace(0,1,20)



In [4]:

    
def f(x,y): return mu - mu*x -beta*x*y
def g(x,y): return -mu*y + beta*x*y-gamma*y



In [ ]:



In [5]:

    
def SIRBmodel(gamma,beta,mu,S0,I0,T=20,nt=1000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I=np.zeros(nt+1),np.zeros(nt+1)
    S[0]=S0;I[0]=I0
    for i in range(len(t)-1):
        S[i+1] = S[i]+dt*(mu-mu*S[i]-beta*S[i]*I[i])
        I[i+1] = I[i]+dt*(-mu*I[i]+beta*S[i]*I[i]-gamma*I[i]) 
    return S,I



In [6]:

    
def SIRVecField(X,Y,xmin=-0.1,xmax=1.1,ymin=-0.1,ymax=1.1):
    xmesh=np.linspace(xmin,xmax,16); ymesh=np.linspace(ymin,ymax,16)
    for xp in xmesh:
        for yp in ymesh:
            vx=X(xp,yp); vy=Y(xp,yp) # direction vector v=(vx,vy)
            plt.quiver(xp, yp, vx, vy, pivot='middle', headwidth=4, headlength=4,color="red")



In [8]:

    
S1,I1=SIRBmodel(gamma,beta,mu,S10,I10)
S2,I2=SIRBmodel(gamma,beta,mu,S20,I20)
S3,I3=SIRBmodel(gamma,beta,mu,S30,I30)



In [12]:

    
plt.figure(figsize=[8,8])
plt.xlim(-0.1,1.1)
plt.ylim(-0.1,1.1)
plt.plot(S1,I1,'b',S2,I2,'b',S3,I3,'brown',xy0,1-xy0,'k')

plt.plot((gamma+mu)/beta*np.ones(len(xy0)),xy0,'k')

plt.plot(xy0,0*xy0,'k',0*xy0,xy0,'k',xy0,np.ones(len(xy0)),'k',np.ones(len(xy0)),xy0,'k')

plt.scatter((0.9,0.1),(0.1,0.9),c='b',marker='o',alpha=0.6)
plt.xlabel('S')
plt.ylabel('I')
plt.title('Strain S-I-R Epidemic-model with Birth',size=16)
SIRVecField(f,g)



In [19]:

    
fig = plt.figure()
ax=Axes3D(fig)
ax.scatter3D(S1,I1,1.-S1-I1)
ax.scatter(S1,I1)

#ax.scatter3D(S2,I2,1.-S2-I2)
#ax.scatter3D(S3,I3,1.-S3-I3)









    Out[19]:





<mpl_toolkits.mplot3d.art3d.Patch3DCollection at 0x1072915d0>



In [ ]:

Seasonality Model



In [ ]:

$\beta(t)=\beta_0 (1+\kappa\cos(2\pi t))$, $0\le\kappa\le1$



In [52]:

    
r=100.
beta=1800.
mu=0.02
ka=0.07

S10,I10=0.7,0.3



In [53]:

    
def f(x,y,t): return mu - mu*x -betat(t)*x*y
def g(x,y,t): return -mu*y + betat(t)*x*y-r*y
def betat (t): return  beta*(1+ka*np.cos(2*np.pi*t))



In [54]:

    
def SIRBSmodel(gamma,betat,mu,S0,I0,T=100,nt=100000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I=np.zeros(nt+1),np.zeros(nt+1)
    S[0]=S0;I[0]=I0
    for i in range(len(t)-1):
        S[i+1] = S[i]+dt*(mu-mu*S[i]-betat(t[i])*S[i]*I[i])
        I[i+1] = I[i]+dt*(-mu*I[i]+betat(t[i])*S[i]*I[i]-gamma*I[i]) 
    return S,I



In [55]:

    
T=100
nt=100000
dt=100/float(nt)
t=np.linspace(0,T,nt+1)



In [56]:

    
S1,I1=SIRBSmodel(r,betat,mu,S10,I10)



In [90]:

    
plt.plot(t[-25000:],-np.log(I1[-25000:]))









    Out[90]:





[<matplotlib.lines.Line2D at 0x10b6d3910>]



In [92]:

    
plt.figure(figsize=[4,4])
plt.plot(-np.log(S1[-25000:]),-np.log(I1[-25000:]))









    Out[92]:





[<matplotlib.lines.Line2D at 0x10e7a4590>]



In [93]:

    
ka=0.08
S1,I1=SIRBSmodel(r,betat,mu,S10,I10)



In [94]:

    
plt.plot(t[-25000:],-np.log(I1[-25000:]))









    Out[94]:





[<matplotlib.lines.Line2D at 0x10ed505d0>]



In [96]:

    
plt.figure(figsize=[4,4])
plt.plot(-np.log(S1[-25000:]),-np.log(I1[-25000:]))
plt.xlabel('-log(S)')
plt.ylabel('-log(I)')









    Out[96]:





<matplotlib.text.Text at 0x10f025410>



In [114]:

    
ka=0.0865
S1,I1=SIRBSmodel(r,betat,mu,S10,I10)



In [115]:

    
plt.plot(t[-25000:],-np.log(I1[-25000:]))









    Out[115]:





[<matplotlib.lines.Line2D at 0x111db4850>]



In [116]:

    
plt.figure(figsize=[4,4])
plt.plot(-np.log(S1[-25000:]),-np.log(I1[-25000:]))
plt.xlabel('-log(S)')
plt.ylabel('-log(I)')









    Out[116]:





<matplotlib.text.Text at 0x111df6690>



In [100]:

    
ka=0.09
S1,I1=SIRBSmodel(r,betat,mu,S10,I10)



In [ ]:



In [101]:

    
plt.plot(t[-25000:],-np.log(I1[-25000:]))









    Out[101]:





[<matplotlib.lines.Line2D at 0x10f4e0f50>]



In [102]:

    
plt.figure(figsize=[4,4])
plt.plot(-np.log(S1[-25000:]),-np.log(I1[-25000:]))
plt.xlabel('-log(S)')
plt.ylabel('-log(I)')









    Out[102]:





<matplotlib.text.Text at 0x10f586d90>



In [ ]:



In [57]:

    
def f(x,y,t): return mu - mu*x -betat(t)*x*y
def g(x,y,t): return -mu*y + betat(t)*x*y-r*y
def betat (t): return  beta*(1+ka*np.cos(2*np.pi*t))



In [57]:



In [58]:

    
def SIRBSmodel(gamma,betat,mu,S0,I0,T=100,nt=100000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I=np.zeros(nt+1),np.zeros(nt+1)
    S[0]=S0;I[0]=I0
    
    for i in range(len(t)-1):
        S[i+1] = S[i]+dt*(mu-mu*S[i]-betat(t[i])*S[i]*I[i])
        I[i+1] = I[i]+dt*(-mu*I[i]+betat(t[i])*S[i]*I[i]-gamma*I[i]) 
    return S,I



In [59]:

    
ka=0.07
T=100
nt=100000

def vis(ka): 
    fig, ax = plt.subplots(figsize=(4,4))
    gamma=100.
    beta=1800.
    mu=0.02
    S10,I10=0.7,0.3
    
    def f(x,y,t): return mu - mu*x -betat(t)*x*y
    def g(x,y,t): return -mu*y + betat(t)*x*y-gamma*y
    def betat (t): return  beta*(1+ka*np.cos(2*np.pi*t))

    
    #S1,I1=SIRBSmodel(r,betat,mu,S10,I10)
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I=np.zeros(nt+1),np.zeros(nt+1)
    S[0]=S10;I[0]=I10
    
    for i in range(len(t)-1):
        S[i+1] = S[i]+dt*(mu-mu*S[i]-betat(t[i])*S[i]*I[i])
        I[i+1] = I[i]+dt*(-mu*I[i]+betat(t[i])*S[i]*I[i]-gamma*I[i]) 
    #ax.plot(t[-25000:],-np.log(I1[-25000:]))
    ax.plot(-np.log(S[-25000:]),-np.log(I[-25000:]))

    return fig



In [ ]:

    
i = interact(vis,
         ka=widgets.FloatSliderWidget(min=0.0, max=0.1, step=0.001, value=.07, description="ka"),
         )



In [ ]:

Two-Strain Influenza Model



In [ ]:

$$ \frac{d}{d t} \left(\begin{array}{c} S_{1}\\ I_{1}\\ R_{1}\\ S_{2}\\ I_{2}\\ R_{2} \end{array}\right) = \left(\begin{array}{cccccc} - \delta_{1} - \beta_{t} I_{1} & 0 & 0 & \delta_{2} & 0 & g_{2}\\ 0 & \beta_{t} S_{1} - \gamma & 0 & 0 & 0 & 0\\ 0 & \gamma & -g_{1} & 0 & 0 & 0\\ \delta_{1} & 0 & g_{1} & - \delta_{2} - \beta_{t} I_{2} & 0 & 0\\ 0 & \beta_{t} S_{_{2}} - \gamma & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \gamma & -g_{2} \end{array}\right) \left(\begin{array}{c} S_{1}\\ I_{1}\\ R_{1}\\ S_{2}\\ I_{2}\\ R_{2} \end{array}\right) $$

where

$(S_i,I_i,R_i), i=1,2,$ are the number of susceptible, infectious and recovered/immunized individuals with repective to strain $i=1,2$.
$\beta$ is the transmission rate,
$\gamma$ ($\gamma^{-1}=3$ days) is the recovery rate.
Parameters $\delta_1,\delta_2$ reflect the evolution and competition of these two strains.

Assumptions

$N=S_1+I_1+R_1+S_2+I_2+R_2=1$
$g_1=g_2=1$

Stability Analysis



In [16]:

    
from sympy import Symbol , dsolve , Function , Derivative , Eq, lambdify, symbols,Matrix,diff,solve,pprint
from sympy.abc import mu,gamma,beta
import sympy as sym



In [17]:

    
S1,I1,R1,S2,I2,R2 = symbols("S1 I1 R1 S2 I2 R2")
delta1,delta2=symbols("delta1 delta2")

DS1=gamma*I2-beta*S1*I1-delta1*S1+delta2*S2
DI1=(beta*S1-gamma)*I1
#DR1=gamma*I1-R1

DS2=gamma*I1-beta*S2*I2-delta2*S2+delta1*S1
DI2=(beta*S2-gamma)*I2
C1=I2
#DR2=gamma*I2-R2



In [3]:

    
solve([I1,DS1,DI1,DS2,DI2],[S1,I1,S2,I2])









    Out[3]:





[]

Sympy.solve can not solve the two strain SIR system. However, it could be solved manually by reducing the system of equations

Case

$(\beta*S_1-\gamma)*I_1=0$, $(\beta*S_2-\gamma)*I_2=0$:

a) $I_1,I_2=0$

b) $\beta*S_1-\gamma=0,I_2=0\Rightarrow S_1=\gamma/\beta,I_2=0$

c) $\beta*S_2-\gamma=0,I_1=0\Rightarrow S_2=\gamma/\beta,I_1=0$

d) $\beta*S_1-\gamma=0,I_2=0,\beta*S_2-\gamma=0\Rightarrow S_1,S_2=\gamma/\beta$



In [21]:

    
# case a)

DS11=-delta1*S1+delta2*S2
DS21=-delta2*S2+delta1*S1
sum= S1+S2-1



In [23]:

    
EquiPnt1=solve([DS11,DS21,sum],[S1,S2]);pprint(EquiPnt1)









    



⎧       δ₂           δ₁  ⎫
⎨S₁: ───────, S₂: ───────⎬
⎩    δ₁ + δ₂      δ₁ + δ₂⎭



In [38]:

    
# case b)
# I2,R2=0, S1=gamma/beta, R1=gamma*I1

DS12=-gamma*I1-delta1*gamma/beta+delta2*S2
DI12=(gamma-gamma)*I1
DS22=gamma*I1-delta2*S2+delta1*gamma/beta
#cond12=S1-gamma/beta
cond22=gamma/beta+I1+gamma*I1+S2-1



In [39]:

    
EquiPnt2=solve([DS12,DS22,cond22],[I1,S2]);pprint(EquiPnt2)









    



⎧    β⋅δ₂ - δ₁⋅γ - δ₂⋅γ      γ⋅(β + δ₁⋅(γ + 1) - γ)⎫
⎨I₁: ──────────────────, S₂: ──────────────────────⎬
⎩    β⋅(δ₂⋅γ + δ₂ + γ)         β⋅(δ₂⋅(γ + 1) + γ)  ⎭



In [ ]:



In [42]:

    
# case c)
# I1,R1=0, S2=gamma/beta, R2=gamma*I2

DS13=-gamma*I2+delta1*S1-delta2*gamma/beta
DS23=gamma*I2-delta1*S1+delta2*gamma/beta

#cond12=S1-gamma/beta
cond23=gamma/beta+I2+gamma*I2+S1-1



In [43]:

    
EquiPnt3=solve([DS13,DS23,cond23],[I2,S1]);pprint(EquiPnt3)gamma/beta+R1/gamma+R1









    



⎧    β⋅δ₁ - δ₁⋅γ - δ₂⋅γ      γ⋅(β + δ₂⋅(γ + 1) - γ)⎫
⎨I₂: ──────────────────, S₁: ──────────────────────⎬
⎩    β⋅(δ₁⋅γ + δ₁ + γ)         β⋅(δ₁⋅(γ + 1) + γ)  ⎭



In [44]:

    
# case d)
# S1=S2=gamma/beta, R1=gamma*I1, R2=gamma*I2

DS14=R2-R1-delta1*gamma/beta+delta2*gamma/beta
DS24=R1-R2-delta2*gamma/beta+delta1*gamma/beta
cond24=gamma/beta+R1/gamma+R1+gamma/beta+R2/gamma+R2-1



In [45]:

    
EquiPnt4=solve([DS14,DS24,cond24],[R1,R2]);pprint(EquiPnt4)









    



⎧    γ⋅(β - δ₁⋅γ - δ₁ + δ₂⋅γ + δ₂ - 2⋅γ)      γ⋅(β - 2⋅γ + (δ₁ - δ₂)⋅(γ + 1))⎫
⎨R₁: ───────────────────────────────────, R₂: ───────────────────────────────⎬
⎩                2⋅β⋅(γ + 1)                            2⋅β⋅(γ + 1)          ⎭



In [ ]:



In [27]:

    
def betat (t,ka): return  beta*(1+ka*np.cos(2*np.pi*t))
def SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20,T=50,nt=10000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I,R=np.zeros([2,nt+1]),np.zeros([2,nt+1]),np.zeros([2,nt+1])
    S[0,0]=S10;I[0,0]=I10;R[0,0]=R10
    S[1,0]=S20;I[1,0]=I20;R[1,0]=R20
    
    for i in range(len(t)-1):
        S[0,i+1] = S[0,i]+dt*(R[1,i]-betat(t[i],ka)*S[0,i]*I[0,i]-delta1*S[0,i]+delta2*S[1,i])
        S[1,i+1] = S[1,i]+dt*(R[0,i]-betat(t[i],ka)*S[1,i]*I[1,i]-delta2*S[1,i]+delta1*S[0,i])
        I[0,i+1] = I[0,i]+dt*(betat(t[i],ka)*S[0,i]*I[0,i]-gamma*I[0,i]) 
        I[1,i+1] = I[1,i]+dt*(betat(t[i],ka)*S[1,i]*I[1,i]-gamma*I[1,i]) 
        R[0,i+1] = R[0,i]+dt*(gamma*I[0,i]-R[0,i])
        R[1,i+1] = R[1,i]+dt*(gamma*I[1,i]-R[1,i])
    return S,I,R,t



In [20]:

    
delta1=0.9706
delta2=1.8996
beta=33.6607
S10 = .2; I10 =0.01; R10 = 0.4; S20 = .2; I20 = 1./30. ; R20 = 1-S10-I10-R10-S20-I20 ;

gamma=3.4675



In [29]:

    
# ka=0,
ka=0.
S10 = .2; I10 =0.01; R10 = 0.4; S20 = .2; I20 = 1./30. ; R20 = 1-S10-I10-R10-S20-I20 ;

delta1=0.7801
delta2=0.9553
beta=27.7309

gamma=17.5047



In [30]:

    
S,I,R,t=SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20)



In [52]:

    
fig = plt.figure(figsize=[12,4])

ax = fig.add_subplot(131)
ax.set_ylim(0.2,0.46)
ax.plot(S[0,:],S[1,:])
ax = fig.add_subplot(132)
ax.plot(I[0,:],I[1,:])
ax = fig.add_subplot(133)
ax.plot(R[0,:],R[1,:])









    Out[52]:





[<matplotlib.lines.Line2D at 0x10c22a090>]



In [52]:



In [ ]:



In [53]:

    
# ka=0,
ka=0.
S10 = .2; I10 =0.01; R10 = 0.4; S20 = .2; I20 = 1./30. ; R20 = 1-S10-I10-R10-S20-I20 ;

delta1=0.5402
delta2=0.8004
beta=18.4108

gamma=9.7299



In [54]:

    
S,I,R,t=SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20)



In [55]:

    
fig = plt.figure(figsize=[12,4])

ax = fig.add_subplot(131)
ax.set_ylim(0.2,0.46)
ax.plot(S[0,:],S[1,:])
ax = fig.add_subplot(132)
ax.plot(I[0,:],I[1,:])
ax = fig.add_subplot(133)
ax.plot(R[0,:],R[1,:])









    Out[55]:





[<matplotlib.lines.Line2D at 0x10c6efad0>]



In [60]:

    
fig = plt.figure()
ax=Axes3D(fig)
ax.scatter3D(S[0,:],I[0,:],R[0,:])









    Out[60]:





<mpl_toolkits.mplot3d.art3d.Patch3DCollection at 0x10cd36390>



In [59]:

    
fig = plt.figure()

ax1=Axes3D(fig)
ax1.scatter3D(S[1,:],I[1,:],R[1,:])









    Out[59]:





<mpl_toolkits.mplot3d.art3d.Patch3DCollection at 0x102066410>



In [ ]:



In [ ]:



In [21]:

    
def SIRmodel(gamma,beta,delta1,delta2,S10,S20,I10,I20,R10,R20,T=50,nt=1000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I,R=np.zeros([2,nt+1]),np.zeros([2,nt+1]),np.zeros([2,nt+1])
    S[0,0]=S10;I[0,0]=I10;R[0,0]=R10
    S[1,0]=S20;I[1,0]=I20;R[1,0]=R20
    
    for i in range(len(t)-1):
        S[0,i+1] = S[0,i]+dt*(R[1,i]-beta*S[0,i]*I[0,i]-delta1*S[0,i]+delta2*S[1,i])
        S[1,i+1] = S[1,i]+dt*(R[0,i]-beta*S[1,i]*I[1,i]-delta2*S[1,i]+delta1*S[0,i])
        I[0,i+1] = I[0,i]+dt*(beta*S[0,i]*I[0,i]-gamma*I[0,i]) 
        I[1,i+1] = I[1,i]+dt*(beta*S[1,i]*I[1,i]-gamma*I[1,i]) 
        R[0,i+1] = R[0,i]+dt*(gamma*I[0,i]-R[0,i])
        R[1,i+1] = R[1,i]+dt*(gamma*I[1,i]-R[1,i])
    return S,I,R,t



In [22]:

    
S,I,R,t=SIRmodel(gamma,beta,delta1,delta2,S10,S20,I10,I20,R10,R20)



In [23]:

    
plt.ylim(0,0.3)
plt.plot(t,I[0,:],'b',t,I[1,:],'r')









    Out[23]:





[<matplotlib.lines.Line2D at 0x107455c50>,
 <matplotlib.lines.Line2D at 0x102059950>]



In [ ]:

Seasoning Model

$\beta(t)=\beta_0 (1+\kappa\cos(2\pi t))$, $0\le\kappa\le1$



In [80]:

    
ka=0.07

gamma=4.6681; delta1 =0.0347; delta2 =0.0538; beta0 =11.5110;
S10 = .2; I10 =0.01; R10 = 0.4; S20 = .2; I20 = 1./30. ; R20 = 1-S10-I10-R10-S20-I20 ;



In [45]:

    
def betat (t,ka): return  beta*(1+ka*np.cos(2*np.pi*t))
def SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20,T=50,nt=1000):
    t=np.linspace(0,T,nt+1)
    dt=t[1]-t[0]
    S,I,R=np.zeros([2,nt+1]),np.zeros([2,nt+1]),np.zeros([2,nt+1])
    S[0,0]=S10;I[0,0]=I10;R[0,0]=R10
    S[1,0]=S20;I[1,0]=I20;R[1,0]=R20
    
    for i in range(len(t)-1):
        S[0,i+1] = S[0,i]+dt*(R[1,i]-betat(t[i],ka)*S[0,i]*I[0,i]-delta1*S[0,i]+delta2*S[1,i])
        S[1,i+1] = S[1,i]+dt*(R[0,i]-betat(t[i],ka)*S[1,i]*I[1,i]-delta2*S[1,i]+delta1*S[0,i])
        I[0,i+1] = I[0,i]+dt*(betat(t[i],ka)*S[0,i]*I[0,i]-gamma*I[0,i]) 
        I[1,i+1] = I[1,i]+dt*(betat(t[i],ka)*S[1,i]*I[1,i]-gamma*I[1,i]) 
        R[0,i+1] = R[0,i]+dt*(gamma*I[0,i]-R[0,i])
        R[1,i+1] = R[1,i]+dt*(gamma*I[1,i]-R[1,i])
    return S,I,R,t



In [40]:

    
gamma=3.4675
delta1=0.9706
delta2=1.8996
beta0=33.6607
S10=0.4
I10=0.1
S20=0.4
I20=0.1
R10=0.
R20=0.



In [68]:

    
S10 = .2; I10 = 10/1e3; R10 = 2*S10; 
S20 = .2; I20 = 10/(3e2) ; R20 = 1-S10-I10-R10-S20-I20 ;
gamma=1.5499; delta1 =0.8473; delta2 =1.4234; beta0 =4.2359
ka=0.4
S,I,R,t=SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20)



In [70]:

    
plt.ylim(0,.3)
plt.plot(t,I[0,:],'b',t,I[1,:],'r')









    Out[70]:





[<matplotlib.lines.Line2D at 0x109005310>,
 <matplotlib.lines.Line2D at 0x1092617d0>]



In [73]:

    
S10 = .2; I10 = 10/1e3; R10 = 2*S10; 
S20 = .2; I20 = 10/(3e2) ; R20 = 1-S10-I10-R10-S20-I20 ;
gamma=4.6681; delta1 =0.0347; delta2 =0.0538; beta0 =11.5110;

ka=0.1
S,I,R,t=SIRmodelP(gamma,betat,delta1,delta2,ka,S10,S20,I10,I20,R10,R20)



In [74]:

    
plt.ylim(0,0.5)
plt.plot(t,I[0,:],'b',t,I[1,:],'r')









    Out[74]:





[<matplotlib.lines.Line2D at 0x109494050>,
 <matplotlib.lines.Line2D at 0x109574f50>]



In [84]:

    
gamma=1.5499; delta1 =0.8473; delta2 =1.4234; beta0 =4.2359;



In [87]:

    
ka=0.3
S,I,R=SIRmodelP(gamma,betat,ka,delta1,delta2,S10,S20,I10,I20,R10,R20)



In [88]:

    
t=np.linspace(0,20,1001)
plt.ylim(0,1)
plt.plot(t,I[0,:],t,I[1,:])









    Out[88]:





[<matplotlib.lines.Line2D at 0x107aa6b50>,
 <matplotlib.lines.Line2D at 0x107f94410>]



In [ ]: