1）90日で感染者Iピーク(280万人)、200日後に感染者ゼロ、快復者R930万人、未感染者S70万人で収束。
2）320日で感染者Iピーク(50万人)、500日後に感染者ゼロ、快復者R530万人、未感染者S470万人で収束。
■「SIR感染モデル」のシミュレーションとそのグラフ化Gnuplotプログラム↓
#参考
#1) https://jbpress.ismedia.jp/articles/-/59947
#2) http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/index.html
#
# 実行(terminalから)：gnuplot -p -c "SIR.gp"
# 連立微分方程式
# dx/dt = Fx(x,y,z,t)
# dy/dt = Fy(x,y,z,t)
# dz/dt = Fz(x,y,z,t)
# SIRモデル
Fx(x,y,z,t)=-beta*x*y
Fy(x,y,z,t)= beta*x*y -gamma*y
Fz(x,y,z,t)= gamma*y
#パラメータ
# N:人口
# beta：一人当たり感染率
# gamma：一人当たり回復率（隔離率） 回復日数14日なら、回復率1/14=0.07
# x：感染可能者 S
# y；感染者 I
# ｚ:除去者(死亡、免疫獲得、隔離) R
N=1000*10000
gamma=0.07
A=5
P=2
beta= A*(P/100.0)/N #一人が平均A(人/日)に接触、接触ごとに確率P(%)で感染すると仮定。
# ステップ数と時間ステップ幅
T=700
dt=1
# 初期条件
y0=200
x0=N-y0
z0=0
# Runge-Kuttaステップの変数
Dx1=0.0
Dx2=0.0
Dx3=0.0
Dx4=0.0
Dy1=0.0
Dy2=0.0
Dy3=0.0
Dy4=0.0
Dz1=0.0
Dz2=0.0
Dz3=0.0
Dz4=0.0
M=int(N/10000.0)
TITLE='SIRモデルによる感染シミュレーション(人口' .M. '万人 1日の接触者数' .A. '(人/日) 1接触あたりの感染率' .P. '(%))'
set samples T
set title TITLE
set xrange [0:dt*(T-1)]
set grid
#set tmargin 5
set xlabel "日数"
set ylabel "人数"
# 別のファイルに書き出す
output="out.dat"
set table output
plot xold=x0,yold=y0,zold=z0,\
"+" using ($1+dt):\
(Dx1= Fx(xold, yold, zold,$1), \
Dy1 = Fy(xold,yold,zold,$1), \
Dz1 = Fz(xold,yold,zold,$1), \
Dx2 = Fx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dy2 = Fy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dz2 = Fz(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dx3 = Fx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dy3 = Fy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dz3 = Fz(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dx4 = Fx(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
Dy4 = Fy(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
Dz4 = Fz(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
xnew=xold+(dt/6.0)*(Dx1+2.0*Dx2+2.0*Dx3+Dx4), \
ynew=yold+(dt/6.0)*(Dy1+2.0*Dy2+2.0*Dy3+Dy4), \
znew=zold+(dt/6.0)*(Dz1+2.0*Dz2+2.0*Dz3+Dz4), \
xold=xnew, \
yold=ynew, \
zold=znew, \
xnew):(ynew):(znew):(0):(0) with xyerrorbars
unset table
set autoscale x
set autoscale y
set autoscale z
plot output using 2 w l title '感染可能者S', \
'' using 3 w l lc rgb "red" title '感染者I' , \
'' using 4 w l title '除去者R'
2）320日で感染者Iピーク(50万人)、500日後に感染者ゼロ、快復者R530万人、未感染者S470万人で収束。
■SIR感染モデル（人口一定都市(閉鎖都市)を仮定, つまりS(t)+I(t)+R(t)=N）
- 感受性保持者（Susceptible）
- 感染者（Infected）
- 免疫保持者（Recovered、あるいは隔離者 Removed）
- β > 0 は感染率、γ > 0 は回復率（隔離率）を表す（逆数 1/γ は平均感染期間を表す）
- 基本再生産数R0=βS0/γ
■「SIR感染モデル」のシミュレーションとそのグラフ化Gnuplotプログラム↓
#参考
#1) https://jbpress.ismedia.jp/articles/-/59947
#2) http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/index.html
#
# 準備
#Homebrew(https://brew.sh/index_ja)でgnuplotをインストール:brew install gnuplot
#
# 本ファイル：SIR.gp# 実行(terminalから)：gnuplot -p -c "SIR.gp"
# 連立微分方程式
# dx/dt = Fx(x,y,z,t)
# dy/dt = Fy(x,y,z,t)
# dz/dt = Fz(x,y,z,t)
# SIRモデル
Fx(x,y,z,t)=-beta*x*y
Fy(x,y,z,t)= beta*x*y -gamma*y
Fz(x,y,z,t)= gamma*y
#パラメータ
# N:人口
# beta：一人当たり感染率
# gamma：一人当たり回復率（隔離率） 回復日数14日なら、回復率1/14=0.07
# x：感染可能者 S
# y；感染者 I
# ｚ:除去者(死亡、免疫獲得、隔離) R
N=1000*10000
gamma=0.07
A=5
P=2
beta= A*(P/100.0)/N #一人が平均A(人/日)に接触、接触ごとに確率P(%)で感染すると仮定。
# ステップ数と時間ステップ幅
T=700
dt=1
# 初期条件
y0=200
x0=N-y0
z0=0
# Runge-Kuttaステップの変数
Dx1=0.0
Dx2=0.0
Dx3=0.0
Dx4=0.0
Dy1=0.0
Dy2=0.0
Dy3=0.0
Dy4=0.0
Dz1=0.0
Dz2=0.0
Dz3=0.0
Dz4=0.0
M=int(N/10000.0)
TITLE='SIRモデルによる感染シミュレーション(人口' .M. '万人 1日の接触者数' .A. '(人/日) 1接触あたりの感染率' .P. '(%))'
set samples T
set title TITLE
set xrange [0:dt*(T-1)]
set grid
#set tmargin 5
set xlabel "日数"
set ylabel "人数"
# 別のファイルに書き出す
output="out.dat"
set table output
plot xold=x0,yold=y0,zold=z0,\
"+" using ($1+dt):\
(Dx1= Fx(xold, yold, zold,$1), \
Dy1 = Fy(xold,yold,zold,$1), \
Dz1 = Fz(xold,yold,zold,$1), \
Dx2 = Fx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dy2 = Fy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dz2 = Fz(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
Dx3 = Fx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dy3 = Fy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dz3 = Fz(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
Dx4 = Fx(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
Dy4 = Fy(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
Dz4 = Fz(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
xnew=xold+(dt/6.0)*(Dx1+2.0*Dx2+2.0*Dx3+Dx4), \
ynew=yold+(dt/6.0)*(Dy1+2.0*Dy2+2.0*Dy3+Dy4), \
znew=zold+(dt/6.0)*(Dz1+2.0*Dz2+2.0*Dz3+Dz4), \
xold=xnew, \
yold=ynew, \
zold=znew, \
xnew):(ynew):(znew):(0):(0) with xyerrorbars
unset table
set autoscale x
set autoscale y
set autoscale z
plot output using 2 w l title '感染可能者S', \
'' using 3 w l lc rgb "red" title '感染者I' , \
'' using 4 w l title '除去者R'