« 4回目のコロナウイルスワクチン(2価のファイザーのワクチン、オミクロンBA.1対応)を打った。発熱はほぼなし。ただ打った腕を中心にあちこちが痛い、、、が1日で治った。 | トップページ | アンディ・ウォーホル・キョウト / ANDY WARHOL KYOTO@京都市京セラ美術館を観てきた。京都に2回来たことがあったのか!しかも徹子の部屋に出てた。三つのマリリン、最後の晩餐が見られます。あと坂本龍一さんが急に出てきたのにびっくりする。 »

2022年10月17日 (月)

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

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

日本語プログラミング言語なでしこで数値計算(1) まずは何はともあれ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の初期値やβ、γのパラメータをいろいろいじって遊べる。

Sir1

リストはこんな感じ。

Sir2

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

 


# コロナウイルスのような感染症流行についての微分方程式
# 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座標変換]まで線描画
ここまで

 

 

 


 

« 4回目のコロナウイルスワクチン(2価のファイザーのワクチン、オミクロンBA.1対応)を打った。発熱はほぼなし。ただ打った腕を中心にあちこちが痛い、、、が1日で治った。 | トップページ | アンディ・ウォーホル・キョウト / ANDY WARHOL KYOTO@京都市京セラ美術館を観てきた。京都に2回来たことがあったのか!しかも徹子の部屋に出てた。三つのマリリン、最後の晩餐が見られます。あと坂本龍一さんが急に出てきたのにびっくりする。 »

パソコン・インターネット」カテゴリの記事

学問・資格」カテゴリの記事

日記・コラム・つぶやき」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 4回目のコロナウイルスワクチン(2価のファイザーのワクチン、オミクロンBA.1対応)を打った。発熱はほぼなし。ただ打った腕を中心にあちこちが痛い、、、が1日で治った。 | トップページ | アンディ・ウォーホル・キョウト / ANDY WARHOL KYOTO@京都市京セラ美術館を観てきた。京都に2回来たことがあったのか!しかも徹子の部屋に出てた。三つのマリリン、最後の晩餐が見られます。あと坂本龍一さんが急に出てきたのにびっくりする。 »

最近の記事

最近のコメント

2022年12月
        1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 31
フォト
無料ブログはココログ