# Program To Create a Mathematical Model in Python Language Assignment Solution.

## Instructions

Objective
Write a program to create a mathematical model in python language.

## 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)