import numpy as np
from matplotlib import rcParamsimport matplotlib.pyplot as plt
from scipy. integrate import odeint#sirモデルの常微分方程式
def SIR_equation(v,t,β,γ):
S, I, R =v
dSdt=-βSI
dIdt=βSI-γI
dRdt=γI
return [dSdt,dIdt,dRdt]
for β in np.linspace(0.10,0.20,11):
for γ in np.linspace(0.10,0.20,11):#初期値設定
S0=.999999928
I0=.000000071
R0=0.0
#係数決定#0 ≦t ≦100のあいだを1000分割
time_list=np.linspace(1,150,150)var_list=odeint( SIR_equation, [S0,I0,R0], time_list, args=(β,γ)
)
#微分方程式をここまでで解いた
#var_listにはS、I、Rのそれぞれの時間に対する結果が配列で入っているので、
#それらを分けて配列(_1ist)に収納
S_list = var_list [:,0]
I_list = var_list [:,1]
R_list = var_list [:,2]
tot_list = S_list + I_list + R_list#精題にプロットできるよう色々工夫してみる
plt.xlabel("time") plt.xlim(1,150) plt.ylim(0,0.0004) plt.ylabel("population ratio") plt.fill_between(time_list, tot_list,label ="population",alpha = 0.5) #plt.fill_between(time_list, S_list,label = "Susceptible",alpha = 0.5) #plt.fill_between(time_list, R_list,label = "Recovered",alpha = 0.5) plt.fill_between(time_list, I_list,label="Infected",alpha = 0.5) plt.fill_between(time_list, I_list+R_list,label="Infected(cumulative)",alpha = 0.5) plt.legend(loc ='center right') plt.title(r"SIR model ($\beta=%.2f,\gamma=%.2f$)"%(β,γ)) plt.savefig(r"SIR_model_beta=%.2f_gamma=%.2f.png"%(β,γ)) plt.close()
あなたの回答
tips
プレビュー