Aritalab:Lecture/Bioinformatics/Lotka - Revision history

Adm: Undo revision 304420 by Adm (talk)

2013-04-20T16:57:53Z

Undo revision 304420 by Adm (talk)

Adm: /* 競争方程式 */

2013-04-20T05:42:03Z

‎競争方程式

Adm: /* 安定性解析 */

2013-04-20T05:41:21Z

‎安定性解析

Adm: /* リアプノフ関数 */

2013-04-17T22:29:46Z

‎リアプノフ関数

Adm at 21:45, 17 April 2013

2013-04-17T21:45:09Z

Show changes

Adm: Created page with "==ロトカ・ヴォルテラ方程式== ロトカ・ヴォルテラ (Lotka-Volterra) 方程式は、被食者 (prey) ・捕食者 (predator) の関係を表す代表例で..."

2013-04-15T23:34:12Z

Created page with "==ロトカ・ヴォルテラ方程式== ロトカ・ヴォルテラ (Lotka-Volterra) 方程式は、被食者 (prey) ・捕食者 (predator) の関係を表す代表例で..."

New page

==ロトカ・ヴォルテラ方程式==

ロトカ・ヴォルテラ (Lotka-Volterra) 方程式は、被食者 (prey) ・捕食者 (predator) の関係を表す代表例です。<ref>ロトカは人口モデルのロトカと同一人物で、ヴォルテラ (Vito Volterra) はイタリアの数学者です。それぞれが独立に発見した理論なのでこう呼ばれます。</ref>
このモデルでは、被食者 x は放っておくと増加し (ax)、捕食者によって減少します (bxy) 。捕食者 y は捕食によって増加し (dxy)、放っておくと減少します (cy)。ここで a, b, c, d は全て正のパラメータです。

:<math>\begin{align}
\frac{dx}{dt} &= ax - bxy = bx (a/b - y)\\
\frac{dy}{dt} &= dxy - cy = dy (x - c/d)
\end{align}</math>

これ以降、 x , y を適当なスケーリングによって無次元化しておきます。問題の本質は変わりません。

:<math>\begin{align}
\frac{dx}{dt} &= x ( b - y)\\
\frac{dy}{dt} &= y (-a + x)
\end{align}</math>

この解は
* <math>x(t) = y(t) = 0</math>
* <math>x(t) = 0, y(t) = y(0) e^{-at}</math>
* <math>y(t) = 0, x(t) = x(0) e^{bt} </math>
これらがxy平面において x > 0, y > 0 象限の境界をなしており、(x , y) の動く範囲を中に包含しています。平衡点は原点および (a , b) です。この点を境に dx/dt, dy/dt の符号が変化します。その符号から x = a, y = b を境に4つのエリアに分割され、半時計方向に回転する軌道の存在がわかります。

===リアプノフ関数===
方程式は
:<math>(\frac{a}{x} - 1) \frac{dx}{dt} + (\frac{b}{y} - 1) \frac{dy}{dt} = 0</math>
を満たすので
:<math>\frac{d}{dt}[ a\log x - x + b \log y - y] = 0</math>
が成立します。すなわち <math>H(x,y) = a\log x - x + b\log y - y = const.</math> です。
Hを x で微分すると (a/x - 1), 二階微分は - a / x2 なので x は a で最大値をとります。同様に y は b で最大値をとります。H は (a , b) を中心に - ∞ まで減少し、その等高線は閉曲線を成します。つまり軌道は周期的です。
このように、平衡点の安定性を示すのに使う関数をリアプノフ関数と呼びます。
<ref>リアプノフ関数は全ての平衡点について必ず見つかるわけでもなく、特に決まった定義はありません。</ref>

周期の振幅と振動数は a , b に依存します。平均値を調べてみましょう。ちょっと技巧的ですが log x(t) の時間微分を積分します。
:<math>\begin{align}
\frac{1}{T}\int^T_0 \frac{d}{dt} \log x(t) dt &= \frac{1}{T}\int^T_0 \frac{dx}{dt} \frac{d}{dx} \log x(t) dt = \frac{1}{T}\int^T_0 x (b - y) \frac{1}{x} dt \\
\log x(T) - \log x(0) &= 0 = b - \frac{1}{T}\int^T_0 y(t) dt
\end{align}</math>
つまり y の時間平均は b です。同様に x の時間平均は a です。

