|
|
コロナウィルス感染者数をSIRモデルでシュミレーションしてみました。
http://www.math.twcu.ac.jp/ogita/lec/sim_infection.pdf
https://rad-it21.com/サイエンス/michikoshi-shugo_20200331/
https://wagtail.cds.tohoku.ac.jp/coda/python/p-6-application-sup-ode-sir-model.html
https://ex-gram.com/infection_vs_vaccination2/
Susceptible:未感染者数
Infected :感染者数
Recovered :免疫保持者数、回復した人数
このS,I,Rを３元連立微分方程式としてルンゲクッタで解いてグラフ表示しています。
このプログラムではこれらパラメータをスライドバーで与えています。
Sを緑線
Iを赤線
Rを青線でグラフ表示しています。
縦軸を人数、横軸が日数です。
PUBLIC NUMERIC B,G
LOCATE VALUE NOWAIT(1) ,RANGE 10000 TO 10000000,AT 1000000 : N !'人口
LOCATE VALUE NOWAIT(2) ,RANGE 1 TO 10000,AT 10 : I !'初期感染者数
LOCATE VALUE NOWAIT(3) ,RANGE .01 TO 1,AT .2 : BETA !'感染率
LOCATE VALUE NOWAIT(4) ,RANGE 1 TO 30,AT 14 :GAMMA !'回復までの日数
DO
LOCATE VALUE NOWAIT(1) : N !'人口
LOCATE VALUE NOWAIT(2) : I !'初期感染者数
LOCATE VALUE NOWAIT(3) : BETA !'感染率
LOCATE VALUE NOWAIT(4) :GAMMA !'回復までの日数
LET N=INT(N)
LET I=INT(I)
LET R=0 !'初期回復者数
LET GAMMA=INT(GAMMA)
LET S=N-I !'初期未感染者数
LET T=0 !'経過日数
LET G=1/GAMMA !'回復率
LET B=BETA/N !'感染率
LET H=.5
LET DAY=180 !'期間(日)
SET WINDOW -15,DAY,-N/40,N*1.02
SET COLOR MIX(15) 0,0,0
DRAW GRID(30,N/5)
DO
LET K1=F1(T,S,I,R)
LET L1=F2(T,S,I,R)
LET M1=F3(T,S,I,R)
LET K2=F1(T+H/2,S+H/2*K1,I+H/2*L1,R+H/2*M1)
LET L2=F2(T+H/2,S+H/2*K1,I+H/2*L1,R+H/2*M1)
LET M2=F3(T+H/2,S+H/2*K1,I+H/2*L1,R+H/2*M1)
LET K3=F1(T+H/2,S+H/2*K2,I+H/2*L2,R+H/2*M2)
LET L3=F2(T+H/2,S+H/2*K2,I+H/2*L2,R+H/2*M2)
LET M3=F3(T+H/2,S+H/2*K2,I+H/2*L2,R+H/2*M2)
LET K4=F1(T+H,S+H*K3,I+H*L3,R+H*M3)
LET L4=F2(T+H,S+H*K3,I+H*L3,R+H*M3)
LET M4=F3(T+H,S+H*K3,I+H*L3,R+H*M3)
LET TT=T+H
LET SS=S+H*(K1+2*K2+2*K3+K4)/6
LET II=I+H*(L1+2*L2+2*L3+L4)/6
LET RR=R+H*(M1+2*M2+2*M3+M4)/6
SET LINE COLOR 3 ! 緑
PLOT LINES:TT,SS;T,S
SET LINE COLOR 4 ! 赤
PLOT LINES:TT,II;T,I
SET LINE COLOR 2 ! 青
PLOT LINES:TT,RR;T,R
LET S=SS
LET I=II
LET R=RR
LET T=TT
LOOP UNTIL T>=DAY
SET DRAW MODE EXPLICIT
WAIT DELAY .5
SET DRAW MODE HIDDEN
CLEAR
LOOP
END
EXTERNAL FUNCTION F1(T,S,I,R) !'dS/dt=f1(t,S,I,R)
LET F1=-B*S*I
END FUNCTION
EXTERNAL FUNCTION F2(T,S,I,R) !'dI/dt=f2(t,S,I,R)
LET F2=B*S*I-G*I
END FUNCTION
EXTERNAL FUNCTION F3(T,S,I,R) !'dR/dt=f3(t,S,I,R)
LET F3=G*I
END FUNCTION
|
|