2.3 解軌道

(3a), (3b) は非線形方程式ということもあって、 解を解析的に求めることは(おそらく)出来ないが、 解軌道の方程式を求めることは出来る。 この辺は有名な Lotka-Volterraロトカ・ヴォルテラ の方程式 ¹や、Kepler運動²と事情が似ている。

(3a), (3b) より

$$\frac{\D I}{\D S} =\frac{\frac{\D I}{\D t}}{\;\frac{\D S}{\D t}\;} =\frac{\beta SI-\gamma I}{-\beta SI}=-1+\frac{\rho}{S}.$$

ただし

$$\rho:=\frac{\gamma}{\beta}$$

とおいた。

$$\D I=\left(-1+\frac{\rho}{S}\right)\D S$$

を積分して、と書いてある本が多いが、 $\frac{\D I}{\D t}=\left(-1+\frac{\rho}{S}\right)\frac{\D S}{\D t}$ から、という方が安心する人が多いかな、

$$\int_0^t \frac{\D I}{\D t}\D t =\int_0^t \left(-1+\frac{\rho}{S}\right)\frac{\D S}{\D t}\D t.$$

これから

$$I(t)-I(0)=-S(t)+S(0)+\rho\left(\log S(t)-\log S(0)\right).$$

ゆえに

$$I(t)=I(0)+S(t)-S(0)+\rho\log\frac{S(t)}{S(0)}.$$

$S(0)$, $I(0)$, $S(t)$, $I(t)$ をそれぞれ $S_0$, $I_0$, $S$, $I$ と書くと

$$I=I_0+S_0-S+\rho\log\frac{S}{S_0}.$$ (4)

これは $SI$ 平面における解軌道の方程式である。

$S_0$, $I_0$ を与えて、この曲線を描くのはコンピューターを使えば簡単である。

$I$ は $S$ の関数とみなすとき、

$$\frac{\D I}{\D S}=-1+\frac{\rho}{S}=\frac{\rho-S}{S}$$

であるから、$0<S<\rho$ で増加関数、$S>\rho$ で減少関数であり、 $S=\rho$ で最大値を取る。 また

$$\frac{\D^2 I}{\D S^2}=-\frac{\rho}{S^2}<0$$

であるから、グラフは上に凸である。

$S\to+0$ や $S\to\infty$のとき $I\to-\infty$.

(力学系は、$S$ 軸上の点 $(S,0)$ を平衡点としていて、 任意の解について $\dsp \lim_{t\to+\infty}I(t)=0$ が成り立つ、それと矛盾はしないが、 混乱しないように注意が必要である。)

図 1: (22)の軌道 $\rho =1,000$ (佐藤[8] p. 176 図4)

初期時点 $(S(0),I(0))$ を出発する点 $(S(t),I(t))$ は、 時間がたつにつれて、この軌道の上を右から左に動いてゆく。そのわけは

$$\frac{\D S}{\D t}=-\beta S I<0.$$

$t\to+\infty$ のときに $(S(t),I(t))$ が $S$ 軸上の点に収束することを示すのは、 解析学の練習問題である。

(ゼミ中に細かいところまで調べたような記憶があるが、 記録が残っていないのは残念なことだ…以下は記憶を手繰って)

例えば、 $(S,I)$ を通る曲線と $S$ 軸との交点 $(S_\ast,0)$ を求めるには、 $S_\ast$ についての方程式

$$I=0+S_\ast-S+\rho\log\frac{S}{S_\ast}$$

を解けば良い。これは次のようなMathematica コードで実行できそうだ。

findS[S_, I_, rho_] := 
 Module[{S0}, 
  S0 /. FindRoot[I == 0 + S0 - S + rho Log[S/S0], {S0, rho/2}]]

S0=findS[2000, 10, 1000]

これで $399.626$ という値が得られる (反復法の初期値 $\rho/2$ が適切かは分からないが… 2階導関数が定符号なので、危なく解ける保証が出来そうだ)。 $S_0$ が求まれば、グラフを描くには Plot[S0 - S + rho Log[S/S0], {S, S0, 2000}] のようにすれば良い。

このようにして、 佐藤[8] p. 176の図4を再現してみた (図2. やってみて気付いたが、 $(1000,10)$ という座標は $(1100,10)$ の誤植であろう)。

図 2: 佐藤[8]図4の再現