第5回 2010年5月14日

東京女子大学のトップページ   情報処理センターのページ   東京女子大学の Gmail のページ   東京女子大学の図書館のページ   浅川のホームページ   授業のホームページへ戻る
ロトカ・ヴォルテラ方程式の続き

ヴォルテラの原理

捕食者と被食者の個体密度は周期的に振動し,その振幅と振動数は初期値に依存して決まります。 しかし,個体数の時間平均は一定であり休止点 LVinline007 に等しくなります。すなわち

LVeq003

ここで T は解の周期とします。

LVeq004

ですから,両辺を積分して

LVeq005

すなわち

LVeq006

を得ます。x(T)=x(0) より

LVeq007

となります。x の時間平均についても同様に LVinline008 が成り立ちます。

以上により,ダンコナの疑問,なぜ戦争はサメを好むのか, サメなどの捕食者の個体数が, 戦争前よりもかなり高くなっていたことに対する答えの準備が整いました。

漁師による漁業活動は,被食者の増加率を減少させます。 a の代わりに a-k になるわけです。 同時に捕食者の減少率を増加させる c の代わりに少し大きな値 c+m になります。 しかし相互作用を表す定数 b, d は不変です。 従って,漁師による漁業活動が中断されている場合と比べて, 捕食者の個体数の時間平均は (a-k)/b で少し小さくなり, 被食者の平均は (c+m)/d で少し大きくなります。 戦争により漁業活動が中断されている間,この逆のことが起こり,捕食者は増加し, 被食者は減少したというわけです。

リヤプノフ関数

ロトカ・ヴォルテラ方程式が安定であるための1つの十分条件は リアプノフ関数 Lyapunov function, V(x,y) が存在することであるとされます。 逆に言えば,リアプノフ関数を見つけることができれば,その力学系は安定です。 前出のロトカ・ヴォルテラ方程式

LVeq008

のリアプノフ関数は,以下のとおりです。

LVeq009

この場合のリアプノフ関数 V(x,y) は時間とともに変化せず, 状態点は V(x,y) の等高線をたどりながら動きます。 すなわち一周すればかならず元に戻ります。 このようにシステムの状態の関数であって, 時間変化にともなって増えも減りもしないものを保存量といいます。 平衡点のまわりからずれていたとするとシステムの挙動は, 別の軌道に移り,平衡点には戻らりません。 しかし遠くには慣れてしまうわけでもありません。 このような平衡点は 中立安定 といわれます。

リアプノフ関数を見つけることは必ずしも簡単ではありません。 しかし,見つかればシステムの安定性解析に役立ちます。

Lyapunov
ロトカ・ヴォルテラ方程式のリヤプノフ関数 a=1, b=0.002, c=1, d=0.02

種内競争を持つロトカ・ヴォルテラ方程式

もし,捕食者がいなければロトカ・ヴォルテラ方程式(2) の被食者は LVinline009 となり,その解は指数関数 LVinline010 となって個体数が爆発的に増加するという話は以前しました。 これは,被食者の個体数の増え方にロジスティック方程式 LVinline011 を仮定することで修正できます。 捕食者の種内競争も同じくロジスティック方程式に従うとすると, (2)の代わりに

LVeq017

を考えます。 この場合周期解になるとは限りません。 パラメータのとり方によってはある種が絶滅する場合もありえます。 一例を下図に示します。

LotkaVolterra2
種内競争をもつロトカ・ヴォルテラ方程式

実習

それでは実習です。

java LotkaVolterra2

してください。 上の図を描くシミュレータが起動します。 a, b, c, d, e, f のパラメータ, 及び初期値 x, y を変えて遊んでみてください。

アイソクライン

このときの の平衡点を解析するために少し変形してみましょう。 x=0, y=0 すなわち両個体がいないときには変化がないのはすぐ分かります。 それ以外は, LVinline012 傾き 0 の アイソクライン isocline (iso とは「同じ」という意味であり,cline とは「傾き」という意味)

LVeq018

と,同様に LVinline014 のアイソクライン

LVeq021

を考えます。 この二直線の交点が平衡点となる。平衡点の座標は LVinline013 となります。 このとき x と y との増減は以下の表のようになります。

