# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
N2 = 100000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数（人）100人
Is2_0=0
R2_0 = 0                  # 初期回復者数（人）0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数（人）100人
Is3_0=0
R3_0 = 0                  # 初期回復者数（人）0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

  kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs22 =h * dIs2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR22 = h * dR2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS32 = h * dS3dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm32 =h * dIm3dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm1

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/10
beta23=(beta2+beta3)/10
beta31=(beta3+beta1)/10

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数（人）100人
Is2_0=0
R2_0 = 0                  # 初期回復者数（人）0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数（人）100人
Is3_0=0
R3_0 = 0                  # 初期回復者数（人）0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)



# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0


#c12をきる場合
c12=0
c23=1
c31=1
for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1]（I[j],R[j]も同じ） 要素数n個

  kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

  kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs22 =h * dIs2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR22 = h * dR2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS32 = h * dS3dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm32 =h * dIm3dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化
tha=np.arange
# （2）時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T（日）（250,150,150,700と値を変えてシミュレーションを行う）
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T （T=250のときはn=2500）
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個
tha=np.arange(0,30000,100)
lima=np.arange(0,300,1)
f=0 #flag
cnt=0
th=1000 #１日の感染者がこの値を超えたらlim日間ロックダウンする
lim=30 #ロックダウンの日数
# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数（人）100人
Is2_0=0
R2_0 = 0                  # 初期回復者数（人）0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数（人）100人
Is3_0=0
R3_0 = 0                  # 初期回復者数（人）0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)

peak=np.empty(300)
sumkekka=np.empty(300)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]= Is2_0

R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0

# （4）数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for i in range(300):
  #lim=lima[i]
  lim=lima[i]
  Im1_0 = 100                # 初期感染者数（人）100人
  Is1_0=0
  R1_0 = 0                  # 初期回復者数（人）0人
  S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  Im2_0 = 0                # 初期感染者数（人）100人
  Is2_0=0
  R2_0 = 0                  # 初期回復者数（人）0人
  S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  Im3_0 = 0                # 初期感染者数（人）100人
  Is3_0=0
  R3_0 = 0                  # 初期回復者数（人）0人
  S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1]（I[j],R[j]も同じ） 要素数n個

    if f==0 and j>10 and Is1[j]-Is1[j-10]>th:
      f=1
    
    if f==1:
      cnt=cnt+1
    
    if cnt>0 and cnt<lim:
      c23=1
      c12=0
      c31=0
    else:
      c23=1
      c12=1
      c31=1

    kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

    kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化
tha=np.arange
# （2）時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T（日）（250,150,150,700と値を変えてシミュレーションを行う）
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T （T=250のときはn=2500）
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個
tha=np.arange(0,30000,100)
lima=np.arange(0,300,1)
f=0 #flag
cnt=0
th=1000 #１日の感染者がこの値を超えたらlim日間ロックダウンする
lim=30 #ロックダウンの日数
# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数（人）100人
Is2_0=0
R2_0 = 0                  # 初期回復者数（人）0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数（人）100人
Is3_0=0
R3_0 = 0                  # 初期回復者数（人）0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)

peak=np.empty(300)
sumkekka=np.empty(300)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]= Is2_0

R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0

# （4）数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for i in range(300):
  #lim=lima[i]
  th=tha[i]
  Im1_0 = 100                # 初期感染者数（人）100人
  Is1_0=0
  R1_0 = 0                  # 初期回復者数（人）0人
  S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  Im2_0 = 0                # 初期感染者数（人）100人
  Is2_0=0
  R2_0 = 0                  # 初期回復者数（人）0人
  S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  Im3_0 = 0                # 初期感染者数（人）100人
  Is3_0=0
  R3_0 = 0                  # 初期回復者数（人）0人
  S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

  for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1]（I[j],R[j]も同じ） 要素数n個

    if f==0 and j>10 and Is1[j]-Is1[j-10]>th:
      f=1
    
    if f==1:
      cnt=cnt+1
    
    if cnt>0 and cnt<lim:
      c23=1
      c12=0
      c31=0
    else:
      c23=1
      c12=1
      c31=1

    kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

    kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化



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

kekka11=np.empty(n)
kekka12=np.empty(n)
kekka21=np.empty(n)
kekka22=np.empty(n)
kekka31=np.empty(n)
kekka32=np.empty(n)

sumkekka11=np.empty(n)
sumkekka12=np.empty(n)
sumkekka21=np.empty(n)
sumkekka22=np.empty(n)
sumkekka31=np.empty(n)
sumkekka32=np.empty(n)
# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # 街１の人口
N2 = 10000000               #街２の人口
m1 = 22                   #街１の接触数
m2 = 11                  #街２の接触数
p = 0.003               #5接触ごとに感染が生じる1日あたりの確率
d = 10                   # 感染者の回復平均日数（日）
nu1 = 0               #ワクチン接種率
nu2 = 0   
alpha = 0.01               #重症化率
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta_p=8*10**(-9)

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0=0
Is2_0=0
R2_0=0
S2_0=N2 - Im2_0 - Is2_0 - R2_0 
# 3-3微分方程式
dS1dt =  lambda S1,S2,Im1,Im2,Is1,Is2  ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : - beta1*S1*Im1 - nu1*S1 -beta_p*S1*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2, Im1,Im2,Is1,Is2, R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t :  beta1*S1*Im1 + beta_p*S1*Im2 - alpha*Im1 - gamma*Im1       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2, Im1,Im2,Is1,Is2, R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : alpha*Im1-gamma*Is1
dR1dt =  lambda S1,S2,Im1,Im2,Is1,Is2,  R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : - beta2*S2*Im2 - nu2*S2 -beta_p*S2*Im1             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : beta2*S2*Im2 + beta_p*S2*Im1 - alpha*Im2 - gamma*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t: alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS11 = h * dS1dt(   S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j] )
  kIm11 = h * dIm1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIs11 = h * dIs1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kR11 = h * dR1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])

  kS21 =  h * dS2dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIm21 = h * dIm2dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIs21 = h * dIs2dt( S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kR21 =  h * dR2dt(  S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])


  kS12 = h * dS1dt( S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm12 = h * dIm1dt( S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIs12=h * dIs1dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR12 = h * dR1dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)

  kS22 = h * dS2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIm22 = h * dIm2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIs22=h * dIs2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kR22 = h * dR2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)


  kS13 = h * dS1dt( S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm13 = h * dIm1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIs13 = h * dIs1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR13 = h * dR1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )

  kS23 = h * dS2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm23 = h * dIm2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIs23 = h * dIs2dt( S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR23 = h * dR2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )


  kS14 =  h * dS1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIm14 = h * dIm1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIs14 = h * dIs1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kR14 =  h * dR1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )


  kS24 =  h * dS2dt(  S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIm24 = h * dIm2dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIs24 = h * dIs2dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h  )
  kR24 =  h * dR2dt(  S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )


  S1[j+1]  = S1[j]  + 1/6 * ( kS11  + 2*kS12  + 2*kS13  + kS14 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im1[j+1] = Im1[j] + 1/6 * ( kIm11 + 2*kIm12 + 2*kIm13 + kIm14 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is1[j+1] = Is1[j] + 1/6 * ( kIs11 + 2*kIs12 + 2*kIs13 + kIs14 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R1[j+1]  = R1[j]  + 1/6 * ( kR11  + 2*kR12  + 2*kR13  + kR14 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  S2[j+1]  = S2[j]  + 1/6 * ( kS21  + 2*kS22  + 2*kS23  + kS24 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im2[j+1] = Im2[j] + 1/6 * ( kIm21 + 2*kIm22 + 2*kIm23 + kIm24 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is2[j+1] = Is2[j] + 1/6 * ( kIs21 + 2*kIs22 + 2*kIs23 + kIs24 ) 
  R2[j+1]  = R2[j]  + 1/6 * ( kR21  + 2*kR22  + 2*kR23  + kR24 ) 
  
  sum1[j+1]=sum1[j]+alpha*Im1[j]*h
  sum2[j+1]=sum2[j]+alpha*Im2[j]*h
  
plt.plot(t, S1, color = "blue", label = "S:未感染者", linewidth = 1.0)
plt.plot(t, Im1, color = "red", label = "Im:軽症者", linewidth = 1.0)
plt.plot(t, Is1, color = "green", label = "Is:重症者", linewidth = 1.0)
plt.plot(t, R1, color = "black", label = "R:回復者", linewidth = 1.0)
#plt.plot(t, R1, color= "blue", label = "R:免疫獲得者", li

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N1 = 10000000             # 街１の人口
N2 = 100000               #街２の人口
m1 = 10                   #街１の接触数
m2 = 10                  #街２の接触数
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0   
alpha = 0.01               #重症化率
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta_p=(beta1+beta2)/10

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im1_0 = 100                # 初期感染者数（人）100人
Is1_0=0
R1_0 = 0                  # 初期回復者数（人）0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

Im2_0=0
Is2_0=0
R2_0=0
S2_0=N2 - Im2_0 - Is2_0 - R2_0 
# 3-3微分方程式
dS1dt =  lambda S1,S2,Im1,Im2,Is1,Is2  ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : - beta1*S1*Im1 - nu1*S1 -beta_p*S1*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2, Im1,Im2,Is1,Is2, R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t :  beta1*S1*Im1 + beta_p*S1*Im2 - alpha*Im1 - gamma*Im1       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2, Im1,Im2,Is1,Is2, R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : alpha*Im1-gamma*Is1
dR1dt =  lambda S1,S2,Im1,Im2,Is1,Is2,  R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : - beta2*S2*Im2 - nu2*S2 -beta_p*S2*Im1             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t : beta2*S2*Im2 + beta_p*S2*Im1 - alpha*Im2 - gamma*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t: alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,Im1,Im2,Is1,Is2 ,R1,R2,alpha,beta1,beta2,beta_p,nu1,nu2, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS11 = h * dS1dt(   S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j] )
  kIm11 = h * dIm1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIs11 = h * dIs1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kR11 = h * dR1dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])

  kS21 =  h * dS2dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIm21 = h * dIm2dt(S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kIs21 = h * dIs2dt( S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])
  kR21 =  h * dR2dt(  S1[j],S2[j] ,Im1[j],Im2[j],Is1[j],Is2[j],R1[j],R2[j] ,alpha,beta1,beta2,beta_p,nu1,nu2 ,t[j])


  kS12 = h * dS1dt( S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm12 = h * dIm1dt( S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIs12=h * dIs1dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR12 = h * dR1dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)

  kS22 = h * dS2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIm22 = h * dIm2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kIs22=h * dIs2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)
  kR22 = h * dR2dt(S1[j] + kS11/2 ,S2[j]+kS21/2,Im1[j] + kIm11/2,Im2[j] + kIm21/2 ,Is1[j]+kIs11/2,Is2[j]+kIs21/2,R1[j] + kR11/2 ,R2[j]+kR21/2,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2)


  kS13 = h * dS1dt( S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm13 = h * dIm1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIs13 = h * dIs1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR13 = h * dR1dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22/2,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )

  kS23 = h * dS2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIm23 = h * dIm2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kIs23 = h * dIs2dt( S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )
  kR23 = h * dR2dt(S1[j] + kS12/2 ,S2[j]+kS22/2,Im1[j] + kIm12/2 ,Im2[j] + kIm22/2,Is1[j]+kIs12/2,Is2[j]+kIs22,R1[j] + kR12/2,R2[j]+kR22/2, alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h/2 )


  kS14 =  h * dS1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIm14 = h * dIm1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIs14 = h * dIs1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kR14 =  h * dR1dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )


  kS24 =  h * dS2dt(  S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIm24 = h * dIm2dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kIs24 = h * dIs2dt( S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )
  kR24 =  h * dR2dt(  S1[j] + kS13,S2[j]+kS23,Im1[j] + kIm13,Im2[j] + kIm23 ,Is1[j]+kIs13,Is2[j]+kIs23,R1[j] + kR13 ,R2[j]+kR23,alpha,beta1,beta2,beta_p,nu1,nu2,t[j] + h )


  S1[j+1]  = S1[j]  + 1/6 * ( kS11  + 2*kS12  + 2*kS13  + kS14 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im1[j+1] = Im1[j] + 1/6 * ( kIm11 + 2*kIm12 + 2*kIm13 + kIm14 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is1[j+1] = Is1[j] + 1/6 * ( kIs11 + 2*kIs12 + 2*kIs13 + kIs14 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R1[j+1]  = R1[j]  + 1/6 * ( kR11  + 2*kR12  + 2*kR13  + kR14 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  S2[j+1]  = S2[j]  + 1/6 * ( kS21  + 2*kS22  + 2*kS23  + kS24 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im2[j+1] = Im2[j] + 1/6 * ( kIm21 + 2*kIm22 + 2*kIm23 + kIm24 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is2[j+1] = Is2[j] + 1/6 * ( kIs21 + 2*kIs22 + 2*kIs23 + kIs24 ) 
  R2[j+1]  = R2[j]  + 1/6 * ( kR21  + 2*kR22  + 2*kR23  + kR24 ) 
  
  sum1[j+1]=sum1[j]+alpha*Im1[j]*h
  sum2[j+1]=sum2[j]+alpha*Im2[j]*h
  


print(max(Is1))
print(max(Is2))
print(sum1[n-1])

print(sum2[n-1])
# （5）結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
#plt.plot(t, S1, color = "green", label = "S:未感染者", linewidth = 1.0)
#plt.plot(t, Is1, color = "red", label = "Is:重症", linewidth = 1.0)
plt.plot(t, Is1, color = "blue", label = "beta=2.0*10^-8 max={}".format(max(Is1)), linewidth = 1.0)
plt.plot(t, Is2, color = "red", label = "beta=4.0*10^-8 max={}".format(max(Is2)), linewidth = 1.0)
plt.plot(t, Is1+Is2, color = "green", label = "beta=4.0*10^-8 max={}".format(max(Is2)), linewidth = 1.0)

#plt.plot(t

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

# （2）時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T（日）（250,150,150,700と値を変えてシミュレーションを行う）
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T （T=250のときはn=2500）
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個
ma=np.arange(1,100,1)
f=0
th=np.arange(0,400000,2000) #閾値
lim=np.arange(0,300,1)
# （3）SIRモデル
# 3-1パラメータ
N = 10000000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
m1 = 100
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0.005
nu3 = 0.01
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N           # 接触あたりの感染率
beta2 = m2*p / N           # 接触あたりの感染率
beta3 = m3*p / N           # 接触あたりの感染率

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im_0 = 100                # 初期感染者数（人）100人
Is_0=0
R_0 = 0                  # 初期回復者数（人）0人
S_0 = N - Im_0 - Im_0 - R_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dSdt = lambda S, Im,Is, R,alpha,beta,nu, t : - beta*S*Im - nu*S              # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dImdt = lambda S, Im,Is, R,alpha,beta,nu, t : beta*S*Im - alpha*Im - gamma*Im       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIsdt = lambda S, Im,Is, R,alpha,beta,nu, t : alpha*Im-gamma*Is
dRdt = lambda S, Im,Is, R,alpha,beta,nu, t : gamma*(Im+Is)+nu*S                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
maxIs=np.empty(99)
# 3-5初期値代入
S1[0] = S_0
Im1[0] = Im_0
Is1[0]=Is_0
R1[0] = R_0

S2[0] = S_0
Im2[0] = Im_0
Is2[0]=Is_0
R2[0] = R_0

S3[0] = S_0
Im3[0] = Im_0
Is3[0]=Is_0
R3[0] = R_0
sum1[0]=0
sumkekka=np.empty(300)
for i in range(300):
  c=0
  f=0
  S1[0] = S_0
  Im1[0] = Im_0
  Is1[0]=Is_0
  R1[0] = R_0

  S2[0] = S_0
  Im2[0] = Im_0
  Is2[0]=Is_0
  R2[0] = R_0

  S3[0] = S_0
  Im3[0] = Im_0
  Is3[0]=Is_0
  R3[0] = R_0
  sum1[0]=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個

    if j>10 and Is1[j]-Is1[j-10]>5000 and f==0:
      beta1=20*p/N
      f=1

    if f==1:
      c=c+1
    if c>lim[i]:
      beta1=100*p/N
    kS1 = h * dSdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
    kIm1 = h * dImdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
    kIs1 = h * dIsdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
    kR1 = h * dRdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )

    kS2 = h * dSdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
    kIm2 = h * dImdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
    kIs2=h * dIsdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
    kR2 = h * dRdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )

    kS3 = h * dSdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )
    kIm3 = h * dImdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
    kIs3 = h * dIsdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)    
    kR3 = h * dRdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )

    kS4 = h * dSdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
    kIm4 = h * dImdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
    kIs4 = h * dIsdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
    kR4 = h * dRdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

    S1[j+1] = S1[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
    Im1[j+1] = Im1[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
    Is1[j+1] = Is1[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
    R1[j+1] = R1[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
    sum1[j+1]=sum1[j]+alpha1*Im1[j]*h
      #2個目の計算
  sumkekka[i]=sum1[n-1]

# （5）結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
#plt.plot(t, S1, color = "green", label = "S:未感染者", linewidth = 1.0)
#plt.plot(t, Is1, color = "red", label = "Is:重症", linewidth = 1.0)
plt.plot(lim,sumkekka, color = "blue", linewidth = 1.0)

#plt.plot(t, R1, color= "blue", label = "R:免疫獲得者", linewidth = 1.0)
# グラフの見た目設定
plt.title('緊急事態宣言シミュレーション')  # グラフタイトル パラメータ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.show()

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N = 100000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
m1 = 10
m2 = 20                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 30
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0.005
nu3 = 0.01
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N           # 接触あたりの感染率
beta2 = m2*p / N           # 接触あたりの感染率
beta3 = m3*p / N           # 接触あたりの感染率

gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im_0 = 100                # 初期感染者数（人）100人
Is_0=0
R_0 = 0                  # 初期回復者数（人）0人
S_0 = N - Im_0 - Im_0 - R_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dSdt = lambda S, Im,Is, R,alpha,beta,nu, t : - beta*S*Im - nu*S              # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dImdt = lambda S, Im,Is, R,alpha,beta,nu, t : beta*S*Im - alpha*Im - gamma*Im       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIsdt = lambda S, Im,Is, R,alpha,beta,nu, t : alpha*Im-gamma*Is
dRdt = lambda S, Im,Is, R,alpha,beta,nu, t : gamma*(Im+Is)+nu*S                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
# 3-5初期値代入
S1[0] = S_0
Im1[0] = Im_0
Is1[0]=Is_0
R1[0] = R_0

S2[0] = S_0
Im2[0] = Im_0
Is2[0]=Is_0
R2[0] = R_0

S3[0] = S_0
Im3[0] = Im_0
Is3[0]=Is_0
R3[0] = R_0
sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS1 = h * dSdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIm1 = h * dImdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIs1 = h * dIsdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kR1 = h * dRdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )

  kS2 = h * dSdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  
  kR3 = h * dRdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIm4 = h * dImdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIs4 = h * dIsdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  kR4 = h * dRdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  S1[j+1] = S1[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im1[j+1] = Im1[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is1[j+1] = Is1[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R1[j+1] = R1[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum1[j+1]=sum1[j]+alpha1*Im1[j]*h
  
  #2個目の計算
  kS1 = h * dSdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta2,nu1 ,t[j] )
  kIm1 = h * dImdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta2,nu1 ,t[j] )
  kIs1 = h * dIsdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta2,nu1 ,t[j] )
  kR1 = h * dRdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta2,nu1 ,t[j] )

  kS2 = h * dSdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta2,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta2,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta2,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta2,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta2,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta2,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta2,nu1,t[j] + h/2)
  
  kR3 = h * dRdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta2,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta2,nu1,t[j] + h )
  kIm4 = h * dImdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta2,nu1,t[j] + h )
  kIs4 = h * dIsdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta2,nu1,t[j] + h )
  kR4 = h * dRdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta2,nu1,t[j] + h )

  S2[j+1] = S2[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im2[j+1] = Im2[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is2[j+1] = Is2[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R2[j+1] = R2[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum2[j+1]=sum2[j]+alpha1*Im2[j]*h
  
  #3つめの計算
  kS1 = h * dSdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta3,nu1 ,t[j] )
  kIm1 = h * dImdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta3,nu1 ,t[j] )
  kIs1 = h * dIsdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta3,nu1 ,t[j] )
  kR1 = h * dRdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta3,nu1 ,t[j] )

  kS2 = h * dSdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta3,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta3,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta3,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta3,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta3,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta3,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta3,nu1,t[j] + h/2)
  kR3 = h * dRdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta3,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta3,nu1,t[j] + h )
  kIm4 = h * dImdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta3,nu1,t[j] + h )
  kIs4 = h * dIsdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta3,nu1,t[j] + h )
  kR4 = h * dRdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta3,nu1,t[j] + h )

  S3[j+1] = S3[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im3[j+1] = Im3[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is3[j+1] = Is3[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R3[j+1] = R3[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum3[j+1]=sum3[j]+alpha1*Im3[j]*h
  

print(max(Is1))
print(max(Is2))
print(max(Is3))
print(sum1[n-1])

print(sum2[n-1])
print(sum3[n-1])
# （5）結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
#plt.plot(t, S1, color = "green", label = "S:未感染者", linewidth = 1.0)
#plt.plot(t, Is1, color = "red", label = "Is:重症", linewidth = 1.0)
plt.plot(t, Is1, color = "blue", label = "beta=2.0*10^-8 max={:.3f}".format(max(Is1)/N), linewidth = 1.0)
plt.plot(t, Is2, color = "red", label = "beta=4.0*10^-8 max={:.3f}".format(max(Is2)/N), linewidth = 1.0)
plt.plot(t, Is3, color = "green", label = "beta=6.0*10^-8 max={:.3f}".format(max(Is3)/N), linewidth = 1.0)

#plt.plot(t, R1, color= "blue", label = "R:免疫獲得者", linewidth = 1.0)
# グラフの見た目設定
plt.title('小さい街:βを変更:時間変化')  # グラフタイトル パラメータ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.show()


plt.plot(t, sum1, color = "blue", label = "beta=2.0*10^-8

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N = 100000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 100
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0.005
nu3 = 0.01
alpha1 = 0.01               #重症化率
alpha2 = 0.02
alpha3 = 0.1
beta1 = m1*p / N           # 接触あたりの感染率
beta2 = m2*p / N           # 接触あたりの感染率
beta3 = m3*p / N           # 接触あたりの感染率


gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im_0 = 100                # 初期感染者数（人）100人
Is_0=0
R_0 = 0                  # 初期回復者数（人）0人
S_0 = N - Im_0 - Im_0 - R_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dSdt = lambda S, Im,Is, R,alpha,beta,nu, t : - beta*S*Im - nu*S              # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dImdt = lambda S, Im,Is, R,alpha,beta,nu, t : beta*S*Im - alpha*Im - gamma*Im       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIsdt = lambda S, Im,Is, R,alpha,beta,nu, t : alpha*Im-gamma*Is
dRdt = lambda S, Im,Is, R,alpha,beta,nu, t : gamma*(Im+Is)+nu*S                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)

# 3-5初期値代入
S1[0] = S_0
Im1[0] = Im_0
Is1[0]=Is_0
R1[0] = R_0

S2[0] = S_0
Im2[0] = Im_0
Is2[0]=Is_0
R2[0] = R_0

S3[0] = S_0
Im3[0] = Im_0
Is3[0]=Is_0
R3[0] = R_0

sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS1 = h * dSdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIm1 = h * dImdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIs1 = h * dIsdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kR1 = h * dRdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )

  kS2 = h * dSdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  
  kR3 = h * dRdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIm4 = h * dImdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIs4 = h * dIsdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  kR4 = h * dRdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  S1[j+1] = S1[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im1[j+1] = Im1[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is1[j+1] = Is1[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R1[j+1] = R1[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum1[j+1]=sum1[j]+alpha1*Im1[j]*h
  #2個目の計算
  kS1 = h * dSdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha2,beta1,nu1 ,t[j] )
  kIm1 = h * dImdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha2,beta1,nu1 ,t[j] )
  kIs1 = h * dIsdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha2,beta1,nu1 ,t[j] )
  kR1 = h * dRdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha2,beta1,nu1 ,t[j] )

  kS2 = h * dSdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha2,beta1,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha2,beta1,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha2,beta1,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha2,beta1,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha2,beta1,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha2,beta1,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha2,beta1,nu1,t[j] + h/2)
  
  kR3 = h * dRdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha2,beta1,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha2,beta1,nu1,t[j] + h )
  kIm4 = h * dImdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha2,beta1,nu1,t[j] + h )
  kIs4 = h * dIsdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha2,beta1,nu1,t[j] + h )
  kR4 = h * dRdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha2,beta1,nu1,t[j] + h )

  S2[j+1] = S2[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im2[j+1] = Im2[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is2[j+1] = Is2[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R2[j+1] = R2[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum2[j+1]=sum2[j]+alpha2*Im2[j]*h
  
  #3つめの計算
  kS1 = h * dSdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha3,beta1,nu1 ,t[j] )
  kIm1 = h * dImdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha3,beta1,nu1 ,t[j] )
  kIs1 = h * dIsdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha3,beta1,nu1 ,t[j] )
  kR1 = h * dRdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha3,beta1,nu1 ,t[j] )
  kS2 = h * dSdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha3,beta1,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha3,beta1,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha3,beta1,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha3,beta1,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha3,beta1,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha3,beta1,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha3,beta1,nu1,t[j] + h/2)
  kR3 = h * dRdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha3,beta1,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha3,beta1,nu1,t[j] + h )
  kIm4 = h * dImdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha3,beta1,nu1,t[j] + h )
  kIs4 = h * dIsdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha3,beta1,nu1,t[j] + h )
  kR4 = h * dRdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha3,beta1,nu1,t[j] + h )

  S3[j+1] = S3[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im3[j+1] = Im3[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is3[j+1] = Is3[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R3[j+1] = R3[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum3[j+1]=sum3[j]+alpha3*Im3[j]*h
  
print(R1[n-1])
print(R2[n-1])
print(R3[n-1])

# （5）結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
#plt.plot(t, S1, color = "green", label = "S:未感染者", linewidth = 1.0)
#plt.plot(t, Is1, color = "red", label = "Is:重症", linewidth = 1.0)
plt.plot(t, Is1, color = "blue", label = "alpha=0.01 max={:.2f}".format(max(Is1)/N), linewidth = 1.0)
plt.plot(t, Is2, color = "red", label = "alpha=0.02 max={:.2f}".format(max(Is2)/N), linewidth = 1.0)
plt.plot(t, Is3, color = "green", label = "alpha=0.01 max={:.2f}".format(max(Is3)/N), linewidth = 1.0)

#plt.plot(t, R1, color= "blue", label = "R:免疫獲得者", linewidth = 1.0)
# グラフの見た目設定
plt.title('小さい街:αを変更:時間変化')  # グラフタイトル パラメータ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.show()

plt.plot(t, sum1, color = "blue", label = "alpha=0.01 sum={:.2f}".format(sum1[n-1]/N), linewidth = 1.0)
plt.plot(t, sum2, color = "

# （1）拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

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

# （3）SIRモデル
# 3-1パラメータ
N = 100000             # モデルエリアの人口（人）（東京都1400万人に匹敵するエリアを想定） N=S+I+R＝一定
m1 = 10
m2 = 20                  # 1日1人あたり接触する人数（人）（10,50,100,5と値を変えてシミュレーションを行う）
m3 = 30
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数（日）
nu1 = 0                 #ワクチン接種率
nu2 = 0.005
nu3 = 0.01
alpha1 = 0.01               #重症化率
alpha2 = 0.02
alpha3 = 0.03
beta1 = m1*p / N           # 接触あたりの感染率
beta2 = m2*p / N           # 接触あたりの感染率
beta3 = m3*p / N           # 接触あたりの感染率


gamma = 1/d              # 回復率（隔離率）
# 3-2初期値
Im_0 = 100                # 初期感染者数（人）100人
Is_0=0
R_0 = 0                  # 初期回復者数（人）0人
S_0 = N - Im_0 - Im_0 - R_0      # 初期未感染者数（人）S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dSdt = lambda S, Im,Is, R,alpha,beta,nu, t : - beta*S*Im - nu*S              # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dImdt = lambda S, Im,Is, R,alpha,beta,nu, t : beta*S*Im - alpha*Im - gamma*Im       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIsdt = lambda S, Im,Is, R,alpha,beta,nu, t : alpha*Im-gamma*Is
dRdt = lambda S, Im,Is, R,alpha,beta,nu, t : gamma*(Im+Is)+nu*S                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)
# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum1=np.empty(n)

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum3=np.empty(n)

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個
sum2=np.empty(n)
# 3-5初期値代入
S1[0] = S_0
Im1[0] = Im_0
Is1[0]=Is_0
R1[0] = R_0

S2[0] = S_0
Im2[0] = Im_0
Is2[0]=Is_0
R2[0] = R_0

S3[0] = S_0
Im3[0] = Im_0
Is3[0]=Is_0
R3[0] = R_0
sum1[0]=0
sum2[0]=0
sum3[0]=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個

  kS1 = h * dSdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIm1 = h * dImdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kIs1 = h * dIsdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )
  kR1 = h * dRdt( S1[j] ,Im1[j],Is1[j],R1[j] ,alpha1,beta1,nu1 ,t[j] )

  kS2 = h * dSdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIm2 = h * dImdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kIs2=h * dIsdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )
  kR2 = h * dRdt( S1[j] + kS1/2 ,Im1[j] + kIm1/2 ,Is1[j]+kIs1/2,R1[j] + kR1/2 ,alpha1,beta1,nu1,t[j] + h/2 )

  kS3 = h * dSdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )
  kIm3 = h * dImdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  kIs3 = h * dIsdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2)
  
  kR3 = h * dRdt( S1[j] + kS2/2 ,Im1[j] + kIm2/2 ,Is1[j]+kIs2/2,R1[j] + kR2/2, alpha1,beta1,nu1,t[j] + h/2 )

  kS4 = h * dSdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIm4 = h * dImdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )
  kIs4 = h * dIsdt(S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  kR4 = h * dRdt( S1[j] + kS3 ,Im1[j] + kIm3 ,Is1[j]+kIs3,R1[j] + kR3 ,alpha1,beta1,nu1,t[j] + h )

  S1[j+1] = S1[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im1[j+1] = Im1[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is1[j+1] = Is1[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R1[j+1] = R1[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum1[j+1]=sum1[j]+alpha1*Im1[j]*h
  #2個目の計算
  kS1 = h * dSdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta1,nu2 ,t[j] )
  kIm1 = h * dImdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta1,nu2 ,t[j] )
  kIs1 = h * dIsdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta1,nu2 ,t[j] )
  kR1 = h * dRdt( S2[j] ,Im2[j],Is2[j],R2[j] ,alpha1,beta1,nu2 ,t[j] )

  kS2 = h * dSdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta1,nu2,t[j] + h/2 )
  kIm2 = h * dImdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta1,nu2,t[j] + h/2 )
  kIs2=h * dIsdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta1,nu2,t[j] + h/2 )
  kR2 = h * dRdt( S2[j] + kS1/2 ,Im2[j] + kIm1/2 ,Is2[j]+kIs1/2,R2[j] + kR1/2 ,alpha1,beta1,nu2,t[j] + h/2 )

  kS3 = h * dSdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta1,nu2,t[j] + h/2 )
  kIm3 = h * dImdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta1,nu2,t[j] + h/2)
  kIs3 = h * dIsdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta1,nu2,t[j] + h/2)
  
  kR3 = h * dRdt( S2[j] + kS2/2 ,Im2[j] + kIm2/2 ,Is2[j]+kIs2/2,R2[j] + kR2/2, alpha1,beta1,nu2,t[j] + h/2 )

  kS4 = h * dSdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta1,nu2,t[j] + h )
  kIm4 = h * dImdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta1,nu2,t[j] + h )
  kIs4 = h * dIsdt(S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta1,nu2,t[j] + h )
  kR4 = h * dRdt( S2[j] + kS3 ,Im2[j] + kIm3 ,Is2[j]+kIs3,R2[j] + kR3 ,alpha1,beta1,nu2,t[j] + h )

  S2[j+1] = S2[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im2[j+1] = Im2[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is2[j+1] = Is2[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R2[j+1] = R2[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum2[j+1]=sum2[j]+alpha1*Im2[j]*h
  #3つめの計算
  kS1 = h * dSdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta1,nu3 ,t[j] )
  kIm1 = h * dImdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta1,nu3 ,t[j] )
  kIs1 = h * dIsdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta1,nu3 ,t[j] )
  kR1 = h * dRdt( S3[j] ,Im3[j],Is3[j],R3[j] ,alpha1,beta1,nu3 ,t[j] )

  kS2 = h * dSdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta1,nu3,t[j] + h/2 )
  kIm2 = h * dImdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta1,nu3,t[j] + h/2 )
  kIs2=h * dIsdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta1,nu3,t[j] + h/2 )
  kR2 = h * dRdt( S3[j] + kS1/2 ,Im3[j] + kIm1/2 ,Is3[j]+kIs1/2,R3[j] + kR1/2 ,alpha1,beta1,nu3,t[j] + h/2 )

  kS3 = h * dSdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta1,nu3,t[j] + h/2 )
  kIm3 = h * dImdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta1,nu3,t[j] + h/2)
  kIs3 = h * dIsdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta1,nu3,t[j] + h/2)
  kR3 = h * dRdt( S3[j] + kS2/2 ,Im3[j] + kIm2/2 ,Is3[j]+kIs2/2,R3[j] + kR2/2, alpha1,beta1,nu3,t[j] + h/2 )

  kS4 = h * dSdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta1,nu3,t[j] + h )
  kIm4 = h * dImdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta1,nu3,t[j] + h )
  kIs4 = h * dIsdt(S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta1,nu3,t[j] + h )
  kR4 = h * dRdt( S3[j] + kS3 ,Im3[j] + kIm3 ,Is3[j]+kIs3,R3[j] + kR3 ,alpha1,beta1,nu3,t[j] + h )

  S3[j+1] = S3[j] + 1/6 * ( kS1 + 2*kS2 + 2*kS3 + kS4 )   # 末項 j=n-2 -> S[j+1]=S[n-1]
  Im3[j+1] = Im3[j] + 1/6 * ( kIm1 + 2*kIm2 + 2*kIm3 + kIm4 )   # 末項 j=n-2 -> I[j+1]=I[n-1]
  Is3[j+1] = Is3[j] + 1/6 * ( kIs1 + 2*kIs2 + 2*kIs3 + kIs4 )   # 末項 j=n-2 -> I[j+1]=I[n-1] 
  R3[j+1] = R3[j] + 1/6 * ( kR1 + 2*kR2 + 2*kR3 + kR4 )   # 末項 j=n-2 -> R[j+1]=R[n-1]
  sum3[j+1]=sum3[j]+alpha1*Im3[j]*h

print(max(Is1))
print(max(Is2))
print(max(Is3))

print(sum1[n-1])
print(sum2[n-1])
print(sum3[n-1])

# （5）結果表示 データプロットによるグラフ表示
# 点(t,S),点(t,I),点(t,R) それぞれ要素数n個のプロット
#plt.plot(t, S1, color = "green", label = "S:未感染者", linewidth = 1.0)
#plt.plot(t, Is1, color = "red", label = "Is:重症", linewidth = 1.0)
plt.plot(t, Is1, color = "blue", label = "nu=0. max={:.3f}".format(max(Is1)/N), linewidth = 1.0)
plt.plot(t, Is2, color = "red", label = "nu=0.01 max={:.3f}".format(max(Is2)/N), linewidth = 1.0)
plt.plot(t, Is3, color = "green", label = "nu=0.02 max={:.3f}".format(max(Is3)/N), linewidth = 1.0)


#plt.plot(t, R1, color= "blue", label = "R:免疫獲得者", linewidth = 1.0)
# グラフの見た目設定
plt.title('小さい街:νを変更:時間変化')  # グラフタイトル パラメータ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.show()
#累積のほうのグラフ
plt.plot(t, sum1, color = "blue", label = "nu=0 sum={:.3f}".format(sum1[n-1]/N