第4回 2010年5月7日

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

先週終わらなかったところからー

安定不動点,不安定不動点,カオス

さて,先にも記したとおり,個体成長が,ある法則 f に従うとすると, 現在の状態を x として,次の世代の状態は y は, y = f(x) と表せました。 このとき,f は,個体数が少ない時には,成長し(増加し),多い時には減少する と考えれば,y = f(x) という関数は,直線によって表せるものではなく, コブをもった関数になるでしょう。

一つだけコブを持った関数としては,一番簡単なのは, おそらく中学生で学習する放物線,すなわち x に関する 2 次関数 y = a x2 + b x + c でしょう。 そこで,以下のような単純な関数に関して, 次世代がどのように変動するかを見ていくことにしましょう。

0 と 4 との間の R に対して写像

eq22

は区間 [0, 1] 上で力学系を与えます。 点 0 は明らかに不動点です。 R < 1 ならば,全ての x ∈ (0, 1] に対して xt+1 < xn です。 従って,x の起動は 0 に向かって単調に減少します。 今後は,R > 1 の場合だけを考えます。 F のグラフは点 0 と 1 で x 軸と交わり, 点 1/2 で最大値 R/4 をとる放物線です。 放物線は一辺の長さ 1 の正方形の内部で, 直線 y = x と唯一の点 P で交わります。 P の横座標値は p は pn+1 = pn を満たすので, もう一つの浮動点 p = (R−1)/R が得られます。

平衡状態は xt+1 = xt とおいて計算されるので, 図上では y = x との交点として求められます。 xt に対する次世代 xt+1 の値は曲線を使って求められますが,次に縦軸にある値を横軸に戻すには y = x を用いると直感的でわかりやすいです。このようにして, xt+2, xt+3 と順次グラフを使って, 個体数を求めることができます(下図参照)。

attractor  attractor2

p の安定性に関しては,平均値の定理より,x と p との間に適当な c を置くと

eq23

が成り立ちます。x (従って c ) が p に十分近ければ, inline_eq1 のとき inline_eq2 が成り立ちます。従って,

eq24

となります。 すなわち,x' は x より浮動点 p に近いのです。 このことは,もし x が p のある適当な近傍から選ばれれば,x の軌道 (22 から反復的に得られる数列)は p の近傍に留まり, さらに p に収束することを意味します。 このような意味で,不動点 p は 漸近安定 であると言います。 逆に, inline_eq3 であるならば,

eq25

となるので,x の軌道は不動点から遠ざかります。 このような場合 p は 不安定 であるといいます。

今,F(x) = R x (1-x) なので, inline_eq4 となり,p は

カオスでは, ほんの少しだけ違った初期値ではじめると十分時間が経ったときに, その値は予想できなくなります。 3.5699 はカオスの境界を表しており, ファイゲンバウム点 と呼ばれたりします。

attractor3
R = 4 の時カオスになる

大事な点は,計算可能生は予測可能生を意味しないということです。 決定論的な運動がランダムな動きと区別できない ことが重要です。

R=2.9 のときのグラフを次に示します。

Chaos2

R=4 のときカオスとなります。これを下の図に示します。

Chaos1

実習

端末エミュレータ(Windows ではコマンドプロンプト)から

java LogisticMapping

とタイプしてください。(22) すなわち xt+1 = R xt ( 1 - xt) の計算が始まります。R と初期値を変化させて遊んでみてください。 特に,R がカオス領域のときに,わずかな初期値の変化が, ドラスティックな変化をもたらす点に注目してください。

参考文献

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

ロトカ・ヴォルテラ方程式の発見

第一次世界大戦中, アドリア海はイタリア海軍とオーストリア・ハンガリー帝国との戦いのために, その海域での大規模な漁業は全て中止させられてしまったそうです。


大きな地図で見る
図中Aは,アドリア海の真珠と呼ばれるドブロクニク

その結果,予期せぬ事態が起こりました。 終戦後, 数年してからイタリア人の生態学者ダンコナ D'Ancona は, 漁業市場で漁獲量の統計を取った結果, 戦争後サメなどの肉食魚の割合が, 戦争前よりもかなり高くなっていたことに気づいたのです。 ダンコナは困惑しました。なぜ戦争はサメを好むのだろうか?と。

ダンコナは,ローマの数理物理学の教授であり, その当時第一線の数学者として広く知られていたヴォルテラ Senator Vito Volterra に問い合わせました。 その結果,今日ロトカ・ヴォルテラ方程式として知られる微分方程式を使って, アドリア海における生態系の説明モデルが産まれました。

アドリア海には多数の魚が棲息していますが,単純化すると 2 種類に分けられます。 サメなどの捕食者とその餌となる魚です。 ヴォルテラはこれら 2 種類の魚が相互に依存していることを見抜きました。 すべての捕食者は,その獲物(被食者)を食べることで,獲物の増加率を減少させます。 もし,捕食者の数があまりに増え過ぎると,その獲物である被食者の個体数は減少することになります。 逆に獲物(被食者)の個体数が大きくなると,捕食者は繁栄します。 また,もし被食者の数が少なくなり過ぎると,捕食者は子孫をつくる前に餓死してしまい, その結果捕食者の個体数は減少することになります。

