昨日のブログではネットで集めた情報をアップした．Facebook にも転送したところ，そちらに2つほどコメントをいただいた．受け売りを反省し，自分でコードを書いてみた．
左がテレビでよく見るグラフ．右は学術論文から転載したものだが，I と示した山型グラフだけをふたつ示したものが左のグラフであった．
これらのグラフは上左の連立微分方程式の解である．時間 t の関数 S(t) は Susceptible 感染する可能性がある人数，I(t) はInfectious 感染していて感染性を持つ(他人にうつす)人数，R(t) は Recovered 回復して免疫を持つ人数である．ただし R には死者も含む．β，γ は感染率，回復率である．
右が Mathmatica プログラムで，R(t) は本質的ではないので除外した，時間の単位を「日」とし，回復(あるいは死亡)までの日数を10日と仮定して，γ=1/10 とした．β > γ が感染拡大の条件なので，このプログラムでは β = 16γ すなわち 1/β=10日/16=15時間としてある．初期条件 I(0)=0.001 は，人口1000人なら，そのうちの一人にウイルスが取り付いたことを示す．これがシミュレーションのスタートである．
このモデルは閉じた系を対象とする．すなわちスタート時以外に感染者が外部から侵入することはない．ここでは β は(濃厚)接触率と言い換えてもよさそうだ ; すなわち，接触した時点で感染が起こる．昨日紹介したMASの解説ではβの大小を移動の速い・遅いで表現していたが，なるほどである．
感染者の割合の時間変化 I(t) のパラメータ β への依存性を上に示した． β = 16γ (1/β=10日/16=15時間) では最盛時には8割近くが感染する． β = 2γ (1/β=10日/2=5日) ではこれが2割程度になる...多いと言うべきか，少ないと言うべきか．感染割合が最大となる時期は，β が小さいほど後ろにずれる．ここでは省略するが β が小さいと t→∞ でも S(∞)=0 とならない．すなわち，感染とは無縁のハッピーな人たちが存在しうる．
Mathematica はこの近くの大学ではサイトライセンス契約をしているので，教職員・学生は自由に使うことができる．すでに示したように，このシミュレーションのプログラムは 10 行，パラメータ設定・初期条件・グラフ表示を除けば，本質的な部分はたった２行である．マニュアル検索から始めたが，小一時間で結果が出た...もっともさらにブログを書くために，すでに１時間以上かかっているが．
学生さんたちに，シミュレーションなんて簡単だよ! と言いたい! 冒頭・左のグラフについて，テレビの解説者たちがいかに適当なことを言っているかも，すこしでも自分で計算するとよくわかる!
