# System of Equations Example - SIR Epidemic Model
## CH EN 2450 - Numerical Methods
**Prof. Tony Saad (<a>www.tsaad.net</a>) <br/>Department of Chemical Engineering <br/>University of Utah**
<hr/>


```python
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

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}


```python
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]
```


```python
N  = 3e3 # total population
I0 = 10 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
β  = 0.6
γ  = 1/20.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]
```


```python
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()
```


    

    


# Interactive Simulation


```python
def plot_sir(β,days,tend,semilogy=False):
    N = 10 # 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 - with loss of immunity

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}


```python
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]
```


```python
N = 3e3 # total population
I0 = 10 # 10 sick individuals
R0 = 0 # 0 people recovered
S0 = N - I0 # we need to conserve the total population
β = 0.6
γ = 1/20.0
μ = 1/320.0

tend = 600 # 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
```


```python
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()
```


    

    



```python
plt.plot(S,I)
```




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




    

    



```python

```

0887260d3211a4b82ff353e4e92b1975462d56ec
ipynb
Jupyter Notebook
topics/initial-value-problems/System of Equations Example - SIR Epidemic Model.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-03-27T05:22:34.000Z
2021-01-27T10:49:13.000Z
topics/initial-value-problems/System of Equations Example - SIR Epidemic Model.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
topics/initial-value-problems/System of Equations Example - SIR Epidemic Model.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-12-29T23:31:56.000Z
2021-12-28T19:04:10.000Z
Qwen/Qwen-72B
1. YES
2. YES