このようにして捕食者であるサメの個体数は増えたり減ったりします。 しかし絶滅することはありません。 なぜなら,サメの個体数が減ると必ず餌(被食者)である魚は増殖します。 そうなるとサメ(捕食者)は餌である魚(被食者)を捕まえやすくなるからです。 すなわちサメの個体数が増えるようになります。 また,もし魚(被食者)の個体数が少なければ, サメの個体数も食料不足のために減少して, 結果的に餌である魚(被食者)にとって楽な環境になるというわけです。

つまり,サメと餌である魚の個体数は増えては減り,減っては増えるという具合に循環するわけです。

さてここで,この捕食者と被食者のゲームに漁師という登場人物を加えることにしましょう。 漁師は,捕食者であるサメも被食者である魚もどちらも捕ることができるので, 両者の総数はともに減少します。 しかし,漁の後には,捕食者ー被食者の間のバランスが崩れ, 予期せぬ結果が生じます。 漁によって被食者である魚の数が減少すると,捕食者であるサメの個体数もまた減少します。 なぜなら,サメの餌を漁師に捕られて食料不足になってしまうからである。 つまり結果として,捕食者の個体数が減少するため, 獲物である魚の数がかなり驚異的に増えてしまいます。 もし戦争によって漁が中止されると,その逆のことが起こります。 つまり,捕食者の個体数は増え,獲物である被食者の数は少なくなるわけです。

これが今日「ヴォルテラの法則」として知られるようになった現象です。 すなわち,捕食者ー被食者の両者が互いに相手の個体数を決めるという前提があるなら, 両者の成長率の低下,たとえば漁など,は被食者の増殖と捕食者の減少を招くというものです。

ヴォルテラのモデルは, アメリカ人の化学者ロトカ Alfred J. Lotka によってすでに研究されていたことが分かりました。 それで,今日では,ロトカ・ヴォルテラの方程式と呼ぶようになったのです。

捕食者-被食者方程式

ヴォルテラは,捕食者が存在しない場合,被食者の個体群の成長がある定数 a で与えられ, さらに成長率は捕食者密度 y の関数として線形に減少すると仮定しました。 これから dx/dt = a - by ただし(a,b>0) が得られる。 被食者が存在しない場合,捕食者は死亡せざるをえないので捕食者の成長率は負です。 しかし被食者密度 x が正であるならば, 成長率は回復するので dy/dt = -c + dx (c,d>0) となります。 両者をまとめると次式を得ます。

LVeq001

この式がロトカ・ヴォルテラ方程式と呼ばれるものです。 この式の解は,次の 3 つはすぐに分かるでしょう。

  1. LVinline001
  2. LVinline001
  3. LVinline001

すなわち,捕食者の個体密度 y も被食者の個体密度 x も, ある時刻 t において 0 であるならば,常に 0 です。 被食者が不在の場合 x(t)=0,捕食者は絶滅に向かいます (t → +∞ で y(t) は 0 に収束する。 捕食者がいなければ y(t)=0,被食者の個体数は爆発的に増えます (x(t) → +∞)

捕食者を横軸,被食者を縦軸にして解曲線を描いたものが 下図です。このように変量 x と変量 y とを縦軸,横軸にとって その変化の様子を描いたものを相図 Phase Graph と言います。

LotkaVoltraPhase
ロトカ・ヴォルテラ方程式の相図

捕食者の個体数が増えると被食者個体数が減り,その結果捕食者個体数が減り始めます。 その結果,被食者個体数が増加する,という円環状の関係になっています。 初期値が違えばこの円環の大きさが変化します。しかし,その形は変化しません。 円環はパラメータを変化させることで形状が変わります。

実習

いつものように実習してみよー。

cd
cp -rp ~asakawa/20100423 .
cd 20100423

とすれば準備は完了です。ロトカ・ヴォルテラ方程式のシミュレータを動かすには

java LotkaVolterra

とタイプすれば良いです。

シミュレーションを実行してみると分かりますが, ロトカ・ヴォルテラ方程式には平衡点が存在します。 すなわち (x,y)=(0,0) という捕食者と被食者ともいない場合, 加えて, 両者の個体数が均衡して変化しない点が存在します。 このような平衡点は LVinline004 LVinline005 LVinline006 を満足しなければなりません。従って

LVeq002

です。しかし,この平衡点から少しでも離れれば,振動を永遠に繰り返すことになります。

ロトカ・ヴォルテラ方程式の時間発展

同じロトカ・ヴォルテラ方程式の時間発展を描いたものが 下図です。

LotkaVolterraTime

実習としては,

java LotkaVolterraTime

としてください。

ヴォルテラの原理

捕食者と被食者の個体密度は周期的に振動し,その振幅と振動数は初期値に依存して決まります。 しかし,個体数の時間平均は一定であり休止点 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

してください。 上の図を描くシミュレータが起動します。

アイソクライン

このときの の平衡点を解析するために少し変形してみましょう。 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)の安定性解析を行うことで,現実の生態系の実体に近づくものと思われます。

参考文献