==種内競争を持つ被・捕食者系==
ロトカ・ヴォルテラ方程式は捕食者がいないと被食者が指数的に増えてしまうので、これをロジスティック式で補正します。

:<math>\begin{align}
\frac{dx}{dt} &= x (b - y - ex)\\
\frac{dy}{dt} &= y (-a + x - fy)
\end{align}</math>

このとき <math> \frac{dx}{dt} = 0 , \frac{dy}{dt} = 0 </math> となるのはそれぞれ <math> y = -ex + b, fy = x - a </math> という直線になります。

二本の直線はy切片がそれぞれ正と負です。第一象限で交わる場合もあるし、そうでない場合もあるので、まず交わる場合を考えます。直線<math> y = -ex + b</math>の上では x 方向の傾き dx/dt が 0 になり、y 方向の傾き dy/dt が交点より左は負で減少、交点より右は正で増加します。
同様に <math> fy = x - a</math> 上では y 方向の傾きが 0, 交点より左は x 方向の傾きが正になり増加、交点より右は減少になります。これら時間微分 dx/dt, dy/dt の符号から、反時計回りに回転することがわかります。

第一象限で交わらない場合も議論は同じです。x 軸に相当する部分 (y = 0) を考えると、
<math> y = -ex + b</math> の x 切片 b/e に収束します。 y 軸に相当する部分 (x = 0) を考えると、常に減少です。すなわち、二本の直線が y 座標が負の部分で交わるとき、 (0, b/e) が収束先になります。

===リアプノフ関数===
ロトカ・ヴォルテラ方程式と同様にリアプノフ関数を考えてみましょう。
:<math>H(x,y) = u \log x - x + v \log y - y </math>
ここで二本の直線の交点が
:<math>(u , v) = \Big(\frac{a+bf}{ef+1}, \frac{b-a}{ef+1}\Big)</math>
です。x座標および分母が必ず正になる点に注意します。
交点では<math> v = -eu + b, fv = u - a </math>が成立するのでこれを使うと
:<math>\begin{align}
\frac{H(x,y)}{dt} &= \frac{\partial H}{\partial x} \frac{dx}{dt} +
\frac{\partial H}{\partial y} \frac{dy}{dt} \\
&= (\frac{u}{x} - 1) \frac{dx}{dt} + (\frac{v}{y} - 1) \frac{dy}{dt} \\
&= (\frac{u}{x} - 1) x (b - y - ex) + (\frac{v}{y} - 1)y (-a + x - fy)\\
&= (u-x) [(v + eu) - y - ex] + (v - y) [(-u + fv) + x - fy]\\
&= e(u-x)^2 + f(v-y)^2 \geq 0\\
\end{align}
</math>
つまり H(x,y) は (u , v) においてのみ 0 となり（安定）、それ以外の場合は安定点に向かって上り勾配です。さらに (u , v) は大域安定でもあります。

この定性的解釈は、如何に種内競争の項 e や f が小さくても、存在しさえすれば、それらが大域的安定点を作り出すことを示しています。つまり、ロトカ・ヴォルテラ方程式の周期解は構造的に安定ではありません。

==安定性解析==

一般的な微分方程式系を考えます。

:<math>\begin{align}
\frac{dx}{dt} &= f(x, y)\\
\frac{dy}{dt} &= g(x, y)
\end{align}</math>

定常状態を (x*, y*) とし、そこから少しずれた位置の挙動を調べます。それには (x* + u, y* + v) とした差分 u , v をテーラー展開します。ずれが微小であることから2次以上の部分は無視します。

:<math>\begin{align}
\frac{du}{dt} &= f(x^* + u, y^* + v) = f(x^*, y^*) + f_x(x^*, y^*)u + f_y(x^*, y^*)v + \cdots \\
\frac{dv}{dt} &= g(x^* + u, y^* + v) = g(x^*, y^*) + g_x(x^*, y^*)u + g_y(x^*, y^*)v + \cdots
\end{align}</math>

これを行列の形に書くと

:<math>\frac{d}{dt}\begin{bmatrix} u \\ v \end{bmatrix}
= \begin{bmatrix} f_x & f_y \\ g_x & g_y \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}
= J \begin{bmatrix} u \\ v \end{bmatrix}
</math>