アイソクラインによる x, y の増減
isocline11 isocline11
isocline11 x, y とも減少x は増加,yは減少
isocline11 x は減少,y は増加x,yともに増加

すなわち,アイソクライン線の上下で x と y との増減を考えることができるわけです。 この方法のことをアイソクライン法と呼んでいます。

ここで,パラメータをいじっていろいろ遊んでみましょう。 たとえば,a=8, b=1, c=-6, d=-1, e=1, f=1 を考えます。 この場合, LVinline015 だけが安定な平衡点になります。 パラメータにマイナスを使っているので, もはや捕食者,被食者と呼ぶのは正しくないかも知れません。 そこで x 種 と y 種と表現することにすると,この場合 y 種は滅びる。 最初にどんなに y 種の個体数が多くても LVinline016 は不安定平衡点であるので y 種は生き残ることができないのです。

a=8, b=1, c=-6, d=-0.3, e=1, f=1 を考えます。 この場合,両アイソクライン線の交点が安定な平衡点になって, どのような初期値から出発しても (3,5) に収束します。すなわち x 種と y 種との共存が実現するわけです。

a=8, b=1, c=-10, d=-1, e=0.5, f=1 を考えてみましょう。 この場合,両アイソクライン線の交点は不安定平衡点になる。初期値によって y 種だけが生き残る LVinline017 か,x 種だけが生き残る LVinline018 に収束する双安定になります。

このようにパラメータを変化させると,2 種の共存か, 競争によって 1 種だけが勝ち残ることが生じます。 たとえば,2 種のトカゲが多数の島に棲んでいたとしましょう。 気候条件や栄養状態, 島の大きさなどがほぼ同じであっても, 両種が生息している島がないとすれば, 競争関係のために共存できなかったと考えられるわけです。 そのため,環境の変動や個体数などのランダム性のため島の個体群が滅びることが起こりというわけです。

もう一つの安定なロトカ・ヴォルテラ方程式,リミットサイクル

ロトカ・ヴォルテラ方程式をもう少し現実的にしてみましょう。

LVeq024

というモデルを考えます。 オリジナルのロトカ・ヴォルテラ方程式からの変更点は, 以下のとおりです。

  1. 捕食者のいないときに被食者は無限に増殖することはせず, ロジスティック方程式に従って環境収容力 K に収束する。 これは y=0 とおけば確かめられます。
  2. 捕食者の摂食速度を LVinline019 としたこと。 これは,餌が多いときには捕食速度は餌の量とともに増加します, あまりに多いと食べきれなくなるため飽和し,上限値 a/h を持ちます。

K が有限で h=0 であれば被食者の個体群を安定させる働きによってロトカ・ヴォルテラ方程式 の振動は止まります。

LotkaVolterraAttractor
アトラクター

一方,捕食速度の飽和傾向だけを考慮すると K = ∞, h>0 システムは大域的に不安定になり,振動の振幅が時間とともに大きくなる。

この両方の効果が加わると 極限周期道(リミットサイクル) limit cycle が現れます。

LotkaVolterraLimitCycle
リミットサイクル r=0.18, K=30.0, a=0.02, h=0.1, b=0.0504, c=0.25

実習

java LotkaVolterra3

してみてください。アトラクターやリミットサイクルが現れるパラメータを探してください。

殺虫剤に害虫の駆除は効果があるか

害虫 x と天敵 y を考え,

LVeq025

に殺虫剤の効果 -mx が加わったモデルを考えてみましょう。

LVeq027

殺虫剤を散布すると害虫 x は減るのでしょうか。 もちろん殺虫剤を散布すれば害虫の数は減ります。 ところが(28)の平衡状態を調べると m があっても平衡状態は変化しないことが分かるのです。 つまり害虫の増殖率を下げる効果が天敵 y のレベルを下げ, その結果,害虫の増殖率を上げるように働く効果があるのです。 これらの効果が打ち消しあってしまいます。

実際には,殺虫剤の効果 m は,天敵 y の生存率も下げるでしょう。 すると(28)の右辺にも -my という項を付け加えるモデルが考えられます。

