Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah
import numpy as np
from numpy import *
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
from scipy.integrate import odeint
The SIR model for epidemics consists of three ODEs \begin{equation} \frac{\text{d}S}{\text{d}t} = -\beta I S \end{equation} \begin{equation} \frac{\text{d}I}{\text{d}t} = \beta I S - \gamma I \end{equation} \begin{equation} \frac{\text{d}R}{\text{d}t} = \gamma I \end{equation}
def rhs_sir(PHI, t, N, beta, gamma):
S = PHI[0]
I = PHI[1]
R = PHI[2]
dSdt = -beta/N * S * I
dIdt = beta/N * S * I - gamma * I
dRdt = gamma * I
return [dSdt, dIdt, dRdt]
N = 2e3 # total population
I0 = 1 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
β = 0.6
γ = 1.0/14.0
tend = 100 # days
t = np.linspace(0,tend,int(tend))
# Initial conditions vector
y0 = [S0, I0, R0]
sol = odeint(rhs_sir, y0, t, args=(N, β, γ))
S = sol[:,0]
I = sol[:,1]
R = sol[:,2]
plt.plot(t, S, 'b', alpha=0.65, lw=2, label='Susceptible')
plt.plot(t, I, 'r', alpha=0.65, lw=2, label='Infected')
plt.plot(t, R, 'g', alpha=0.65, lw=2, label='Recovered with immunity')
plt.xlabel('# days')
plt.ylabel('Fraction of Population')
plt.xlim(0,tend)
plt.ylim(0,N)
plt.legend()
plt.title('SIR Model for Pandemics')
plt.grid()
# plt.savefig('SIR model.pdf')
# plt.show()
def plot_sir(β,days,tend,semilogy=False):
N = 1000 # total population
I0 = 1 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
tend = tend*30 # tend is given in months
t = np.linspace(0,tend,int(tend))
γ = 1.0/days
# Initial conditions vector
y0 = [S0, I0, R0]
# Integrate the SIR equations over the time grid, t.
sol = odeint(rhs_sir, y0, t, args=(N, β, γ))
[S, I, R] = sol.T
f = plt.plot
if (semilogy):
f =plt.semilogy
f(t/30, S, 'b', alpha=0.65, lw=2, label='Susceptible')
f(t/30, I, 'r', alpha=0.65, lw=2, label='Infected')
f(t/30, R, 'g', alpha=0.65, lw=2, label='Recovered')
plt.xlabel('# months since outbreak')
plt.ylabel('# Population (1000s)')
# plt.xlim(0,tend/30)
plt.ylim(bottom=1e-4)
plt.legend()
plt.title('SIR Model for Epidemics')
plt.grid()
# plt.savefig('SIR model.pdf')
# plt.show()
import ipywidgets as widgets
from ipywidgets import interact, interact_manual
style = {'description_width': 'initial'}
β = widgets.FloatSlider(value=0.143,min=0.001,max=1,step=0.001,description='$\\beta$:',continuous_update=False)
days = widgets.FloatSlider(value=14,min=1,max=24,step=0.1, description='Recovery (days):' ,style=style,continuous_update=False)
tend = widgets.FloatSlider(value=12,min=1,max=32,step=0.5, description='Simulation Time (months):',style=style,continuous_update=False)
semilogy = widgets.Checkbox(value=False,description='semilogy plot', style=style)
ui1 = widgets.HBox([β,days])
ui2 = widgets.HBox([tend,semilogy])
out = widgets.interactive_output(plot_sir, {'β': β, 'days': days, 'tend': tend, 'semilogy':semilogy})
display(ui1,ui2, out)
HBox(children=(FloatSlider(value=0.143, continuous_update=False, description='$\\beta$:', max=1.0, min=0.001, …
HBox(children=(FloatSlider(value=12.0, continuous_update=False, description='Simulation Time (months):', max=3…
Output()
The SIR model for epidemics consists of three ODEs \begin{equation} \frac{\text{d}S}{\text{d}t} = -\beta I S + \mu R \end{equation} \begin{equation} \frac{\text{d}I}{\text{d}t} = \beta I S - \gamma I \end{equation} \begin{equation} \frac{\text{d}R}{\text{d}t} = \gamma I - \mu R \end{equation}
def rhs_sir_immunity(SIR, t, N, β, γ,μ):
[S, I, R] = SIR #SIR is a vector of solutions
dSdt = -β/N * S * I + μ*R
dIdt = β/N * S * I - γ * I
dRdt = γ * I - μ*R
return [dSdt, dIdt, dRdt]
N = 3e3 # total population
I0 = 1 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
β = 0.6
γ = 1/20.0
μ = 0.001
tend = 700 # days
t = np.linspace(0,tend,tend)
# Initial conditions vector
y0 = [S0, I0, R0]
# Integrate the SIR equations over the time grid, t.
sol = odeint(rhs_sir_immunity, y0, t, args=(N, β, γ,μ))
[S, I, R] = sol.T
plt.plot(t, S/N, 'b', alpha=0.65, lw=2, label='Susceptible')
plt.plot(t, I/N, 'r', alpha=0.65, lw=2, label='Infected')
plt.plot(t, R/N, 'g', alpha=0.65, lw=2, label='Recovered with immunity')
plt.xlabel('# days')
plt.ylabel('Fraction of Population')
plt.xlim(0,tend)
plt.ylim(0,1.1)
plt.legend()
plt.title('SIR Model for Pandemics')
plt.grid()
# plt.savefig('SIR model.pdf')
# plt.show()
plt.plot(S,I)
[<matplotlib.lines.Line2D at 0x7f8a58dc2400>]
def rhs(SIRD, t, N, beta, gamma, δ1, δ2, vaccinateAfter, minI,maxI):
S, I, R, D = SIRD
if (distanceModel == 'Reactive'):
δf2 = 0. #δ2
δf1 = 0. #δ1
if (I>=maxI):
δf1 = 0.0
δf2 = δ2 #* np.sin(2.0*np.pi*t/200)**2
if (I <= minI):
δf1 = δ1 #- beta/N * I #* np.sin(2.0*np.pi*t/200)**2
δf2 = 0.0
else:
δf2 = δ2
δf1 = δ1
vacc = 0.0
if (t>= vaccinateAfter):
vacc = 1.0
δf2 = 0.0 #* np.sin(2.0*np.pi*t/20)
δf1 = 1.0 #* np.sin(2.0*np.pi*t/200)**2
beta = 0.0
dSdt = - beta/N * S * I + δf1*D - δf2*S - vacc*S
dIdt = beta/N * S * I - gamma * I
dRdt = gamma * I + vacc*S
dDdt = - δf1*D + δf2*S
return dSdt, dIdt, dRdt, dDdt
def rhs_sird(PHI, t, N, beta, gamma, d1):
S, I, R, D = PHI
d2 = 1.0 - d1
maxI = 310
minI = 300
dIdt = beta/N * S * I - gamma * I
if (dIdt>=0.8):
d1 = 0.75
elif (dIdt < -0.1):
d1 = 0.2
else:
d1 = 0.5
d2 = 1.0 - d1
dSdt = - beta/N * S * I + d2*D - d1*S
dRdt = gamma * I
dDdt = - d2*D + d1*S
return [dSdt, dIdt, dRdt, dDdt]
N = 1000 # total population
I0 = 1 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
D0 = N - S0 - I0 - R0
β = 0.5
γ = 1/16.0
d1 = 0.5
tend = 60 # days
t = np.linspace(0,tend,1000)
# Initial conditions vector
y0 = [S0, I0, R0, D0]
# Integrate the SIR equations over the time grid, t.
sol = odeint(rhs_sird, y0, t, args=(N, β, γ, d1))
[S, I, R, D] = sol.T
/opt/anaconda3/lib/python3.9/site-packages/scipy/integrate/odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)
plt.plot(t, S, 'b', alpha=0.65, lw=2, label='Susceptible')
plt.plot(t, I, 'r', alpha=0.65, lw=2, label='Infected')
plt.plot(t, R, 'g', alpha=0.65, lw=2, label='Recovered with immunity')
plt.plot(t, D, 'k', alpha=0.65, lw=2, label='Distanced')
plt.xlabel('# days')
plt.ylabel('Fraction of Population')
plt.xlim(0,tend)
plt.ylim(0,N)
plt.legend()
plt.title('SIR Model for Pandemics')
plt.grid()
# plt.savefig('SIR model.pdf')
# plt.show()