行列 J をヤコビ行列 (Jacobian) といいます。線形微分方程式は J の固有値 λ1, λ2 を用いた解があります。

: λ1 ≠ λ2 のとき
:<math>\begin{bmatrix} u \\ v \end{bmatrix}
= \begin{bmatrix} c_{11} \\ c_{21} \end{bmatrix} \exp(\lambda_1 t) +
\begin{bmatrix} c_{12} \\ c_{22} \end{bmatrix} \exp(\lambda_2 t)
</math>

: λ1 ＝ λ2 ＝ λ のとき
:<math>\begin{bmatrix} u \\ v \end{bmatrix}
= \begin{bmatrix} c_{11} \\ c_{21} \end{bmatrix} \exp(\lambda t) +
\begin{bmatrix} c_{12} \\ c_{22} \end{bmatrix} t \exp(\lambda t)
</math>

固有値が実数の時、λ1 < 0, λ2 < 0 であれば、平衡点からのズレ u , v は時間とともに 0 に収束します。どちらかでも正であれば、ズレは発散します。
固有値が複素数の時、２つの固有値は複素共役になり λ1, λ2 = a ± i b です。虚数の exp は絶対値が 1 になることを利用して
:<math>\lim_{t \rightarrow \infty}| \exp(a \pm ib)t | = \lim_{t \rightarrow \infty} \exp(at) \times | \exp(\pm ibt) | = \lim_{t \rightarrow \infty} \exp(at) </math>
です。虚数の部分は影響を与えません。固有値の実部が共に負であれば、平衡点からのズレはやはり 0 に収束することがわかります。
この議論は多元連立の場合もそのまま成立することがわかるでしょう。

まとめると、
* 固有値の実部が全て負であれば、平衡点は局所安定
* 固有値の実部が一つでも正であれば、平衡点は不安定
となります。固有値が複素数の時は近傍で振動がおこります。

2 x 2 行列 <math>A = \begin{bmatrix} a & b\\ c & d\end{bmatrix}</math> の場合は固有方程式 x2 - (a + d) x + (ad - bc) = 0 の解が固有値です。
固有値実部が負であるための必要十分条件は tr(A) = a + d < 0 かつ det(A) = ad - bc > 0 です。

==競争方程式==

種内競争を持つ被・捕食者系とは係数の符号が異なる2種の競争系を考えましょう。
:<math>\begin{align}
\frac{dx}{dt} &= x (1 - ax - by)\\
\frac{dy}{dt} &= y (1 - cx - dy)
\end{align}</math>
このとき安定点は
:<math>(u , v) = \Big(\frac{d-b}{ad-bc}, \frac{a-c}{ad-bc}\Big)</math>
です。
ヤコビ行列は
:<math>J = \begin{bmatrix} -a u & -b u\\ -c v & -d v \end{bmatrix}</math>
です。ロトカ・ヴォルテラ式の場合と異なり、二本の直線はいずれも傾きが負です。直線どうしは第一象限で交わる場合もあるし、交わらない場合もあります。

交わる場合から考えます。ad − bc > 0 の場合、交点 (u , v) の分母は正になり、交点が第一象限にあることから d > b > 0, a > c > 0 です。ヤコビ行列の固有値は負になります。
このとき交点は安定点となり、dx/dt, dy/dt の符号を考えると (u , v) に落ち込むことがわかります。

次に ad − bc < 0 の場合を考えます。交点の分母が負になるので b > d > 0, c > a > 0 です。このとき、x 軸上 (y = 0) の切片 1/a および y 軸上 (x = 0) の切片 1/d に収束します。初期条件に依存してどちらかの種が絶滅することがわかります。

===リアプノフ関数===

競争方程式のリアプノフ関数を考えましょう。
:<math>H(x,y) = u \log x - x + v \log y - y </math>
二本の直線の交点は
:<math>(u , v) = \Big(\frac{d-b}{ad-bc}, \frac{a-c}{ad-bc}\Big)</math>
です。
交点では<math> 1 = au + bv, 1 = cu + dv</math>が成立するので
:<math>\begin{align}
\frac{H(x,y)}{dt} &= \frac{\partial H}{\partial x} \frac{dx}{dt} +
\frac{\partial H}{\partial y} \frac{dy}{dt} \\
&= (\frac{u}{x} - 1) \frac{dx}{dt} + (\frac{v}{y} - 1) \frac{dy}{dt} \\
&= (\frac{u}{x} - 1) x (1-ax-by) + (\frac{v}{y} - 1)y (1-cx-dy)\\
&= (u-x) [a(u-x) +b(v-y)] + (v - y) [c(u-x)+d(v-y)]\\
&= a(u-x)^2 + (b+c)(u-x)(v-y) + d(v-y)^2\\
\end{align}
</math>
二次形式の形をとり (u , v) では 0 になっています。

