お前らに質問（３月１６日） [お前らに質問]
お前らに質問（３月１６日）
たとえば、
という連立常微分方程式（ロトカ・ヴォルテラの方程式という）がある。ここで、a、b、c、dは正の定数だケロ。
そして、ｔとx,yに関する次のデータがある。
この表のデータをもとに、（１）式に出てくる実定数a、b、c、dのおおよその値を求めたいんだけれど、ネムネコは、このデータからどうやってa、b、c、dの推測したらいいのか、この方法をまったく思いつかずに困っている。
そこで、お前らに質問だにゃ。
【質問】
上の表のデータから、どうやって、（１）の連立常微分方程式に登場する定数a、b、c、dの値を推測したらいいでしょうか。
その方法とその方法で推測したa、b、c、dの値をネムネコに教えよ！！
「ｔとx、yに関するデータが７組あるので、この関係をtの６次の多項式で近似して…」とか考えたんだけれど、６次だと近似曲線がウネリすぎて、これは使い物にならない。
ラグランジュ補間などの高次の多項式で実験データを近似すると、このようなウネウネ現象が発生するだケロ。
このことは６次式で近似する前にわかっていたことだけれど、話のネタとして提示したまでだにゃ。
ここ↓を見ると、
http://www.aoni.waseda.jp/yuuka/Sim/Lotka.html
a=1、b=0.2、c=1、d=0.031
と出ていて、この値をもとにシミュレーションしてるけど、このシミュレーションの結果は、どう贔屓目に見ても、データと一致していないにゃ。
引用：http://www.aoni.waseda.jp/yuuka/Sim/Lotka.html
この記事を書いたのは、どうも、早稲田の先生みたいなんだけれど、この先生は、いったい、どうやって、a、b、c、dの推測値を求めたんだろう。
わからないにゃ。困ったケロよ。
この先生の記事には出ていないけれど、参考までに、横軸にウサギ、縦軸にヤマネコの頭数ををとり、位相図（(・・?）を書くと次のようになるにゃ。
反時計回りに、閉曲線上を永遠に回り続けるんだケロよ。
こういう実験データの処理は、何かとその経験を有するddt³さんが得意そうだから、教えてくれないかな(^^ゞ。
できたら、魚と鮫の場合のa、b、c、ｄの値も(^^)。
どうでもいい与太話だけれど、何でも、この微分方程式は鮫と深い関係があるらしいんだケロ。
戦争が起きると、何故か、鮫が増え、鮫に食べられる他の魚が減る。このことに疑問をもった生物学者のダンコナさんが数学者であるヴォルテラに問い合わせ、ロトカ・ヴォルテラ方程式が生み出されたという話。
と少し脱線したところで、話をもとに戻そう。
普通、この手の問題は、（１）式の係数a、b、c、dの値とx(t)、y(t)の初期値が与えられていて、あるいは、自分で勝手にいくつか値を設定して、この微分方程式を（数値的に）解けという問題だけれど、これはデータから係数を求めよという一種の逆問題。
逆問題は難しいんだケロよ。
しかも、シミュレーション結果を見ると、どのような値を設定しても、データをうまく再現出来ないことはわかりきっているし、ネムネコは、この時点で、もう戦意を喪失してしまう。
こんなもの、どう頑張ったって、データと合われることなんてできないニャ。
ホント、困ったもんだケロよ。
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の場合について、オイラー法で数値的に解いた結果を示すにゃ。
これだと、渦巻きを描きながら、パンデミックの定常解に近づくことがよくわかると思うにゃ。
記事の内容とは関係ないけれど、この動画も埋め込んでおこう。
例によって、計算に使用したスプレッドシートを公開しておいたので、興味ある奴はダウンロードし、値をいろいろ変えて計算してみるといいにゃ。
