sirsモデル　コロナに対するシミュレーションについて

こちらはsirモデルについて書いたものをsirsモデルに変えようとしています。
参考サイトを貼ります。
[リンク内容](https://rpubs.com/ktgrstsh/tokyor84 )
こちらを元に助けてもらいたいです。
sirは実際に動いたのでこれにこちらのサイトでいうrhoに関するものを付け加えていくと言った感じです。
初めて質問させてもらうので色々不備あると思いますが、よろしくお願い致します。
質問してください！

import
1from matplotlib import rcParams
2import matplotlib.pyplot as plt
3from scipy. integrate import odeint
4
5
6#sirモデルの常微分方程式
7
8def SIR_equation(v,t,β,γ,):
9    S, I, R =v
10    dSdt=-β*S*I
11    dIdt=β*S*I-γ*I
12    dRdt=γ*I
13    return [dSdt,dIdt,dRdt]
14for β in np.linspace(0.10,0.20,11):
15    for γ in np.linspace(0.10,0.20,11):
16        
17        
18#初期値設定
19          S0=.999999928
20          I0=.000000071
21          R0=0.0
22             
23#係数決定
24
25
26#0 ≦t ≦100のあいだを1000分割
27        time_list=np.linspace(1,150,150)
28
29        var_list=odeint(
30           SIR_equation,
31           [S0,I0,R0],
32           time_list,
33           args=(β,γ)
34)
35#微分方程式をここまでで解いた
36#var_listにはS、I、Rのそれぞれの時間に対する結果が配列で入っているので、
37#それらを分けて配列（_1ist）に収納
38        S_list = var_list [:,0] 
39        I_list = var_list [:,1] 
40        R_list = var_list [:,2] 
41        tot_list = S_list + I_list + R_list
42
43#精題にプロットできるよう色々工夫してみる
44       
45        plt.xlabel("time")
46        plt.xlim(0,100)
47        plt.ylim(0,1)
48        plt.ylabel("population ratio")
49        plt.fill_between(time_list, tot_list,label ="population",alpha = 0.5)
50        #plt.fill_between(time_list, S_list,label = "Susceptible",alpha = 0.5)
51        #plt.fill_between(time_list, R_list,label = "Recovered",alpha = 0.5)
52        plt.fill_between(time_list, I_list,label="Infected",alpha = 0.5)
53        plt.fill_between(time_list, I_list+R_list,label="Infected(cumulative)",alpha = 0.5)
54        
55
56        plt.legend(loc ='center right')
57        plt.title(r"SIR model ($\beta=%.2f,\gamma=%.2f$)"%(β,γ))
58        plt.savefig(r"SIR_model_beta=%.2f_gamma=%.2f.png"%(β,γ))
59        plt.close()
60
61コード

can110

2020/07/14 08:47

前回の質問と同じようですが質問文は修正できます。前回の質問を修正し、この質問は削除依頼をかけてください。

itta0602

2020/07/14 09:06

わかりました！canさんはこちらの質問わかりますでしょうか？

can110

2020/07/14 09:14

いえ。専門外なので分かりません。ひとつアドバイスするとすれば、専門外の人（理系高校生レベルくらい）に分かるように sir、sirsモデルおよびその違いについて、かいつまんで説明してもらうと回答得られやすくなると思います。

itta0602

2020/07/14 09:48

ありがとうございました

Penpen7

2020/07/14 10:04 編集

dSdt=-β*S*I dIdt=β*S*I-γ*I dRdt=γ*I にρR項を文字通り付け足せばいいだけだと思うのですが、違うんですか？

Penpen7

2020/07/14 10:07

def SIR_equation(v,t,β,γ,): は def SIR_equation(v,t,β,γ,ρ): としたかったのですかね？

行動規範の内容に同意します

回答1件

ベストアンサー

「文字通り」ρを加えればいいだけの話でしょう。
(各定数は参考URLに準拠, for文は簡単化のため削除)

python
1from matplotlib import rcParams
2import matplotlib.pyplot as plt
3from scipy. integrate import odeint
4import numpy as np
5
6
7def SIRS_equation(v, t, β, γ, ρ):
8    S, I, R = v
9    dSdt = -β*S*I+ρ*R
10    dIdt = β*S*I-γ*I
11    dRdt = γ*I-ρ*R
12    return [dSdt, dIdt, dRdt]
13
14ρ=0.01
15β=1.0
16γ=0.1
17S0 = .9999999
18I0 = 1-S0
19R0 = 0.0
20time_list = np.linspace(0, 150, 150)
21
22var_list = odeint(
23    SIRS_equation, [S0, I0, R0], time_list, args=(β, γ, ρ))
24S_list = var_list[:, 0]
25I_list = var_list[:, 1]
26R_list = var_list[:, 2]
27tot_list = S_list + I_list + R_list
28
29# 精題にプロットできるよう色々工夫してみる
30
31plt.xlabel("time")
32plt.xlim(0, 100)
33plt.ylim(0, 1)
34plt.ylabel("population ratio")
35plt.fill_between(time_list, tot_list, label="population", alpha=0.5)
36plt.fill_between(time_list, I_list, label="Infected", alpha=0.5)
37plt.fill_between(time_list, I_list+R_list, label="Infected(cumulative)", alpha=0.5)
38plt.legend(loc ='center right')
39plt.title(r"SIRS model ($\beta=%.2f,\gamma=%.2f, \rho=%.2f$)"%(β, γ, ρ))
40plt.savefig(r"SIRS_model_beta=%.2f_gamma=%.2f_rho=%.2f.png" %(β, γ, ρ))
41plt.close()