LVeq029

この場合,殺虫剤の散布はかえって害虫 x の個体数を増加させてしまいます。 このように, 殺虫剤の散布がかえって害虫の個体数を増加させてしまうことはしばしば見られることなのです。 害虫が突然変異によって殺虫剤への耐性を進化させたと考えるよりも, 案外こういうところに原因があるのかも知れません。

ロトカ・ヴォルテラモデルの拡張

捕食者と被食者が多様な生態系では,ロトカ・ヴォルテラ方程式の多変数拡張が 用いられます。

LVeq031

種間の相互作用 aij は,i 種が j 種の捕食者であるなら a >0 であり,aji<0 となる。

たとえば,光合成によって有機物を作る植物は生産者と呼ばれ,それを食べる 一時消費者(草食動物),さらにそれを食べる2次消費者(肉食動物)というように 多様な生態系が形作られます。 この(31)の安定性解析を行うことで,現実の生態系の実体に近づくものと思われます。

参考文献

線形微分方程式系

線形微分方程式系

生理現象,生命現象を扱った微分方程式の多くは非線形です。 したがって,解析的な解が存在しないか,あるいは, 存在しても求めることが難しいことが多いです。 この非線形性がモデルの本質ではあります。 しかし,部分的には線形の微分方程式に還元できる場合もあるのです。 ここでは,そのような線形微分方程式の性質について見てみましょう。

非線形微分方程式の線形化

関数 f, g を未知とする非線形微分方程式系,

NEQ001

があるとします。この式の平衡点を NEQ002 とします。 すなわち,

NEQ003

を満たすものとします。このとき f(x,y) を平衡点の近傍で線形化することを考えます。 f(x,y) を平衡点の近傍でテーラー展開し,

NEQ005

同様に g(x,y) についても,

NEQ006

として,高次の項を無視します。 以上より,平衡点におけるヤコビ行列

NEQ006

が定義できるので,

NEQ006

を得ます。以下では,この連立微分方程式系を解くことを考えてみましょう。

定式化

ロトカ・ヴォルテラ方程式の捕食者と被捕食者との関係を, 連立微分方程式(微分方程式系ともいう)としてとらえます。 x1, x2 という 2 種類の動物を考え, この動物間で共存,競合などの関係が生じると考えます。 するとこれら 2 種の個体数の時間変動は,

EQ01
EQ02

と表わせます。 LEQ01, LEQ02, LEQ03, と行列表現すれば,上式は

EQ04

と書くことができます。

実は,この考え方は重要なのです。今,線形とは限らない微分方程式系

EQ05

を考えます。 f を平衡点, つまりすべての t に対して LEQ04 をみたす点 mbx の近くでテイラー展開し,一次の項までとれば, (4)になるわけです。

さて,(4)の右辺を見てみると, この行例 mbA は, mbR2 を線形変換

EQ06

によって,写像することだと考えることができます。

この様子をみるために,実際に mbx mbx とを実際に, mbR2 上に描いてみましょう。 それには,ベクトル LEQ05 を原点 O から座標 x1, x2 を持つ点 P に向かう矢線と考え, 終点 P からベクトル mbx に対する矢線を惹くと分かりやすいです。 たとえば,

EQ07

としてみましょう。 LEQ06 に対して点 LEQ07 に注目します。次に, LEQ08 に対応する矢線 OQ を描きます。 最後に,OQ の始点が,点 P にくるように OQ を平衡移動します。

FigLV1

このようにしてできた mbR2 上の軌道,曲線のことを LEQ09 の相曲線といったりします。

このようにして, mbR2 上の各点に対してベクトルを描き入れるとベクトル場が描けます。

FigLV2

ところで,原点 LEQ10 は常に微分方程式(2)の解となります。 したがって,原点 LEQ11 も(2)の特別な軌跡と考えることができます。 このような点は,微分方程式の平衡点,または特異点と呼ばれます。

一般に,軌道が描かれる平面 mbR2 は微分方程式(2)の相平面と呼ばれます。 相平面は,微分方程式の解曲線で隙間なく埋め尽くされていますが, その曲線は交わることはありません。

