SIRモデルの改良 [数値解析]
SIRモデルの改良
ケルマックとマッケンドリンクによって１９２７年に提出されたSIRモデル
に出産率b、自然死亡率μを導入し、次のように修正する。
（２）の辺々を足し合わせると、
したがって、全人口
に関して、次の式を得る。
とすると、安定な平衡値として
が得られる。
問 次の微分方程式を解き、
となることを確かめよ。
【解】
両辺に
を掛けると、
両辺をtで積分すると、
μ>0だから、
となり、全人口の平衡値
に一致する。
（解答終）
そこで、はじめから全人口は
であると仮定し、
について考えればよい。
感染症の流行初期の段階では、S(t)≒Nだから、（４）の第２式は次のように線形化可能（マルサスのモデル）。
したがって、基本再生産数
と定義することができ、R₀>1ならば感染は拡大し、R₀<1ならば感染は自然消滅することになる。
（４）式で
になる(S,I)の点――このような点を平衡点という。力学系の話なので、詳しい話は、物理担当のddt³さんに聞くとよい(^^ゞ――を求めると、
という定常解が得られる。
前者の(N,0)は感染者がいない場合の定常解である。後者の定常解はR₀>1のときにのみ意味をもつ解で、これは集団の中に感染者が常に存在するエンデミックを表す定常解である。
――シンプルなSIRモデルでは、R₀の値にかかわらず
となるので、エンデミックな状態を再現できないが、修正した新たなモデルでは、出産やヒトの流入などで、集団の感受人口が常に供給されるので、エンデミックが起きる。――
エンデミックの状態を
とすると、
ネムネコは、感染症の流行モデルの研究者じゃなく、ただのネコなので、こんな抽象的な話をされても具体的なイメージがわかなくて、わかった気分になれない。
――あくまで”わかった”気分で、本当にここまでの内容を理解できているかどうかは、また、別の話！！――
だ・か・ら、
数値的であろうがどのような方法であろうが、（３）を解き、それをグラフにし、視覚化しないと！！
というわけで、例によって、出生率ｂ、自然死亡率μ、感染率β、隔離率（回復率）γに適当な値を入れ、オイラー法を用いて、連立微分方程式（４）を解いてみたにゃ。
――オイラー法ならば、誰でも表計算ソフトを使うことで簡単に計算できる！！。根性を出せば、手計算や、そろばん、タイガー計算機などを使って、この問題を数値的に解くことさえできる。だって、この計算で使うのは、小学校で習う加減乗除の四則演算のみで、ルートや三角関数といった難しい数学は一切不要なんだから――
出生率b=1、自然死亡率μ=0.01とすると、安定な人口の平衡値Nは、
そして、β=0.005、γ=0.1とすると、基本再生産数R₀は
t=0における感受性人口S(0)を９９、感染者人口I(0)=1とし、オイラー法を用いて（４）を解くと、次のような結果が得られる。
横軸にS、縦軸にIを取った位相図は次のとおり。
t=400におけるS(t)、I(t)の値は
で、（７）式を用いて計算したエンデミックな定常状態
を再現していることがわかるにゃ。
「オイラー法は精度が悪い」と、多くのヒトがオイラー法を酷評するけれど、この問題の場合、意外に正確に計算できているにゃ。
だ・か・ら、
定性的に、そして、ある程度、定量的に（３）の解の挙動を知りたないのならば、簡単なオイラー法で十分だケロ。
数値計算の講義を担当する（大学の）先生方などは、決まって、「オイラー法は精度が悪いから使っちゃダメ。４次のルンゲ・クッタを使え」と、オイラー法を悪者扱いするけれど、オイラー法はそうそう捨てたもんじゃないってこと。
ではあるが、
４次のルンゲ・クッタ法を用いて、計算領域をt=500まで広げ、より高精度に計算した結果は次のとおり。
この計算結果によると、少し振動しながら、エンデミック状態に近づいてゆくみたいだね。
このことは、位相図を見ると、よくわかるにゃ。
反時計回りの渦（？）を描きならが、平衡点
に収束するみたいだケロね。
そして、t=500くらいまで計算すると、先に求めたエンデミック状態の解(22.00, 7.09)に到達するのであった。
この図だと、うずまき具合がわかりづらいので、R₀=2.38の場合について、オイラー法で数値的に解いた結果を示すにゃ。
これだと、渦巻きを描きながら、パンデミックの定常解に近づくことがよくわかると思うにゃ。
記事の内容とは関係ないけれど、この動画も埋め込んでおこう。
例によって、計算に使用したスプレッドシートを公開しておいたので、興味ある奴はダウンロードし、値をいろいろ変えて計算してみるといいにゃ。
コメント 0