anonymous No title
No License Python
2021年10月12日
Copy Clone
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

# (2)時間変数tの導入
Time = 10             # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*Time                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,Time,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個

# (3)SIRモデル
# 3-1パラメータ
lamda = 0           # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
d=0                       #標的細胞の死亡率
beta=0.00001157                     #ウィルス粒子の侵入率
delta=3.412                     #感染細胞の死亡率
p=0.020099                         #ウィルス粒子を吐き出す数
c=3.381                        #ウィルス粒子の死亡率
# 3-2初期値
I_0 = 1                # 初期感染細胞
V_0 = 0.961            # 初期ウィルス粒子数
T_0 = 400000000        #初期標的細胞
# 3-3微分方程式
dTdt = lambda T, I, V, t : lamda-d*T-beta*T*V               # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIdt = lambda T, I, V, t : beta*T*V-delta*I       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dVdt = lambda T, I, V, t : p*I-c*V                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
T = np.empty(n)          # T[0],...,T[n-1] 要素数n個
I = np.empty(n)          # I[0],...,I[n-1] 要素数n個
V = np.empty(n)          # V[0],...,R[n-1] 要素数n個
# 3-5初期値代入
T[0] = T_0
I[0] = I_0
V[0] = V_0

# (4)数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個
  
  
  kT1 = h * dTdt( T[j] ,I[j] ,V[j] ,t[j] )
  kI1 = h * dIdt( T[j] ,I[j] ,V[j] ,t[j] )
  kV1 = h * dVdt( T[j] ,I[j] ,V[j] ,t[j] )

  kT2 = h * dTdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )
  kI2 = h * dIdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )
  kV2 = h * dVdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )

  kT3 = h * dTdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )
  kI3 = h * dIdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )
  kV3 = h * dVdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )

  kT4 = h * dTdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )
  kI4 = h * dIdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )
  kV4 = h * dVdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )

  T[j+1] = T[j] + 1/6 * ( kT1 + 2*kT2 + 2*kT3 + kT4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  I[j+1] = I[j] + 1/6 * ( kI1 + 2*kI2 + 2*kI3 + kI4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  V[j+1] = V[j] + 1/6 * ( kV1 + 2*kV2 + 2*kV3 + kV4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]

# (5)結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
plt.plot(t, T, color = "green", label = "T:標的細胞", linewidth = 1.0)
plt.plot(t, I, color = "red", label = "I:感染細胞", linewidth = 1.0)
plt.plot(t, V, color= "blue", label = "V:ウィルス粒子数", linewidth = 1.0)
# グラフの見た目設定
#plt.title('SIRモデル RK4によるシミュレーション(m={},T={})'.format(m,T))  # グラフタイトル パラメータmとTの値表示
#plt.yticks(np.arange(0,N+0.1,N/10))    # y軸 目盛りの配分 0からN(=1000万)までを10等分 N/10(=100万)刻み Nを含めるためNをN+0.1としておく
#plt.gca().set_yticklabels(['{:.0f}%'.format(y/(N/100)) for y in plt.gca().get_yticks()])   # y軸目盛りを%表示に変更
plt.xlabel('時間')                              # 横軸ラベル
plt.ylabel('数(総数に対する割合)')      # 縦軸ラベル
#plt.grid(True)                                        # グリッド表示
#plt.legend()                                          # 凡例表示
# 設定反映しプロット描画
plt.yscale('log')
plt.show()
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

# (2)時間変数tの導入
Time = 10             # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*Time                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,Time,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個

# (3)SIRモデル
# 3-1パラメータ
lamda = 0           # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
d=0                       #標的細胞の死亡率
beta=0.00001157                     #ウィルス粒子の侵入率
delta=3.412                     #感染細胞の死亡率
p=0.020099                         #ウィルス粒子を吐き出す数
c=3.381                        #ウィルス粒子の死亡率
# 3-2初期値
I_0 = 1                # 初期感染細胞
V_0 = 0.961            # 初期ウィルス粒子数
T_0 = 400000000        #初期標的細胞
# 3-3微分方程式
dTdt = lambda T, I, V, t : lamda-d*T-beta*T*V               # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIdt = lambda T, I, V, t : beta*T*V-delta*I       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dVdt = lambda T, I, V, t : p*I-c*V                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
T = np.empty(n)          # T[0],...,T[n-1] 要素数n個
I = np.empty(n)          # I[0],...,I[n-1] 要素数n個
V = np.empty(n)          # V[0],...,R[n-1] 要素数n個
# 3-5初期値代入
T[0] = T_0
I[0] = I_0
V[0] = V_0

# (4)数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個
  
  
  kT1 = h * dTdt( T[j] ,I[j] ,V[j] ,t[j] )
  kI1 = h * dIdt( T[j] ,I[j] ,V[j] ,t[j] )
  kV1 = h * dVdt( T[j] ,I[j] ,V[j] ,t[j] )

  kT2 = h * dTdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )
  kI2 = h * dIdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )
  kV2 = h * dVdt( T[j] + kT1/2 ,I[j] + kI1/2 ,V[j] + kV1/2 ,t[j] + h/2 )

  kT3 = h * dTdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )
  kI3 = h * dIdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )
  kV3 = h * dVdt( T[j] + kT2/2 ,I[j] + kI2/2 ,V[j] + kV2/2, t[j] + h/2 )

  kT4 = h * dTdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )
  kI4 = h * dIdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )
  kV4 = h * dVdt( T[j] + kT3 ,I[j] + kI3 ,V[j] + kV3 ,t[j] + h )

  T[j+1] = T[j] + 1/6 * ( kT1 + 2*kT2 + 2*kT3 + kT4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  I[j+1] = I[j] + 1/6 * ( kI1 + 2*kI2 + 2*kI3 + kI4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  V[j+1] = V[j] + 1/6 * ( kV1 + 2*kV2 + 2*kV3 + kV4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]

# (5)結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
plt.plot(t, T, color = "green", label = "T:標的細胞", linewidth = 1.0)
plt.plot(t, I, color = "red", label = "I:感染細胞", linewidth = 1.0)
plt.plot(t, V, color= "blue", label = "V:ウィルス粒子数", linewidth = 1.0)
# グラフの見た目設定
#plt.title('SIRモデル RK4によるシミュレーション(m={},T={})'.format(m,T))  # グラフタイトル パラメータmとTの値表示
#plt.yticks(np.arange(0,N+0.1,N/10))    # y軸 目盛りの配分 0からN(=1000万)までを10等分 N/10(=100万)刻み Nを含めるためNをN+0.1としておく
#plt.gca().set_yticklabels(['{:.0f}%'.format(y/(N/100)) for y in plt.gca().get_yticks()])   # y軸目盛りを%表示に変更
plt.xlabel('時間')                              # 横軸ラベル
plt.ylabel('数(総数に対する割合)')      # 縦軸ラベル
#plt.grid(True)                                        # グリッド表示
#plt.legend()                                          # 凡例表示
# 設定反映しプロット描画
plt.yscale('log')
plt.show()
No one still commented. Please first comment.
Output