<references/>

@@ Line 138: / Line 138: @@
 ==競争方程式==
 安定性解析をベースに、被・捕食者系を捉え直してみます。
-ここで係数の符号が異なる線形の競争系を考えましょう。
+ここで係数の符号を変化させた競争系を考えましょう。
 :<math>\begin{align}
-\frac{dx}{dt} &= 1 - ax - by\\
+\frac{dx}{dt} &= x (1 - ax - by)\\
-\frac{dy}{dt} &= 1 - cx - dy
+\frac{dy}{dt} &= y (1 - cx - dy)
 \end{align}</math>
 種 x と y の相互作用の関係は対称形で、被・捕食者系とは y の単独の増殖率および c の符号が違っています。安定点は
 :<math>(u , v) = \Big(\frac{d-b}{ad-bc}, \frac{a-c}{ad-bc}\Big)</math>
-ヤコビ行列は
+(u , v) におけるヤコビ行列は
 :<math>J = \begin{bmatrix} -a  & -b \\ -c  & -d  \end{bmatrix}</math>
-です。被・捕食者系と異なり、二本の直線はいずれも傾きが負です。直線どうしは第一象限で交わる場合もあるし、交わらない場合もあります。
+です<ref>微分方程式自体は非線形なのですが x の成長率 x dx/dt および y の成長率 y dy/dt を考えると線形微分方程式になっています。解析結果は線形微分方程式のものと変わりません。</ref>。被・捕食者系と異なり、二本の直線はいずれも傾きが負です。直線どうしは第一象限で交わる場合もあるし、交わらない場合もあります。
 交わる場合から考えます。ad &minus; bc > 0 の場合、交点 (u , v) の分母は正になり、交点が第一象限にあることから d > b > 0, a > c > 0 です。さらに ad - bc > 0, (-a + -d) < 0 を満たすので、ヤコビ行列の固有値は負になります。つまり交点は大域安定な不動点です。

@@ Line 138: / Line 138: @@
 ==競争方程式==
 安定性解析をベースに、被・捕食者系を捉え直してみます。
-ここで係数の符号が異なる競争系を考えましょう。
+ここで係数の符号が異なる線形の競争系を考えましょう。
 :<math>\begin{align}
-\frac{dx}{dt} &= x (1 - ax - by)\\
+\frac{dx}{dt} &= 1 - ax - by\\
-\frac{dy}{dt} &= y (1 - cx - dy)
+\frac{dy}{dt} &= 1 - cx - dy
 \end{align}</math>
 種 x と y の相互作用の関係は対称形で、被・捕食者系とは y の単独の増殖率および c の符号が違っています。安定点は

@@ Line 87: / Line 87: @@
 ==安定性解析==
-ここで話を一般化します。一般的な微分方程式系を考えましょう。
+ここで話を一般化するため、一般的な線形微分方程式系を考えましょう。
 :<math>\begin{align}

@@ Line 26: / Line 26: @@
 :<math>(\frac{a}{x} - 1) \frac{dx}{dt} + (\frac{b}{y} - 1) \frac{dy}{dt} = 0</math>
 が成立します。つまり微分してこの値になる式は x , y に依存しません。
-:<math>\frac{dH(x,y)}{dt} = \frac{d}{dt}[ a\log x - x + b \log y - y] = 0</math>
+:<math>\begin{align}\frac{dH(x,y)}{dt} &= \frac{d}{dt}[ a\log x - x + b \log y - y] = 0 \\
-:<math>H(x,y) = a\log x - x + b\log y - y = const.
+H(x,y) &= a\log x - x + b\log y - y = const.\end{align}
 Hを x で微分すると (a/x - 1), 二階微分は - a / x<sup>2</sup> なので x は a で最大値をとります。同様に y は b で最大値をとります。つまり H は (a , b) を中心に - &infin; まで減少し、その等高線は閉曲線を成します。つまり軌道は周期的です。