Share Code Quickly

What is Harigami?

Harigami is a simple source code sharing service. You can share source code by an URL which is issued after submission.
You can also run your source code online such as Python, Ruby, JavaScript and so on.
For free.Use for code review.
Preferences
12345…326327
フェアリーテイル56~ kadai1
C 
#include <stdio.h>
#define NUMBER 5
#define FAILED -1
int search(int vc[], int key, int no){
    int i = 0;
    vc[no] = key;
    while(1){
        if(vc[i]==key){
            break;
        }
        i++;
    }
    return(i==no ? FAILED : i);
}

int main(void){
    int i,key,index;
    int vx[NUMBER+1];
    for (int i=0; i<NUMBER; i++) {
        printf("vx[%d]: ", i);
        scanf("%d", &vx[i]);
    }
    printf("探す値:");
    scanf("%d", &key);

    index = search(vx, key, NUMBER);
    
    if(index==FAILED){
        printf("探索失敗");
    }
    else {
        printf("%d is in %d\n", key, index);
    }
    return 0;
}
  • 1
  • 0
  • 2022年01月13日
anonymous No title
C 
#include <stdio.h>

int a[11] = {6, 7, 11, 12, 14, 18, 19, 21, 24, 28, 37};
int main(void){
    int key = 11;
    int l = -1;
    int r = sizeof(a)/sizeof(int);
    while (r-l > 1) {
        int m = l + (r - l)/2;
        if (a[m] > key) {
            r = m;
        }
        else {
            l = m;
        }
    }
    printf("%d", r);
    return 0;
}
  • 0
  • 0
  • 2022年01月13日
ごるだっく 3node_origin
Python 
# （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
  • 0
  • 0
  • 2022年01月11日
ごるだっく 3node_case2検証
Python 
# （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
  • 0
  • 0
  • 2022年01月11日
ごるだっく 3node_rockdown_封鎖期間変化
Python 
# （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/
  • 0
  • 0
  • 2022年01月11日
ごるだっく 3node_rockdown_閾値変化
Python 
# （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,
  • 0
  • 0
  • 2022年01月11日
ごるだっく 2node_安定性
Python 
# （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
  • 0
  • 0
  • 2022年01月11日
ごるだっく 2node_beta'変化
Python 
# （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
  • 0
  • 0
  • 2022年01月11日
ごるだっく SIIRmodel_感染制御
Python 
# （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()
  • 0
  • 0
  • 2022年01月11日
ごるだっく SIIRモデル,beta3値変化
Python 
# （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
  • 0
  • 0
  • 2022年01月11日
12345…326327