import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
from scipy.optimize import root
from scipy.linalg import eig
#standard values of parameters
N = 1.0
b = 0.018
alpha = 0.33
gamma = 0.2
e = 0.2
lamb = 0.1
beta0 = 1.68
mu = 0.1
kappa = 1117
#equations
def beta(P):
return beta0*(1-mu)*np.exp(-kappa*P/N)
def Sdot(S, I, P):
return b*(N-S) - (beta(P)*S*I)/N
def Edot(S, E, I, P):
return beta(P)*S*I/N - (alpha + b)*E
def Idot(E, I):
return alpha*E - (gamma + b)*I
def Rdot(I, R):
return gamma*I - b*R
def Pdot(I, P):
return e*gamma*I - lamb*P
# Runge Kutta 4 algorithm
def RK4(S0, E0, I0, R0, P0, n, h=0.01):
#n = amount of integration steps
#h = size of time step
assert abs(S0+E0+I0+R0 - N) < 1e-8, 'S0+E0+I0+R0 should be equal to N' #check if the starting conditions are valid
#prepare arrays that will be returned and set first values equal to starting conditions
S = np.zeros(n)
S[0] = S0
E = np.zeros(n)
E[0] = E0
I = np.zeros(n)
I[0] = I0
R = np.zeros(n)
R[0] = R0
P = np.zeros(n)
P[0] = P0
#prepare arrays for derivatives in the calculated points and set first values equal to derivatives in the starting conditions
dS = np.zeros(n)
dS[0] = Sdot(S0, I0, P0)
dE = np.zeros(n)
dE[0] = Edot(S0, E0, I0, P0)
dI = np.zeros(n)
dI[0] = Idot(E0, I0)
dR = np.zeros(n)
dR[0] = Rdot(I0, R0)
dP = np.zeros(n)
dP[0] = Pdot(I0, P0)
#prepare time array
t = np.zeros(n)
#runge kutta algorithm
for i in range(1, n):
k1 = h*Sdot(S[i-1], I[i-1], P[i-1])
l1 = h*Edot(S[i-1], E[i-1], I[i-1], P[i-1])
m1 = h*Idot(E[i-1], I[i-1])
n1 = h*Rdot(I[i-1], R[i-1])
o1 = h*Pdot(I[i-1], P[i-1])
k2 = h*Sdot(S[i-1]+k1/2, I[i-1]+m1/2, P[i-1]+o1/2)
l2 = h*Edot(S[i-1]+k1/2, E[i-1]+l1/2, I[i-1]+m1/2, P[i-1]+o1/2)
m2 = h*Idot(E[i-1]+l1/2, I[i-1]+m1/2)
n2 = h*Rdot(I[i-1]+m1/2, R[i-1]+n1/2)
o2 = h*Pdot(I[i-1]+m1/2, P[i-1]+o1/2)
k3 = h*Sdot(S[i-1]+k2/2, I[i-1]+m2/2, P[i-1]+o2/2)
l3 = h*Edot(S[i-1]+k2/2, E[i-1]+l2/2, I[i-1]+m2/2, P[i-1]+o2/2)
m3 = h*Idot(E[i-1]+l2/2, I[i-1]+m2/2)
n3 = h*Rdot(I[i-1]+m2/2, R[i-1]+n2/2)
o3 = h*Pdot(I[i-1]+m2/2, P[i-1]+o2/2)
k4 = h*Sdot(S[i-1]+k3, I[i-1]+m3, P[i-1]+o3)
l4 = h*Edot(S[i-1]+k3, E[i-1]+l3, I[i-1]+m3, P[i-1]+o3)
m4 = h*Idot(E[i-1]+l3, I[i-1]+m3)
n4 = h*Rdot(I[i-1]+m3, R[i-1]+n3)
o4 = h*Pdot(I[i-1]+m3, P[i-1]+o3)
S[i] = S[i - 1] + (k1+(k2*2)+(k3*2)+k4)/6
E[i] = E[i - 1] + (l1+(l2*2)+(l3*2)+l4)/6
I[i] = I[i - 1] + (m1+(m2*2)+(m3*2)+m4)/6
R[i] = R[i - 1] + (n1+(n2*2)+(n3*2)+n4)/6
P[i] = P[i - 1] + (o1+(o2*2)+(o3*2)+o4)/6
t[i] = t[i - 1] + h
dS[i] = Sdot(S[i], I[i], P[i])
dE[i] = Edot(S[i], E[i], I[i], P[i])
dI[i] = Idot(E[i], I[i])
dR[i] = Rdot(I[i], R[i])
dP[i] = Pdot(I[i], P[i])
return S, E, I, R, P, t, dS, dE, dI, dR, dP
#using the RK4 algorithm for certain starting conditions
#evolution towards a stable point
S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)
fig = plt.figure(figsize=(13,10))
ax1 = fig.add_subplot(321)
ax1.plot(t, S, color='r')
ax1.set_xlabel('t')
ax1.set_ylabel('S')
ax2 = fig.add_subplot(322)
ax2.plot(t, E, color='g')
ax2.set_xlabel('t')
ax2.set_ylabel('E')
ax3 = fig.add_subplot(323)
ax3.plot(t, I, color='k')
ax3.set_xlabel('t')
ax3.set_ylabel('I')
ax4 = fig.add_subplot(324)
ax4.plot(t, R, color='y')
ax4.set_xlabel('t')
ax4.set_ylabel('R')
ax5 = fig.add_subplot(325)
ax5.plot(t, P, color='c')
ax5.set_xlabel('t')
ax5.set_ylabel('P')
ax6 = fig.add_subplot(326)
ax6.plot(t, S+E+I+R-N, color='m')
ax6.set_xlabel('t')
ax6.set_ylabel('S+E+I+R-N')
plt.show()
#evolution towards a limit cycle
alpha = 0.33
beta0 = 100
S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.82025, E0=0.0075, I0=0.01125, R0=0.1610, P0=0.0057, n=75000, h=0.01)
fig = plt.figure(figsize=(13,10))
ax1 = fig.add_subplot(321)
ax1.plot(t, S, color='r')
ax1.set_xlabel('t')
ax1.set_ylabel('S')
ax2 = fig.add_subplot(322)
ax2.plot(t, E, color='g')
ax2.set_xlabel('t')
ax2.set_ylabel('E')
ax3 = fig.add_subplot(323)
ax3.plot(t, I, color='k')
ax3.set_xlabel('t')
ax3.set_ylabel('I')
ax4 = fig.add_subplot(324)
ax4.plot(t, R, color='y')
ax4.set_xlabel('t')
ax4.set_ylabel('R')
ax5 = fig.add_subplot(325)
ax5.plot(t, P, color='c')
ax5.set_xlabel('t')
ax5.set_ylabel('P')
plt.show()
plt.plot(t, S+E+I+R-N)
plt.title(r'Error on the condition $N=S+E+I+R$')
plt.xlabel('t')
plt.ylabel(r'S+E+I+R-N')
plt.show()
alpha = 0.33
beta0 = 1.68
mu=0
S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)
mu=0.5
S1, E1, I1, R1, P1, t1, dS1, dE1, dI1, dR1, dP1 = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)
mu=1.0
S2, E2, I2, R2, P2, t2, dS2, dE2, dI2, dR2, dP2 = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)
plt.plot(t, I, label=r'$\mu=0$', color='r')
plt.plot(t1, I1, label=r'$\mu=0.5$', color='g')
plt.plot(t2, I2, label=r'$\mu=1$', color='b')
plt.title('Influence of governmental measures')
plt.xlabel('t')
plt.ylabel('I')
plt.legend()
plt.show()
mu = 0.1 #reset mu
alpha = 0.33
beta0 = 70
S0, E0, I0, R0, P0, t0, dS0, dE0, dI0, dR0, dP0 = RK4(S0=0.8483, E0=0.0031, I0=0.0058, R0=0.1428, P0=0.0047, n=25000, h=0.01)
beta0 = 55
S1, E1, I1, R1, P1, t1, dS1, dE1, dI1, dR1, dP1 = RK4(S0=0.8561, E0=0.0050, I0=0.0075, R0=0.1314, P0=0.0044, n=25000, h=0.01)
beta0 = 60
S2, E2, I2, R2, P2, t2, dS2, dE2, dI2, dR2, dP2 = RK4(S0=0.8479, E0=0.0032, I0=0.0089, R0=0.1400, P0=0.0054, n=25000, h=0.01)
beta0 = 50
S3, E3, I3, R3, P3, t3, dS3, dE3, dI3, dR3, dP3 = RK4(S0=0.8567, E0=0.0055, I0=0.0093, R0=0.1285, P0=0.0046, n=25000, h=0.01)
beta0 = 75
S4, E4, I4, R4, P4, t4, dS4, dE4, dI4, dR4, dP4 = RK4(S0=0.8258, E0=0.0048, I0=0.0164, R0=0.1530, P0=0.0067, n=25000, h=0.01)
beta0 = 80
S5, E5, I5, R5, P5, t5, dS5, dE5, dI5, dR5, dP5 = RK4(S0=0.8123, E0=0.0228, I0=0.0221, R0=0.1428, P0=0.0047, n=25000, h=0.01)
beta0 = 85
S6, E6, I6, R6, P6, t6, dS6, dE6, dI6, dR6, dP6 = RK4(S0=0.8064, E0=0.0236, I0=0.0240, R0=0.1460, P0=0.0049, n=25000, h=0.01)
beta0 = 90
S7, E7, I7, R7, P7, t7, dS7, dE7, dI7, dR7, dP7 = RK4(S0=0.8314, E0=0.0016, I0=0.0074, R0=0.1596, P0=0.0061, n=25000, h=0.01)
beta0 = 95
S8, E8, I8, R8, P8, t8, dS8, dE8, dI8, dR8, dP8 = RK4(S0=0.8110, E0=0.0048, I0=0.0188, R0=0.1654, P0=0.0075, n=25000, h=0.01)
beta0 = 100
S9, E9, I9, R9, P9, t9, dS9, dE9, dI9, dR9, dP9 = RK4(S0=0.7932, E0=0.0227, I0=0.0288, R0=0.1553, P0=0.0057, n=25000, h=0.01)
#print final positions and use these as starting conditions until we stay on a limit cycle
'''
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S0[-1], E0[-1], I0[-1], R0[-1], P0[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S1[-1], E1[-1], I1[-1], R1[-1], P1[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S2[-1], E2[-1], I2[-1], R2[-1], P2[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S3[-1], E3[-1], I3[-1], R3[-1], P3[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S4[-1], E4[-1], I4[-1], R4[-1], P4[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S5[-1], E5[-1], I5[-1], R5[-1], P5[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S6[-1], E6[-1], I6[-1], R6[-1], P6[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S7[-1], E7[-1], I7[-1], R7[-1], P7[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S8[-1], E8[-1], I8[-1], R8[-1], P8[-1]))
print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S9[-1], E9[-1], I9[-1], R9[-1], P9[-1]))
'''
#determine mins and maxes to use in bifurcation diagram
beta0_values = [70, 55, 62, 50, 75, 80, 85, 90, 95, 100]
maxes_beta = [np.max(I0), np.max(I1), np.max(I2), np.max(I3), np.max(I4), np.max(I5), np.max(I6), np.max(I7), np.max(I8), np.max(I9)]
mins_beta = [np.min(I0), np.min(I1), np.min(I2), np.min(I3), np.min(I4), np.min(I5), np.min(I6), np.min(I7), np.min(I8), np.min(I9)]
#extra integrations to show spiraling towards the limit cycle/stable point
beta0 = 100
S89, E89, I89, R89, P89, t89, dS89, dE89, dI89, dR89, dP89 = RK4(S0=0.82025, E0=0.0075, I0=0.01125, R0=0.1610, P0=0.0057, n=75000, h=0.01)
S91, E91, I91, R91, P91, t91, dS91, dE91, dI91, dR91, dP89 = RK4(S0=0.7475, E0=0.0427, I0=0.0488, R0=0.1610, P0=0.0057, n=75000, h=0.01)
beta0 = 20
Ss, Es, Is, Rs, Ps, ts, dSs, dEs, dIs, dRs, dPs = RK4(S0=0.7932, E0=0.0227, I0=0.0288, R0=0.1553, P0=0.0057, n=75000, h=0.01)
plt.plot(E89, I89, color='b')
plt.plot(E91, I91, color='r')
plt.plot(0.0075, 0.01125, 'ko')
plt.plot(0.0427, 0.0488, 'ko')
plt.plot(E9, I9, color='k')
plt.ylabel('I')
plt.xlabel('E')
plt.title(r'$\beta_0=100$')
plt.show()
plt.plot(Es, Is, color='b')
plt.plot(0.0227, 0.0288, 'ko')
plt.ylabel('I')
plt.xlabel('E')
plt.title(r'$\beta_0=20$')
plt.show()
plt.plot(E3, I3, color='b', label=r'$\beta_0=50$')
plt.plot(E1, I1, color='r', label=r'$\beta_0=55$')
plt.plot(E2, I2, color='k', label=r'$\beta_0=60$')
plt.plot(E0, I0, color='g', label=r'$\beta_0=70$')
plt.plot(E4, I4, color='m', label=r'$\beta_0=75$')
plt.ylabel('I')
plt.xlabel('E')
plt.legend()
plt.show()
plt.plot(E5, I5, color='g', label=r'$\beta_0=80$')
plt.plot(E6, I6, color='r', label=r'$\beta_0=85$')
plt.plot(E7, I7, color='k', label=r'$\beta_0=90$')
plt.plot(E8, I8, color='b', label=r'$\beta_0=95$')
plt.plot(E9, I9, color='m', label=r'$\beta_0=100$')
plt.ylabel('I')
plt.xlabel('E')
plt.legend()
plt.show()
beta0 = 50
alpha = 0.20
S10, E10, I10, R10, P10, t10, dS10, dE10, dI10, dR10, dP10 = RK4(S0=0.853318, E0=0.011345, I0=0.009508, R0=0.125829, P0=0.004264, n=50000, h=0.005)
alpha = 0.25
S11, E11, I11, R11, P11, t11, dS11, dE11, dI11, dR11, dP11 = RK4(S0=0.851354, E0=0.005885, I0=0.010236, R0=0.132525, P0=0.005148, n=50000, h=0.005)
alpha = 0.30
S12, E12, I12, R12, P12, t12, dS12, dE12, dI12, dR12, dP12 = RK4(S0=0.849711, E0=0.006169, I0=0.012282, R0=0.131838, P0=0.005212, n=50000, h=0.005)
alpha = 0.22
S13, E13, I13, R13, P13, t13, dS13, dE13, dI13, dR13, dP13 = RK4(S0=0.855633, E0=0.009476, I0=0.008403, R0=0.126488, P0=0.004150, n=50000, h=0.005)
alpha = 0.27
S14, E14, I14, R14, P14, t14, dS14, dE14, dI14, dR14, dP14 = RK4(S0=0.855109, E0=0.005119, I0=0.008625, R0=0.131147, P0=0.004823, n=50000, h=0.005)
alpha = 0.34
S15, E15, I15, R15, P15, t15, dS15, dE15, dI15, dR15, dP15 = RK4(S0=0.853360, E0=0.007786, I0=0.011758, R0=0.127096, P0=0.004554, n=50000, h=0.005)
alpha = 0.24
S16, E16, I16, R16, P16, t16, dS16, dE16, dI16, dR16, dP16 = RK4(S0=0.856571, E0=0.008556, I0=0.008048, R0=0.126825, P0=0.004107, n=50000, h=0.005)
alpha = 0.195
S17, E17, I17, R17, P17, t17, dS17, dE17, dI17, dR17, dP17 = RK4(S0=0.852537, E0=0.011212, I0=0.009925, R0=0.126326, P0=0.004423, n=50000, h=0.005)
alpha = 0.28
S18, E18, I18, R18, P18, t18, dS18, dE18, dI18, dR18, dP18 = RK4(S0=0.845182, E0=0.014548, I0=0.014180, R0=0.126090, P0=0.004276, n=50000, h=0.005)
alpha = 0.32
S19, E19, I19, R19, P19, t19, dS19, dE19, dI19, dR19, dP19 = RK4(S0=0.848053, E0=0.008059, I0=0.014062, R0=0.129826, P0=0.005016, n=50000, h=0.005)
#print final positions and use these as starting conditions until we stay on a limit cycle
'''
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S10[-1], E10[-1], I10[-1], R10[-1], P10[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S11[-1], E11[-1], I11[-1], R11[-1], P11[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S12[-1], E12[-1], I12[-1], R12[-1], P12[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S13[-1], E13[-1], I13[-1], R13[-1], P13[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S14[-1], E14[-1], I14[-1], R14[-1], P14[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S15[-1], E15[-1], I15[-1], R15[-1], P15[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S16[-1], E16[-1], I16[-1], R16[-1], P16[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S17[-1], E17[-1], I17[-1], R17[-1], P17[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S18[-1], E18[-1], I18[-1], R18[-1], P18[-1]))
print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S19[-1], E19[-1], I19[-1], R19[-1], P19[-1]))
'''
#determine mins and maxes to use in bifuration diagram
alpha_values = [0.20, 0.25, 0.30, 0.22, 0.27, 0.34, 0.24, 0.19, 0.28, 0.32]
maxes_alpha = [np.max(I10), np.max(I11), np.max(I12), np.max(I13), np.max(I14), np.max(I15), np.max(I16), np.max(I17), np.max(I18), np.max(I19)]
mins_alpha = [np.min(I10), np.min(I11), np.min(I12), np.min(I13), np.min(I14), np.min(I15), np.min(I16), np.min(I17), np.min(I18), np.min(I19)]
beta0 = 50
alpha = 0.22
S89, E89, I89, R89, P89, t89, dS89, dE89, dI89, dR89, dP89 = RK4(S0=0.846612, E0=0.0115, I0=0.01125, R0=0.130638, P0=0.004150, n=300000, h=0.01)
S91, E91, I91, R91, P91, t91, dS91, dE91, dI91, dR91, dP89 = RK4(S0=0.813562, E0=0.027, I0=0.0288, R0=0.130638, P0=0.004150, n=300000, h=0.01)
alpha = 0.12
Ss, Es, Is, Rs, Ps, ts, dSs, dEs, dIs, dRs, dPs = RK4(S0=0.813562, E0=0.027, I0=0.0288, R0=0.130638, P0=0.004150, n=300000, h=0.01)
plt.plot(E89, I89, color='b')
plt.plot(E91, I91, color='r')
plt.plot(E13, I13, color='k')
plt.plot(0.0115, 0.01125, 'ko')
plt.plot(0.027, 0.0288, 'ko')
plt.ylabel('I')
plt.xlabel('E')
plt.title(r'$\alpha=0.220$')
plt.show()
plt.plot(Es, Is, color='b')
plt.plot(0.027, 0.0288, 'ko')
plt.ylabel('I')
plt.xlabel('E')
plt.title(r'$\alpha=0.120$')
plt.show()
plt.plot(E17, I17, color='m', label=r'$\alpha=0.195$')
plt.plot(E10, I10, color='b', label=r'$\alpha=0.200$')
plt.plot(E13, I13, color='r', label=r'$\alpha=0.220$')
plt.plot(E16, I16, color='k', label=r'$\alpha=0.240$')
plt.plot(E11, I11, color='g', label=r'$\alpha=0.250$')
plt.ylabel('I')
plt.xlabel('E')
plt.legend()
plt.show()
plt.plot(E14, I14, color='b', label=r'$\alpha=0.270$')
plt.plot(E18, I18, color='r', label=r'$\alpha=0.280$')
plt.plot(E12, I12, color='k', label=r'$\alpha=0.300$')
plt.plot(E19, I19, color='g', label=r'$\alpha=0.320$')
plt.plot(E15, I15, color='m', label=r'$\alpha=0.340$')
plt.ylabel('I')
plt.xlabel('E')
plt.legend()
plt.show()
def SEIRP_deriv(x): #x = [S,E,I,R,P]
dS = b*(N-x[0]) - (beta(x[4])*x[0]*x[2])/N
dE = beta(x[4])*x[0]*x[2]/N - (alpha + b)*x[1]
dI = alpha*x[1] - (gamma + b)*x[2]
dR = gamma*x[2] - b*x[3]
dP = e*gamma*x[2] - lamb*x[4]
return np.array([dS, dE, dI, dR, dP])
alpha=0.33 #fixed alpha
method, options = 'hybr', {'maxfev': 100000} #choose method for finding roots
start_point = (0.8567, 0.0055, 0.015, 0.1285, 0.0046) #define starting conditions
betas = np.linspace(0, 100, 10000) #beta0 will be varied
#prepare lists to save the solutions
S_roots_beta = []
E_roots_beta = []
I_roots_beta = []
R_roots_beta = []
P_roots_beta = []
for beta0 in betas: #walk over all beta0-values
sol = root(SEIRP_deriv, start_point, method=method, options=options) #find the roots
if sol.success: #if a solution was found, add it to the lists
S_roots_beta.append(sol.x[0])
E_roots_beta.append(sol.x[1])
I_roots_beta.append(sol.x[2])
R_roots_beta.append(sol.x[3])
P_roots_beta.append(sol.x[4])
else:
print('solution not found')
plt.plot(betas[:4866], I_roots_beta[:4866], color='k')
plt.plot(betas[4867:], I_roots_beta[4867:], 'k--') #The location of the Hopf-point, will be determined later on in the stability analysis
plt.plot(beta0_values, maxes_beta, 'ko')
plt.plot(beta0_values, mins_beta, 'ko')
plt.xlim(0, 100)
plt.ylim(-0.002, 0.025)
plt.xlabel(r'$\beta_0$')
plt.title(r'Bifurcation diagram for $\beta_0$')
plt.ylabel('I')
plt.plot(betas, np.zeros(betas.size), 'k--') #trivial solution (see paper)
plt.show()
plt.plot(betas, I_roots_beta, color='k')
plt.xlim(0, 1)
plt.ylim(-0.002, 0.004)
plt.xlabel(r'$\beta_0$')
plt.ylabel('I')
plt.title(r'Zoom on the bifurcation diagram for $\beta_0$')
plt.plot(betas, np.zeros(betas.size), 'k--') #trivial solution (see paper)
plt.show()
beta0 = 50 #keep beta0 fixed
method, options = 'hybr', {'maxfev': 100000} #choose method for finding roots
start_point = (0.8000, 0.0050, 0.0050, 0.1850, 0.0050) #define starting conditions
alphas = np.linspace(0, 0.5, 1000) #define range of alpha-values
#prepare lists to save the solutions
S_roots_alpha = []
E_roots_alpha = []
I_roots_alpha = []
R_roots_alpha = []
P_roots_alpha = []
for alpha in alphas: #walk over all alpha-values
sol = root(SEIRP_deriv, start_point, method=method, options=options) #find the roots
if sol.success: #if a solution was found, add it to the lists
S_roots_alpha.append(sol.x[0])
E_roots_alpha.append(sol.x[1])
I_roots_alpha.append(sol.x[2])
R_roots_alpha.append(sol.x[3])
P_roots_alpha.append(sol.x[4])
else:
print('solution not found')
plt.plot(alphas[:380], I_roots_alpha[:380], color='k')
plt.plot(alphas[381:686], I_roots_alpha[381:686], 'k--') #The location of the first Hopf-point, will be determined later on in the stability analysis
plt.plot(alphas[687:], I_roots_alpha[687:], color='k') #The location of the second Hopf-point, will be determined later on in the stability analysis
plt.plot(alpha_values, maxes_alpha, 'ko')
plt.plot(alpha_values, mins_alpha, 'ko')
plt.xlim(0, 0.5)
plt.ylim(-0.002, 0.017)
plt.xlabel(r'$\alpha$')
plt.title(r'Bifurcation diagram for $\alpha$')
plt.ylabel('I')
plt.plot(alphas, np.zeros(alphas.size), 'k--') #trivial solution (see paper)
plt.show()
import sympy as sym
#define SEIRP-model in sympy-format
def SEIRP_sym(N, b, alpha, gamma, e, lamb, beta0, mu, kappa, S, E, I, R, P):
eq1 = b*(N-S) - ((beta0*(1-mu)*sym.exp(-kappa*P/N))*S*I)/N
eq2 = (beta0*(1-mu)*sym.exp(-kappa*P/N))*S*I/N - (alpha + b)*E
eq3 = alpha*E - (gamma + b)*I
eq4 = gamma*I - b*R
eq5 = e*gamma*I - lamb*P
return sym.Matrix([eq1, eq2, eq3, eq4, eq5])
#define symbols
S_s = sym.Symbol('S')
E_s = sym.Symbol('E')
I_s = sym.Symbol('I')
R_s = sym.Symbol('R')
P_s = sym.Symbol('P')
N_s = sym.Symbol('N')
b_s = sym.Symbol('b')
alpha_s = sym.Symbol('alpha')
gamma_s = sym.Symbol('gamma')
e_s = sym.Symbol('e')
lamb_s = sym.Symbol('lambda')
beta0_s = sym.Symbol('beta_0')
mu_s = sym.Symbol('mu')
kappa_s = sym.Symbol('kappa')
Y = sym.Matrix([S_s, E_s, I_s, R_s, P_s])
#calculate Jacobian
SEIRP_sym(N_s, b_s, alpha_s, gamma_s, e_s, lamb_s, beta0_s, mu_s, kappa_s, S_s, E_s, I_s, R_s, P_s).jacobian(Y)
#numpy-implementation of the Jacobian (easier to calculate eigenvalues)
def jacobian(S, E, I, R, P):
j = np.array([[-b-beta(P)*I/N, 0, -beta(P)*S/N, 0, I*S*beta(P)*kappa/N**2],[beta(P)*I/N, -alpha-b, beta(P)*S/N, 0, -I*S*beta(P)*kappa/N**2],[0, alpha, -gamma-b, 0, 0],[0, 0, gamma, -b, 0],[0, 0, e*gamma, 0, -lamb]])
return j
alpha = 0.33 #fixed alpha
eigenvalues = np.zeros([np.size(betas), 5], dtype=float) #prepare an array to store the real parts of the eigenvalues
for i, beta0 in enumerate(betas): #loop over beta0-values
j = jacobian(N, 0, 0, 0, 0) #define jacobian for trivial solution
eigs, _ = eig(j) #calculate eigenvalues with scipy (imaginary part is ignored because dtype of array is float)
eigs = np.array(eigs) #make numpy array from tuple of eigenvalues
eigenvalues[i] = np.real(eigs) #save in array
plt.plot(betas, eigenvalues[:, 0], 'b')
plt.plot(betas, eigenvalues[:, 1], 'g')
plt.plot(betas, eigenvalues[:, 2], 'r')
plt.plot(betas, eigenvalues[:, 3], 'm')
plt.plot(betas, eigenvalues[:, 4], 'k')
plt.xlabel(r'$\beta_0$')
plt.ylabel('eigenvalues')
plt.title(r'Eigenvalues of the Jacobian when varrying $\beta_0$')
plt.show()
plt.plot(betas[0:100], eigenvalues[:, 0][0:100], 'b')
plt.plot(betas[0:100], eigenvalues[:, 1][0:100], 'g')
plt.plot(betas[0:100], eigenvalues[:, 2][0:100], 'r')
plt.plot(betas[0:100], eigenvalues[:, 3][0:100], 'm')
plt.plot(betas[0:100], eigenvalues[:, 4][0:100], 'k')
plt.xlabel(r'$\beta_0$')
plt.ylabel('eigenvalues')
plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\beta_0$')
plt.show()
print('The black eigenvalues turns from negative to positive between beta0=0.25 and beta0=0.26 --> trivial solution is stable for beta0<=0.25')
print('beta0={} --> black eigenvalue={}'.format(betas[25], eigenvalues[25, 4]))
print('beta0={} --> black eigenvalue={}'.format(betas[26], eigenvalues[26, 4]))
The black eigenvalues turns from negative to positive between beta0=0.25 and beta0=0.26 --> trivial solution is stable for beta0<=0.25
beta0=0.25002500250025006 --> black eigenvalue=-0.002852849840348559
beta0=0.26002600260026004 --> black eigenvalue=0.0023992340078670127
beta0 = 50 #fixed beta0
eigenvalues = np.zeros([np.size(alphas), 5], dtype=float) #prepare an array to store the real parts of the eigenvalues
for i, alpha in enumerate(alphas): #loop over alpha-values
j = jacobian(N, 0, 0, 0, 0) #define jacobian for trivial solution
eigs, _ = eig(j) #calculate eigenvalues with scipy (imaginary part is ignored because dtype of array is float)
eigs = np.array(eigs) #make numpy array from tuple of eigenvalues
eigenvalues[i] = np.real(eigs) #save in array
plt.plot(alphas, eigenvalues[:, 0], 'b')
plt.plot(alphas, eigenvalues[:, 1], 'g')
plt.plot(alphas, eigenvalues[:, 2], 'r')
plt.plot(alphas, eigenvalues[:, 3], 'm')
plt.plot(alphas, eigenvalues[:, 4], 'k')
plt.xlabel(r'$\alpha$')
plt.ylabel('eigenvalues')
plt.title(r'Eigenvalues of the Jacobian when varrying $\alpha$')
plt.show()
plt.plot(alphas[0:100], eigenvalues[:, 0][0:100], 'b')
plt.plot(alphas[0:100], eigenvalues[:, 1][0:100], 'g')
plt.plot(alphas[0:100], eigenvalues[:, 2][0:100], 'r')
plt.plot(alphas[0:100], eigenvalues[:, 3][0:100], 'm')
plt.plot(alphas[0:100], eigenvalues[:, 4][0:100], 'k')
plt.xlabel(r'$\alpha$')
plt.ylabel('eigenvalues')
plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\alpha$')
plt.show()
print('There is always a positive eigenvalue --> trivial solution is unstable for all alphas')
There is always a positive eigenvalue --> trivial solution is unstable for all alphas
alpha = 0.33 #fixed alpha
eigenvalues = np.zeros([np.size(betas), 5], dtype=float) #prepare array to store real parts of eigenvalues
for i, beta0 in enumerate(betas): #loop over beta0-values and make indices available
#use index to select a certain fixed point that was calulated earlier (see part on bifurcation diagrams)
SS = S_roots_beta[i]
EE = E_roots_beta[i]
II = I_roots_beta[i]
RR = R_roots_beta[i]
PP = P_roots_beta[i]
j = jacobian(SS, EE, II, RR, PP) #plug fixed point into jacobian
eigs, _ = eig(j) #calculate eigenvalues (imaginary part is ignored because dtype of array is float)
eigs = np.array(eigs) #make numpy array of tuple
eigenvalues[i] = np.real(eigs) #store in array
plt.plot(betas, eigenvalues[:, 0], 'b')
plt.plot(betas, eigenvalues[:, 1], 'g')
plt.plot(betas, eigenvalues[:, 2], 'r')
plt.plot(betas, eigenvalues[:, 3], 'm')
plt.plot(betas, eigenvalues[:, 4], 'k')
plt.xlabel(r'$\beta_0$')
plt.ylabel('eigenvalues')
plt.title(r'Eigenvalues of the Jacobian when varrying $\beta_0$')
plt.show()
plt.plot(betas, eigenvalues[:, 0], 'k', label='eigenvalue')
plt.plot(betas, eigenvalues[:, 1], 'k')
plt.plot(betas, eigenvalues[:, 2], 'k')
plt.plot(betas, eigenvalues[:, 3], 'k')
plt.plot(betas, eigenvalues[:, 4], 'k')
plt.xlim(40, 100)
plt.ylim(-0.01, 0.01)
plt.hlines(0, 40, 100, label=r'$y=0$', color='r')
plt.xlabel(r'$\beta_0$')
plt.ylabel('eigenvalues')
plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\beta_0$')
plt.legend()
plt.show()
print('eigenvalue = {}'.format(eigenvalues[:, 2][4866]))
print('eigenvalue = {}'.format(eigenvalues[:, 2][4867]))
print('beta_0 = {}'.format(betas[4866]))
eigenvalue = -5.412377913002753e-07
eigenvalue = 9.008643967363394e-07
beta_0 = 48.66486648664867
beta0=50 #fixed beta0
eigenvalues = np.zeros([np.size(alphas), 5], dtype=float) #prepare array to store real parts of eigenvalues
for i, alpha in enumerate(alphas): #loop over alpha-values and make indices available
#use index to select a certain fixed point that was calulated earlier (see part on bifurcation diagrams)
SS = S_roots_alpha[i]
EE = E_roots_alpha[i]
II = I_roots_alpha[i]
RR = R_roots_alpha[i]
PP = P_roots_alpha[i]
j = jacobian(SS, EE, II, RR, PP) #plug fixed point into jacobian
eigs, _ = eig(j) #calculate eigenvalues (imaginary part is ignored because dtype of array is float)
eigs = np.array(eigs) #make numpy array of tuple
eigenvalues[i] = np.real(eigs) #store in array
plt.plot(alphas, eigenvalues[:, 0], 'b')
plt.plot(alphas, eigenvalues[:, 1], 'g')
plt.plot(alphas, eigenvalues[:, 2], 'r')
plt.plot(alphas, eigenvalues[:, 3], 'm')
plt.plot(alphas, eigenvalues[:, 4], 'k')
plt.xlabel(r'$\alpha$')
plt.ylabel('eigenvalues')
plt.title(r'Eigenvalues of the Jacobian when varrying $\alpha$')
plt.show()
plt.plot(alphas, eigenvalues[:, 0], 'k', label='eigenvalue')
plt.plot(alphas, eigenvalues[:, 1], 'k')
plt.plot(alphas, eigenvalues[:, 2], 'k')
plt.plot(alphas, eigenvalues[:, 3], 'k')
plt.plot(alphas, eigenvalues[:, 4], 'k')
plt.xlim(0.18, 0.36)
plt.ylim(-0.001, 0.001)
plt.hlines(0, 0.18, 0.36, color='r', label=r'$y=0$')
plt.legend()
plt.xlabel(r'$\alpha$')
plt.ylabel('eigenvalues')
plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\alpha$')
plt.show()
print('eigenvalue = {}'.format(eigenvalues[:, 2][380]))
print('eigenvalue = {}'.format(eigenvalues[:, 2][381]))
print('alpha = {}'.format(alphas[380]))
print('\n')
print('eigenvalue = {}'.format(eigenvalues[:, 2][686]))
print('eigenvalue = {}'.format(eigenvalues[:, 2][687]))
print('alpha = {}'.format(alphas[686]))
eigenvalue = -1.0873164097474977e-05
eigenvalue = 1.341644201019107e-06
alpha = 0.19019019019019018
eigenvalue = 1.4305858619251266e-08
eigenvalue = -7.494583171244185e-06
alpha = 0.3433433433433433
alpha = 0.33
beta0 = 70
S_array_grid = np.linspace(0, 0.0500, 25)
E_array_grid = np.linspace(0, 0.0500, 25)
I_array_grid = np.linspace(0, 0.0500, 25)
R_array_grid = np.linspace(0, 0.0500, 25)
P_array_grid = np.linspace(0, 0.0500, 25)
B1 , C1 = np.meshgrid(E_array_grid, I_array_grid) #create grid
DA1, DB1, DC1, DD1, DE1 = SEIRP_deriv([0.8304, B1, C1, 0.1497, 0.0065]) #change cross-section here
M = (np.hypot(DB1, DC1)) # Norm of the growth rate
M[ M == 0] = 0.0001 # Avoid zero division errors
DB1 /= M # Normalize the arrows
DC1 /= M
#-------------------------------------------------------
fig = plt.figure(figsize=(17,11))
plt.ylabel('I')
plt.xlabel('E')
Q = plt.quiver(B1, C1, DB1, DC1, M, pivot='middle', cmap=plt.cm.jet)
plt.plot(E0,I0)
plt.xlim(0, 0.0500)
plt.ylim(0, 0.0500)
plt.title('Velocity vector field in the E,I-plane, and a limit cycle that intersects this plane')
plt.colorbar()
plt.show()
plt.plot(E0, I0)
for i in range(22000, 25000, 100):
plt.quiver(E0[i], I0[i], dE0[i], dI0[i], pivot='tail', color='r')
plt.xlabel('E')
plt.ylabel('I')
plt.title('Limit cycle with (E,I) velocity vectors')
plt.show()