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