行列の指数関数

微分方程式 dx/dx = ax の解は,初期値を x0 とすれば, x(t) = x0eat と表わせました。 今度は微分方程式系 LEQ12 の解で LEQ13 を満たす解を求めることを考えます。 この2つは同じパターンですから, この解を

EQ08

とおいてみます。指数関数 LEQ14 は,テーラー展開を使って,

EQ09

でしたね。

ちょっと脱線,テーラー展開は
EQ10
でしたねー。

とにかく, このように指数関数はべき級数に展開できるのですから, (9)のまねをして,

EQ12

とおいてみます。ここで mbI は単位行列です。

右辺の級数は次ようのな行列です。いま

EQ13

LEQ15 の要素を

EQ14

LEQ16 と書いておきます。

(12)の右辺第 k 項は 2×2 行列

EQ15

ですから,この式の LEQ17 要素は t のベキ級数

EQ16

ただし

EQ17

です。 この級数は任意の実数 t に対して収束することが示されます。 そこでこの LEQ18 要素を行列 LEQ19 と定義することにします。 これを行列 LEQ20 の指数関数と呼ぶことにします。 これは次のような性質を持ちます。

EQ18

練習問題

EQ22

の解を求めてみましょう。

LEQ21 ですから,

EQ23
LEQ22

などとなるので,

EQ24

これにより解 LEQ23 は,

EQ28

となります。

ちなみに,

EQ29
EQ30

とおくと,

EQ31

となります。つまり,

EQ32

と表せるわけです。行列 LEQ24 は,

EQ33

となるので,オイラーの公式

EQ34

の行列バージョンということができます。

固有値と固有ベクトル

線形微分方程式系

EQ35

を実際に解いてみましょう。たとえば,

EQ36

は,任意の実数 λ によって

EQ37

が成り立ちます。一般に線形変換 mbx に対して(37)を みたすベクトル LEQ25 LEQ25 のなかに存在するとき, 実数 λ を線形変換 mbA の固有値,そのときの mbx を固有値 λ に対する mbA の固有ベクトルと言いましたよねー。

一般に,行列 LEQ27 が 2 次の正方行列であれば,

EQ38

の場合

EQ39

が必要十分条件です。この場合, (a - λ)(d - λ) - bc = 0 すなわち,

EQ40

の根を固有値 λ と言いました。(40)を線形変換 mbA の特性多項式と言ったのでしたね。

行列 mbA の固有値は次のようにして求めたのでした。

EQ41
故に, (λ - 1)(λ - 3) = 0 ということになります。 λ1=1 に対応する固有ベクトルを (p11, p21)T とすると,

EQ42

だから,固有ベクトルは方向さえ表せれば良いので,たとえば LEQ28 となります。 同様にして,λ2=3 に対応する固有ベクトルは, LEQ29 となります。

固有値は,線形微分方程式系の安定性と関係があります。 固有値は,その固有値に対応する固有ベクトルの方向に状態点の動きの速さを示すことになるからです。 上の例では,固有値が正なので,軌道は原点から離れていく不安定ノードになります。

FigLV3

実習

次の各式を解いて,固有値を求め,対応する固有ベクトルを求めてください。 同時に,グラフを描いて固有値とグラフとの関係を考察してください。

問題 1

EQ43

鞍点: 特性多項式は,LEQ30 ですから,固有値 LEQ31, LEQ31, を得ます。

問題 2

安定結節点

EQ44

問題 3

EQ45

問題 4

EQ46

問題 5

EQ47

問題 6

EQ48

問題 7

EQ49

問題 8

EQ50

問題 9

EQ51

ただし,λ < 0

問題 10

EQ52

ただし,λ > 0

問題 11

EQ53

問題 12

EQ54

問題 13

EQ55

問題 14

EQ56

問題 15

EQ57

問題 16

EQ58

問題 17

EQ59

問題 18

EQ60

以上をまとめると, 固有値が 2 つとも正のとき,2 つとも負のとき,正と負のとき, 複素数で実部が正のとき,実部が負のとき,実部が 0 のとき, などに分けて考えることができます。