import numpy as np import matplotlib.pyplot as plt from scipy import integrate N = 100000. I0 = 10. beta = 0.3 mu = 0.1 # Question 1 def new_cases_SIR(pop, t, beta, mu, N): S = pop[0] I = pop[1] R = pop[2] dS = -beta*S*I/N dI = beta*S*I/N-mu*I dR = mu*I return (dS, dI, dR) # Question 2 initial_pop = (N-I0, I0, 0) time = np.arange(0, 365, 1) SIR = integrate.odeint(new_cases_SIR, initial_pop, time, args = (beta, mu, N)) plt.plot(time, SIR) plt.legend(('S','I','R')) plt.savefig('SIR.png') plt.close() # Question 3 + 4 betas = np.arange(0.15, 0.31, 0.01) maxI = [] maxR = [] for beta in betas: SIR = integrate.odeint(new_cases_SIR, initial_pop, time, args = (beta, mu, N)) maxI.append(max(SIR.T[1])) maxR.append(max(SIR.T[2])) plt.plot(time, SIR.T[1]) plt.legend(betas) plt.savefig('I-beta.png') plt.close() # Question 5 plt.plot(betas/mu, maxI) plt.plot(betas/mu, maxR) plt.legend(('maxI', 'maxR')) plt.savefig('max.png') plt.close() # Question 6 l = 1e5 d = 0.1 a = 0.5 u = 5 k = 100 beta = 2.e-6 N = 1e6 V0 = 1 # Question 6 def new_cases_viral(pop, t, l, d, beta, a, k, u): S = pop[0] I = pop[1] V = pop[2] dS = l-d*S-beta*S*V dI = beta*S*V-a*I dV = k*I-u*V return (dS, dI, dV) initial_pop = (N, 0, N-V0) time = np.arange(0, 10, .1) SIV = integrate.odeint(new_cases_viral, initial_pop, time, args = (l, d, beta, a, k, u)) plt.plot(time, SIV) plt.legend(('S','I','V')) plt.savefig('SIV.png') plt.close()