SIR Model

The example describes how, given a contagious illness, the population varies among the three groups of the indeividuals that are susceptible to the illness ($S$), infected by the illness ($I$) and recovered from the illness ($R$).

Nonlinear Ordinary Differential Equations for $$ S\overset{\beta I}{\longrightarrow}I\overset{\gamma}{\longrightarrow}R$$

\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*}



In [1]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy import log
from scipy import integrate
from scipy.optimize import fsolve

import scipy.linalg as la
import scipy.integrate as spi

from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
#from matplotlib import rc
#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
#rc('font',**{'family':'serif','serif':['Palatino']})
#rc('text', usetex=True)



In [2]:

    
S0=0.995;I0=.005;R0=0.
beta=0.;gamma=0.35
n=1000
T=20
dt=T/float(n)
t=np.linspace(0,T,n+1)



In [3]:

    
plt.xlim(0,20)
plt.ylim(0,1.1)

beta=0.75
S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;I[0]=I0;R[0]=R0
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]) 
    #R[i+1] = 1.-S[i+1]-I[i+1]
R=1.-S-I
plt.plot(t,S,t,I,t,R) 
plt.legend(['S','I','R'])
plt.grid(True)
plt.title(r'$\beta=0.75,\gamma=0.35$',size=20)









    Out[3]:





<matplotlib.text.Text at 0x109166050>

Solve system of ODE's by Scipy

Procedures

define time partition, t_range, initials, INPUT:(S0,I0);
define Governed equation,diff_eqs;
solved by scipy.integrate.odeint(diff_eqs, INPUT, t_range)



In [4]:

    
plt.xlim(0,20)
plt.ylim(0,1.1)

beta=0.75
S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;I[0]=I0;R[0]=R0
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]) 
    #R[i+1] = 1.-S[i+1]-I[i+1]
R=1.-S-I

def diff_eqs(INP,t):  
    '''The main set of equations'''
    Y=np.zeros((2))
    V = INP    
    Y[0] = - beta * V[0] * V[1]
    Y[1] = beta * V[0] * V[1] - gamma * V[1] 
    return Y   # For odeint



t_start = 0.0; t_end = T; t_inc = dt
t_range = np.arange(t_start, t_end+t_inc, t_inc)
INPUT=(S0,I0)
RES = spi.odeint(diff_eqs,INPUT,t_range)

Rec=1. - (RES[:,0]+RES[:,1])


plt.plot(t,RES[:,0],t,RES[:,1],t,Rec) 
plt.legend(['S','I','R'])
plt.grid(True)
plt.title(r'$\beta=0.75,\gamma=0.35$',size=20)









    Out[4]:





<matplotlib.text.Text at 0x1092b33d0>



In [16]:



In [5]:

    
beta=2.
def vis2(beta,gamma):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    S0=0.995;I0=.005;R0=0.
    n=1000
    T=20
    dt=T/float(n)
    t=np.linspace(0,T,n+1)
    
    ax.set_xlim(0,T)
    ax.set_ylim(0,1.1)

    S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
    S[0]=S0;I[0]=I0;R[0]=R0
    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]) 
        R[i+1] = 1.-S[i+1]-I[i+1]
    ax.plot(t,S)
    ax.plot(t,I)
    ax.plot(t,R)
    ax.set_title("SIR model")

    ax.set_xlabel("Time (days)")
    ax.set_ylabel("Percent of Population")
    ax.legend(('S','I','R'))
    ax.grid(True)
    return fig



In [6]:

    
i = interact(vis2,
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=5, step=0.1, value=1.5, description="gamma"),
         )



In [ ]:



In [24]:

    
beta=2.
gamma=1.5
S0=0.95
I0=0.05
R0=0.
def vis2(beta,gamma,S0,I0):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    n=1000
    T=40
    R0=1-S0-I0
    dt=T/float(n)
    
    t=np.linspace(0,T,n+1)
    R0=gamma/beta
    ax.set_xlim(0,T)
    ax.set_ylim(0,1.1)

    S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
    S[0]=S0;I[0]=I0;R[0]=R0
    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]) 
        R[i+1] = 1.-S[i+1]-I[i+1]
    ax.plot(t,S)
    ax.plot(t,I)
    ax.plot(t,R)
    ax.set_title('SIR model $R_0=$%f' %R0)

    ax.set_xlabel("Time (days)")
    ax.set_ylabel("Percent of Population")
    ax.legend(('S','I','R'))
    ax.grid(True)
    return fig



In [25]:

    
i = interact(vis2,
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=2.0, step=0.1, value=1.5, description="gamma"),
         S0=widgets.FloatSliderWidget(min=0.9, max=1.0, step=0.01, value=0.95, description="S0"),
         I0=widgets.FloatSliderWidget(min=0., max=0.1, step=0.01, value=0.01, description="I0"),
         )



In [29]:

Reproduction Ratio

$\frac{d I}{d t}$ : rate of changes (growth velocity) of infectives.
The number of infectives grows as long as the proportion of susceptibles, $S(t)$, in the population exceeds $\frac{\gamma}{\beta}$, i.e. $$ \frac{d I}{d t}=\beta S I -\gamma I\ge(\beta\frac{\gamma}{\beta}-\gamma)\cdot I\ge0$$

Defintion

$R_0=\frac{\gamma}{\beta}$ is called repoduction ratio, and

$R_0>1$: epidemic is out of control.
$R_0<1$: epidemic is under control.

What relation between $S$ and $I$

Theoretically, it's hard to solve from $S-I-R$ ODE's.

Fortunately, we can reduce the $S-I$ DE's into: $$ \frac { d I}{d S}=\frac{\frac{d I}{d S}}{\frac{d S}{d I}}=\frac{\beta I S-\gamma I}{-\beta I S} =-1+\frac{\gamma}{\beta S}$$ And $$ I=\int\left(\frac{\gamma}{\beta S} -1\right) d S=\frac{\gamma}{\beta} \ln S -S+C$$ Replcing the initial values, $S(0),I(0)$, gets the constant $C$: $$ C=I(0)-\frac{\gamma}{\beta}\ln(S(0))+S(0)$$



In [2]:



In [11]:

    
beta=2.
gamma=1.5
S0=0.3
I0=0.05
R0=0.
def SI3(beta,gamma):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    n=1000
    St=np.linspace(S0,1,n+1)
    
    ax.set_xlim(0,1)
    ax.set_ylim(0,1.1)
    I=gamma*log(St)/beta-St+I0+S0-gamma*log(S0)/beta

    ax.plot(St,I)
    ax.set_title("S-I Relation")

    ax.set_xlabel("S")
    ax.set_ylabel("I")

    ax.grid(True)
    return fig



In [12]:

    
i = interact(SI3,
         beta=widgets.FloatSliderWidget(min=0.0, max=10.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=10., step=0.1, value=1.5, description="gamma"),
         )

Note

Facts:

$I=0\to S,I,R$ are contstants: $I$ is an equivalent state, i.e. nothing will happen.
If case starts out at $I(0)>0$, $\frac{ d S}{d t}<0$ implies $S$ being decreasing.
$ I=\frac{\gamma}{\beta}\ln S -S +C \to -\infty$ as $S\to0$.

Conclusion

Epidemic will end as $I\to 0$ while $S$ approaches some positive term.
While $I$ is at equilibium, $S(t)=S_T$ for some $t\ge t_S$

Questions

Observe the relation between $S$ and $I$ by making plots of $S-I$

Problem: How long is $t_S$?

To get more precisely, i.e. when would it be in control, duration time for approaching $S(t)$ has to to be estimested. Replacing $I=I(S)$ above gets:

\begin{eqnarray} \frac{d S}{d t} &= &-\beta I S \cr &= &-\gamma S \ln S +\beta S^2-\beta C S \end{eqnarray}

This implies $$\int^{S_T}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S} =\int^T_0 d t =T$$

Solution Processes

