Instructions
Requirements and Specifications
Source Code
## Required Imports
import numpy as np
import matplotlib.pyplot as plt
## Brief Explanation of the method applied in this assignment
In this assignment, we will simulate the Spread of an Infectious Disease, using a model knwon as Susceptible, Infected and Recovered (SIR) model. The model assumes that a population P contains:
* s: Number of people who are susceptible to infection
* i: The number of people who are already infected
* r: The number of people who have recovered from the infection
So, these three parameters conform a population P, which means that the following relation holds:
$P = s + i + r$
The model is described by the following differential equations:
$\frac{ds}{dt}=-asi$
$\frac{di}{dt}=asi-bi$
$\frac{dr}{dr}=bi$
Where **a** is the rate of infection, and **b** the rate of recovery
The model will be solved using Euler's Method for Integration. Given an initial state $s_{0}, i_{0}, r_{0}$ which defines the initial values of the parameters, and given the following functions:
$x(s_{n}) = -asi$
$y(i_{n}) = asi-bi$
$z(r_{n}) = bi$
And given a time-step $h$, the new state $n+1$ can be calculated as:
$s_{n+1} = s_{n} + hx(s_{n})$
$i_{n+1} = i_{n} + hy(s_{n})$
$r_{n+1} = r_{n} + hz(s_{n})$
## Implementation
### Function that stores all three differential equations
def f(s, i, r, a, b):
df = np.zeros((3,1))
df[0] = -a*s*i
df[1] = a*s*i - b*i
df[2] = b*i
return df
### Case 1: 0% chance of recovery. All population must end infected
a = 0.00001
b = 0.0
s0 = 44358
i0 = 3
r0 = 0
n_days = 100
h = 0.01 # step
N = int(n_days/h) # number of steps
t = np.linspace(0, n_days, N)
# Now, simulate
s = np.zeros(N)
i = np.zeros(N)
r = np.zeros(N)
s[0] = s0
i[0] = i0
r[0] = r0
for n in range(N-1):
x_old = np.array([s[n], i[n], r[n]]).reshape((3,1))
x_new = x_old + h*f(s[n], i[n], r[n], a, b)
s[n+1] = x_new[0]
i[n+1] = x_new[1]
r[n+1] = x_new[2]
# Now, plot
plt.figure()
plt.plot(t, s, label = 'S')
plt.plot(t, i, label = 'I')
plt.plot(t, r, label = 'R')
plt.legend()
plt.grid(True)
plt.show()
We see that, because there is a 0% chance of recovery, all the population ends infected
## Case 2: 0% chance of infection, but assume that there are 20000 persons infected in the beginning
a = 0.0000
b = 0.025
s0 = 24358
i0 = 20000
r0 = 0
n_days = 100
h = 0.01 # step
N = int(n_days/h) # number of steps
t = np.linspace(0, n_days, N)
# Now, simulate
s = np.zeros(N)
i = np.zeros(N)
r = np.zeros(N)
s[0] = s0
i[0] = i0
r[0] = r0
for n in range(N-1):
x_old = np.array([s[n], i[n], r[n]]).reshape((3,1))
x_new = x_old + h*f(s[n], i[n], r[n], a, b)
s[n+1] = x_new[0]
i[n+1] = x_new[1]
r[n+1] = x_new[2]
# Now, plot
plt.figure()
plt.plot(t, s, label = 'S')
plt.plot(t, i, label = 'I')
plt.plot(t, r, label = 'R')
plt.legend()
plt.grid(True)
plt.show()
We see that, because there is a 0% chance of infection, all the susceptibles ones remains in that state, while the infected gets recovered
## Effect of the infection rate a, on the sensivity of the model, performing 100 experiments
# Define b
b = 0.025
# Define control value
a_control = 0.00001
a_std = 0.2*a_control # 20%
# Initial values
s0 = 24358
i0 = 20000
r0 = 0
n_days = 100
h = 0.01 # step
N = int(n_days/h) # number of steps
t = np.linspace(0, n_days, N)
# Generate the values for a
a_values = np.random.normal(a_control, a_std, 100)
# List to store the peaks of infection for each simulation
i_peaks = list()
#a_values = np.abs(a_values)
# Simulate for each a
plt.figure()
for i in range(len(a_values)):
a = a_values[i]
# Now, simulate
s = np.zeros(N)
i = np.zeros(N)
r = np.zeros(N)
s[0] = s0
i[0] = i0
r[0] = r0
for n in range(N-1):
x_old = np.array([s[n], i[n], r[n]]).reshape((3,1))
x_new = x_old + h*f(s[n], i[n], r[n], a, b)
s[n+1] = x_new[0]
i[n+1] = x_new[1]
r[n+1] = x_new[2]
# Now, plot
plt.plot(t, i, color = 'grey', linewidth=0.5)
i_peaks.append(max(i))
# Now simulate for the mean
a = np.mean(a_values)
s = np.zeros(N)
i = np.zeros(N)
r = np.zeros(N)
s[0] = s0
i[0] = i0
r[0] = r0
for n in range(N-1):
x_old = np.array([s[n], i[n], r[n]]).reshape((3,1))
x_new = x_old + h*f(s[n], i[n], r[n], a, b)
s[n+1] = x_new[0]
i[n+1] = x_new[1]
r[n+1] = x_new[2]
# Now, plot
plt.plot(t, i, linewidth = 2, color = 'red', label = 'Mean')
plt.legend()
plt.grid(True)
plt.show()
## Histogram showing Frequency of Peaks of Infection
plt.figure()
plt.hist(i_peaks)
plt.xlabel('Peak Infections')
plt.ylabel('Frequency')
plt.grid(True)