Select Git revision
Forked from
orestis.malaspin / math_tech_info
Source project has a limited visibility.
covid.py 1.88 KiB
import numpy as np
import matplotlib.pyplot as plt
swiss = np.array([2.100350058, 3.150525088, 4.900816803, 6.534422404, 10.501750292, 13.302217036, 24.970828471, 31.271878646, 39.323220537, 43.640606768, 57.292882147, 76.079346558, 100.116686114, 131.271878646, 158.576429405, 256.709451575, 256.709451575, 268.378063011, 309.218203034, 353.325554259, 453.675612602])
swiss = np.array([18, 27, 42, 56, 90, 114, 214, 268, 337, 374, 491, 652, 858, 1125, 1359, 2200, 2200, 2300, 2650, 3028, 3888]) / 0.02
days = np.array(range(1,len(swiss)+1))
def s(dt, beta, S0, I0, N):
return S0 - dt * (beta * S0 * I0 / N)
def ii(dt, beta, lamb, S0, I0, N):
return I0 + dt * (beta * S0 * I0 / N - lamb * I0)
def r(dt, lamb, R0, I0, N):
return R0 + dt * lamb * I0
def compute_n(S, R, I):
return S + R + I
def timestep(S0, I0, R0, dt, beta, lamb, N):
S1 = s(dt, beta, S0, I0, N)
I1 = ii(dt, beta, lamb, S0, I0, N)
R1 = r(dt, lamb, R0, I0, N)
return S1, I1, R1
S0 = 8000000
I0 = swiss[0]
R0 = 0
max_t = 3*days[len(swiss)-1]
n_steps = 100000
dt = max_t / n_steps
N = compute_n(S0, I0, R0)
lamb = 1.0 / 14.0
beta_1 = 0.35
beta_2 = beta_1 / 10
beta = beta_1
R0 = beta / lamb
print(R0)
s_list = [S0]
r_list = [R0]
i_list = [I0]
t_list = [1]
for i in range(0, n_steps):
S1, I1, R1 = timestep(s_list[i], i_list[i], r_list[i], dt, beta, lamb, N)
s_list.append(S1)
i_list.append(I1)
r_list.append(R1)
t_list.append(t_list[i]+dt)
# if (R1 >= 0.9 * N):
# print(t_list[i]+dt)
# break
# if (t_list[i+1] < 27) :
# beta = beta_1
# elif (I1 > 0.01*N) :
# beta = beta_2
# else:
# beta = beta_1
s = np.array(s_list)
r = np.array(r_list)
ii = np.array(i_list)
t = np.array(t_list)
s
plt.semilogy(t, s, 'b')
plt.semilogy(t, r, 'r')
plt.semilogy(t, ii, 'k')
plt.semilogy(days, swiss, 'k*')
plt.legend(['S', 'I', 'R', 'swiss'])
plt.show()