It's true: epidemic will explose and then enventually be under control.
It is more interesting to know when $I$ come back to $I_0$ again than expect to get $I=0$.

As we note above, when to get $S=S_1$ while $I=I_0$ come back again?

Answer

This is the key point:

\begin{eqnarray} C&=&I(t)-\frac{\gamma}{\beta}\ln S(t)+S(t) \\ &=& I(0)-\frac{\gamma}{\beta}\ln(S(0))+S(0) \\ &=& I(0)-\frac{\gamma}{\beta}\ln S_1 +S_1 \end{eqnarray}

Replacing the upper limit in the integral and get the time, $T$, when the epidemic will be under control: \begin{eqnarray} T&=&\int^T_0 d t \\ &=&\int^{S_1}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S}\qquad(\text{ where } S_1<S_0) \end{eqnarray}



In [13]:



In [2]:

    
b=8.
g=2.3
S0=0.99
I0=0.01
R0=0.
C0= I0-g/b*log(S0)+S0



In [3]:

    
print C0









    



1.00288947156



In [4]:

    
def func(x):
        return -g/b*np.log(x)+x+I0-C0
SS=fsolve(func,0.01);print(SS)









    



[ 0.03583275]



In [18]:

    
def xv0(x):
    return 0*x
def intf(x):
    return 1/(b*x*x-b*C0*x-g*x*log(x))



In [28]:

    
xv=np.linspace(SS,S0,100)
xv2=np.linspace(0,1,101)

plt.ylim(-50,5)
plt.text(0.2,-20,r'$\int^{S_1}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S}$',size=20)
plt.text(0.2,-30,r'$S_0=0.99,S_1=0.0358$',size=12)

plt.plot(xv2,xv0(xv2))
plt.plot(xv,intf(xv),'b--')
plt.fill_between(xv, intf(xv),0, facecolor='red', alpha=0.1)









    Out[28]:





<matplotlib.collections.PolyCollection at 0x10b10b690>



In [11]:

    
x2 = lambda v: 1/(b*v*v-b*C0*v-g*v*log(v))
integrate.quad(x2,S0,SS)









    Out[11]:





(3.0758505917569066, 1.4204760932643643e-09)

Conclusion

After $3.0758505917569066$ units of time, it would be under control.

Computer Practice

Make an interactive demo with input, $\beta,\gamma,R_0,I_0$, and $T$ output.



In [ ]:

Zombie Apocalypse

Philip Munz, Ioan Hudea, Joe Imad and Robert J. Smith 2009

\begin{eqnarray*} \frac {d S}{d t} &=& -\beta S Z\\ \frac {d Z}{d t} &=& (\beta-\alpha) S Z +\zeta R\\ \frac {d R}{d t} &=& \alpha S Z -\zeta R\ \end{eqnarray*}

Encounters between zombies and live humans can result in the human becoming a zombie, $\beta S Z$ ;
the zombie becoming dead, $\alpha S Z$;
The dead sometimes spontaneously turn into zombies, $\zeta R$.



In [18]:

    
beta=1.
alpha=4.
zeta=0.2
S0=.99
Z0=0.01
R0=0.0
T=100.
n=10000
t=np.linspace(0,T,n+1)
dt=t[1]-t[0]
S,Z,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;Z[0]=Z0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*Z[i]
    Z[i+1] = Z[i]+dt*((beta-alpha)*S[i]*Z[i]+zeta*R[i]) 
    R[i+1] = 1.-S[i+1]-Z[i+1]
plt.plot(t,S,t,Z,t,R)
plt.legend(('S','Z','R'))
print('S[100]','Z[100]','R[100]:',S[-1],Z[-1],R[-1])









    



('S[100]', 'Z[100]', 'R[100]:', 5.1107327881136156e-11, 0.99341788079421345, 0.0065821191546792113)



In [7]:

    
t[:5]









    Out[7]:





array([ 0. ,  0.1,  0.2,  0.3,  0.4])



In [ ]:

Ebola: SEIR Model

Ref:

G.Chowell, N.W.Hengartner,C.Castillo-Chavez,P.W.Fenimore,J.M.Hyman, The basic reproductive number of Ebola and the effects of public health measures: the cases of Congo and Uganda, Journal of Theoretical Biology, 2004
Anderson, R.M., May, R.M., 1991. Infectious Diseases of Humans. Oxford University Press, Oxford.
Brauer, F., Castillo-Chavez, C., 2000. Mathematical Models in Sartwell, P.E., 1966. The distribution of incubation periods and the PopulationBiologyandEpidemiology.Springer,NewYork.
Matt J. Keeling & Pejman Rohani, Modeling Infectious Diseases in Humans and Animals

Susceptible individuals in class $S$ in contact with the virus enter the exposed class $E$ at the per-capita rate $\beta I/N$; where $\beta$ is transmission rate per person per day, $N$is the total effective population size, and $I/N$ is the probability that a contact is made with a infectious individual (i.e. uniform mixing is assumed). Exposed individuals undergo an average incubation period (assumed asymptomatic and uninfectious) of $1/k$ days before progressing to the infectious class $I$: Infectious individuals move to the $R$-class (death or recovered) at the per-capita rate $\gamma$, as figure:

Simple SEIR Simulation ( with births and deaths)

\begin{eqnarray} \frac{dS}{dt} &= \mu-(\beta I+\mu) S \\ \frac{dE}{dt} &= \beta S I-(\mu+k) E \\ \frac{dI}{dt} &= k E-(\mu+\gamma) I \\ \frac{dR}{dt} &= \gamma I-\mu E \end{eqnarray}



In [16]:

    
mu=1/(70*365.0)
beta=520/365.0
k=1/14.0
gamma=1/7.0
ND=10*365.0
TS=1.0
S0=0.1
E0=1e-4
I0=1e-4
INPUT = (S0, E0, I0)

def diff_eqs(INP,t):  
	'''The main set of equations'''
	Y=np.zeros((3))
	V = INP    
	Y[0] = mu - beta * V[0] * V[2] - mu * V[0]
	Y[1] = beta * V[0] * V[2] - k * V[1] - mu * V[1]
	Y[2] = k * V[1] - gamma * V[2] - mu * V[2]
	return Y   # For odeint



t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

Rec=1. - (RES[:,0]+RES[:,1]+RES[:,2])
#print RES

#Ploting
plt.subplot(311)
plt.plot(RES[:,0], '-g', label='Susceptibles')
plt.title('SEIR with Births and Deaths')
plt.xlabel('Time')
plt.ylabel('Susceptibles')
plt.subplot(312)
plt.plot(RES[:,1], '-m', label='Exposed')
plt.plot(RES[:,2], '-r', label='Infectious')
plt.legend(loc=0)
plt.xlabel('Time')
plt.ylabel('Infected')
plt.subplot(313)
plt.plot(Rec, '-k', label='Recovereds')
plt.xlabel('Time')
plt.ylabel('Recovereds')









    Out[16]:





<matplotlib.text.Text at 0x7fe883341b50>



