日本語プログラミング言語なでしこで数値計算（3） コロナウイルスのような感染症流行を表す微分方程式、SIRモデルをルンゲクッタ法で計算する。西浦博さんの「感染症疫学のためのデータ分析入門」を参考に。

さて1回目と2回目はこういうのをやってみた。

日本語プログラミング言語なでしこで数値計算（１）　まずは何はともあれ4段4次のルンゲクッタ法でローレンツ方程式を計算してみる。

日本語プログラミング言語なでしこで数値計算（2）　カオス（ローレンツ方程式）の次はフラクタルだということでマンデルブロ集合を描いてみる。

今回はせっかくルンゲクッタ法の計算ルーチンを作ったので、このコロナ禍で有名になったSIRモデル

https://ja.wikipedia.org/wiki/SIR%E3%83%A2%E3%83%87%E3%83%AB

を計算してみよう。

説明はプログラムにも書いていますが

# コロナウイルスのような感染症流行についての微分方程式
# SIRモデル
# dS/dt = -βSI
# dI/dt = βSI-γI
# dR/dt = γI
# を4段4次のルンゲクッタ法で計算する。
# S: まだ未感染で免疫がない感受性人口
# I: 感染している感染性人口
# R: 回復した治癒人口
# N:全人口=S+I+R (1に規格化）
# β：感染伝達率, γ：回復率

で、これは西浦博さんの「感染症疫学のためのデータ分析入門」を参考にした。

ではプログラムのリンクはこちら：

感染症流行のSIRモデルをルンゲクッタ法で計算

結果の例はこちら。S,I,Rの初期値やβ、γのパラメータをいろいろいじって遊べる。

リストはこんな感じ。

テキストデータでのソースコード：

# コロナウイルスのような感染症流行についての微分方程式

# SIRモデル

# dS/dt = -βSI

# dI/dt = βSI-γI

# dR/dt = γI

# を4段4次のルンゲクッタ法で計算する。

# S: まだ未感染で免疫がない感受性人口

# I: 感染している感染性人口

# R: 回復した治癒人口

# N:全人口=S+I+R (1に規格化）

# β：感染伝達率, γ：回復率

全描画クリア。

#　画面設定

画面幅＝描画中キャンバスの「width」をDOM属性取得。

画面高＝描画中キャンバスの「height」をDOM属性取得。 

「18px sans-serif」に描画フォント設定。

[80,20]に「SIRモデル(ルンゲクッタ法で計算)」を文字描画。

[80,40]に「S:赤, I:青, R:緑」を文字描画。

時間に0を代入する。

時間ステップに0.01を代入する。

最大時間に140を代入する。

x最大値に最大時間を代入する。

x最小値に-10を代入する。

y最大値に1.2を代入する。

y最小値に-0.2を代入する。

x最小値と0からx最大値と0まで黒色で1で直線プロット。

0とy最小値から0とy最大値まで黒色で1で直線プロット。

#初期値

S0に0.9999を代入する。

I0に0.0001を代入する。

R0に0を代入する。

ベータに0.3を代入する。

ガンマに0.1を代入する。

時間が最大時間以下の間

#ルンゲクッタ法のメイン計算

　KS1にS0とI0とR0の関数fSを代入する。

   KI1にS0とI0とR0の関数fIを代入する。

   KR1にS0とI0とR0の関数fRを代入する。

   KS2に(S0+時間ステップ×KS1/2)と(I0＋時間ステップ×KI1/2)と(R0＋時間ステップ×KR1/2)の関数fSを代入する。

   KI2に(S0+時間ステップ×KS1/2)と(I0＋時間ステップ×KI1/2)と(R0＋時間ステップ×KR1/2)の関数fIを代入する。

   KR2に(S0+時間ステップ×KS1/2)と(I0＋時間ステップ×KI1/2)と(R0＋時間ステップ×KR1/2)の関数fRを代入する。

   KS3に(S0+時間ステップ×KS2/2)と(I0＋時間ステップ×KI2/2)と(R0＋時間ステップ×KR2/2)の関数fSを代入する。

   KI3に(S0+時間ステップ×KS2/2)と(I0＋時間ステップ×KI2/2)と(R0＋時間ステップ×KR2/2)の関数fIを代入する。

   KR3に(S0+時間ステップ×KS2/2)と(I0＋時間ステップ×KI2/2)と(R0＋時間ステップ×KR2/2)の関数fRを代入する。

   KS4に(S0+時間ステップ×KS3)と(I0＋時間ステップ×KI3)と(R0＋時間ステップ×KR3)の関数fSを代入する。

   KI4に(S0+時間ステップ×KS3)と(I0＋時間ステップ×KI3)と(R0＋時間ステップ×KR3)の関数fIを代入する。

   KR4に(S0+時間ステップ×KS3)と(I0＋時間ステップ×KI3)と(R0＋時間ステップ×KR3)の関数fRを代入する。

   S1 にS0 + 時間ステップ×(KS1+2*KS2+2*KS3+KS4)/6を代入する。

   I1 にI0 + 時間ステップ×(KI1+2*KI2+2*KI3+KI4)/6を代入する。

   R1 にR0 + 時間ステップ×(KR1+2*KR2+2*KR3+KR4)/6を代入する。

   時間とS0から（時間+時間ステップ）とS1まで赤色で5で直線プロット。

   時間とI0から（時間+時間ステップ）とI1まで青色で5で直線プロット。

   時間とR0から（時間+時間ステップ）とR1まで緑色で5で直線プロット。

   S0=S1

   I0=I1

   R0=R1

　 時間を時間ステップだけ増やす。

ここまで。

#SIRモデル

●（SとIとRの）関数fSとは

　-ベータ×S×Iを戻すこと

ここまで

●（SとIとRの）関数fIとは

　ベータ×S×I - ガンマ×Iを戻すこと

ここまで

●（SとIとRの）関数fRとは

　ガンマ×Iを戻すこと

ここまで

●(xの）x座標変換とは

    それ=画面幅 * (x - x最小値）/(x最大値 - x最小値)

ここまで

●(yの）y座標変換とは

    それ=画面高 * (y - y最大値）/(y最小値 - y最大値)

ここまで

●(x1とy1からx2とy2までcでwで)直線プロットとは

　　cに線色設定

　　wに線太設定

　　[x1のx座標変換, y1のy座標変換]から[x2のx座標変換, y2のy座標変換]まで線描画

ここまで