In [15]:

    
beta=520/365.
gamma=1/7.0
#mu=1/(70*365.0)
mudays=70
def vis3(mudays,beta,gamma):
    #mu=1/(70*365.0)
    #beta=520/365.0
    k=1/14.0
    #gamma=1/7.0
    mu=1/float(mudays*365)
    ND=10*365.0
    TS=1.0
    S0=0.1
    E0=1e-4
    I0=1e-4
    INPUT = (S0, E0, I0)

    def diff_eqs(INP,t):  
        '''The main set of equations'''
        Y=np.zeros((3))
        V = INP    
        Y[0] = mu - beta * V[0] * V[2] - mu * V[0]
        Y[1] = beta * V[0] * V[2] - k * V[1] - mu * V[1]
        Y[2] = k * V[1] - gamma * V[2] - mu * V[2]
        return Y   # For odeint
    
    t_start = 0.0; t_end = ND; t_inc = TS
    t_range = np.arange(t_start, t_end+t_inc, t_inc)
    RES = spi.odeint(diff_eqs,INPUT,t_range)
    Rec=1. - (RES[:,0]+RES[:,1]+RES[:,2])
    
    fig, ax = plt.subplots()
    plt.subplot(311)
    plt.plot(RES[:,0], '-g', label='Susceptibles')
    plt.title('SEIR with Births and Deaths,mu=1/mudys/365')
    plt.xlabel('Time')
    plt.ylabel('Susceptibles')
    plt.subplot(312)
    plt.plot(RES[:,1], '-m', label='Exposed')
    plt.plot(RES[:,2], '-r', label='Infectious')
    plt.legend(loc=0)
    plt.xlabel('Time')
    plt.ylabel('Infected')
    plt.subplot(313)
    plt.plot(Rec, '-k', label='Recovereds')
    plt.xlabel('Time')
    plt.ylabel('Recovereds')
    
    return ax



In [17]:

    
i = interact(vis3,
         mudays  =widgets.FloatSliderWidget(min=1, max=365, step=1, value=70, description="mudays"),    
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=520/365., description="beta"),
         gamma=widgets.FloatSliderWidget(min=0.0, max=5, step=0.1, value=1/7., description="gamma"),
         )

Governed Differential Equations

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

where $C(t)$ is not an epidemiological state but serves to keep track of the cumulative number of Ebola cases from the time of onset of symptoms.

$1/k$ Mean incubation period (days)
$1/\gamma$ Mean infectious period (days)
$\beta=\frac{\beta_0+\beta_1}{2}$

Case for Congo 1995

N=564500
$\beta_0=0.33,\beta_1=0.09$
$1/k=5.3,1/\gamma=5.61$



In [ ]:



In [2]:

    
N=64500.
S0=0.9*N;E0=0.0;I0=N-S0-E0;R0=0.;C0=0.;
beta0=0.33;beta1=0.09;
beta=beta0;gamma=1/5.61

k=1/5.3
n=10000
T=100
dt=T/float(n)
t=np.linspace(0,T,n+1)



In [3]:

    
S,E,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;E[0]=E0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]/N
    E[i+1] = E[i]+dt*(beta*S[i]*I[i]/N-k*E[i])
    I[i+1] = I[i]+dt*(k*E[i]-gamma*I[i]) 
    R[i+1] = N-S[i+1]-I[i+1]-E[i+1]
plt.plot(t,S/N,t,I/N,'r--',t,E/N,'y',t,R/N,'k-') 
#plt.ylim(0,1)
plt.legend(['Suspective','Infected','Exposed','Recoverd'])









    Out[3]:





<matplotlib.legend.Legend at 0x106470190>



In [26]:

    
beta/gamma









    Out[26]:





1.8513000000000002



In [27]:

    
J=np.matrix([[-k, beta],[k, -gamma]]);J









    Out[27]:





matrix([[-0.18867925,  0.33      ],
        [ 0.18867925, -0.17825312]])



In [9]:

    
import scipy.linalg as la



In [29]:

    
r=np.max(np.real(la.eigvals(J)));r









    Out[29]:





0.06611612264241562



In [32]:

    
#r=0.07
1+(r*r+(k+gamma)*r)/k/gamma









    Out[32]:





1.9093917

Intervention Case

After diease breakout, there ahould be some intervention strategies made to control the spread of Ebola.

The transmission rate, $\beta$, could be modeled reasonably as follows:

$$ \beta ( t ) = \left\{\begin{array}{l} \beta_{0},\qquad\qquad\qquad\qquad t \leqslant t_{I}\\ \beta_{1} + ( \beta_{0} - \beta_{1} ) e^{-q ( t-t_{I} )}, \quad t \geqslant t_{I} \end{array}\right. $$

i.e. decreases fom $\beta_0$ to $\beta_1$ at $t_I$, the time at which interventions start, and continus to be decreaing exponentially.

Control of Transition rate, $q$

Another interpretation of the parameter q can be given in terms $$t_h=t_I+ln(2)/q$$ where $t_h$ is the time at which $\beta(t)=\frac{\beta_1+\beta_0}{2}$



In [4]:

    
q=0.5
tI=5
th=tI+np.log(2)/q
beta1+(beta0-beta1)*np.exp(-q*(th-tI))









    Out[4]:





0.20999999999999999



In [5]:

    
def bt(b0,b1,tI,q,T,n):
    beta0=b0
    beta1=b1
    dt=T/float(n)
    t=np.linspace(0,T,n+1)
    bI=np.zeros(n+1)
    nI=np.ceil(n*tI/T)
    bI[:nI]=beta0*np.ones(nI)
    bI[nI:]=beta1+(beta0-beta1)*np.exp(-q*(t[nI:]-tI))
    return bI



In [9]:

    
bI=bt(beta0,beta1,5,0.29,T,n)



In [10]:

    
plt.xlim(0,100)
plt.plot(t,bI,t,(beta0+beta1)/2.*np.ones(np.size(t)))









    Out[10]:





[<matplotlib.lines.Line2D at 0x106862650>,
 <matplotlib.lines.Line2D at 0x106ac2b10>]



In [11]:

    
def BRN(b):
    s=np.size(b)
    R=np.zeros(s)
    i=0
    for i in np.arange(s):
        J=np.matrix([[-k, b[i]],[k, -gamma]]);
        r=np.max(np.real(la.eigvals(J)));
        R[i]=1+(r*r+(k+gamma)*r)/k/gamma
    return R



In [12]:

    
R00=BRN(bI)



In [13]:

    
plt.plot(t,R00)









    Out[13]:





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



In [26]:

    
S,E,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;E[0]=E0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]/N
    E[i+1] = E[i]+dt*(beta*S[i]*I[i]/N-k*E[i])
    I[i+1] = I[i]+dt*(k*E[i]-gamma*I[i]) 
    R[i+1] = N-S[i+1]-I[i+1]-E[i+1]
plt.plot(t,S/N,t,I/N,'r--',t,E/N,'y',t,R/N,'k-',t,R00) 
#plt.ylim(0,1)
plt.legend(['Suspective','Infected','Exposed','Recoverd','$R_0$'])









    Out[26]:





<matplotlib.legend.Legend at 0x106a5da90>

Symboilc Solution dor ODE's

Up to now, we only solve ODE's by numeric schemes (approximately); can we solve them theoretically? Just as well-know fact:

 Not all integrals can be integrated!

Yes, "Definite No"!

But we still have chance to solve ODE's by "sympy"!

Example (Exponential Growth)

The solution of $y'(t) = y(t), y(0) = y0$ is $y(t)=y0 \exp(t)$



In [35]:

    
from sympy import Symbol , dsolve , Function , Derivative , Eq, lambdify, symbols



In [20]:

    
x=Symbol('x')
y0=Symbol('y0')
y=Function('y')



In [21]:

    
f_=Derivative(y(x),x)-y(x)



In [22]:

    
sol=dsolve(Derivative(y(x),x)-y(x),y(x),ics={y(0):y0})



In [23]:

    
print(sol)









    



y(x) == C1*exp(x)



In [24]:

    
C1=Symbol('C1')



In [25]:

    
ysol=sol.subs({C1:y0})



In [26]:

    
print(ysol)









    



y(x) == y0*exp(x)



In [32]:

    
z= lambdify(x,ysol.rhs)



In [49]:

    
z= lambdify(x,y0*exp(x))



In [33]:

    
z(1)









    Out[33]:





2.71828182845905*y0



In [36]:

    
A,N=symbols("A N")
g_=Derivative(y(x),x)-A*y(x)*(N-y(x))



In [37]:

    
sol2=dsolve(g_,y(x))



In [38]:

    
sol2









    Out[38]:





y(x) == N*exp(N*(A*x + C2))/(exp(N*(A*x + C1)) - 1)



In [ ]:



In [ ]: