ロジスティック方程式

ロジスティック方程式とは，数理生態学の最も基本的な常微分方程式である．方程式は，以下で与え られる．

dN

dt (t) =αN(1−λN(t))  (α, λ≥0)

ロジスティック方程式の意味は次のようにして考える．すなわち，N(t)をある生物種の個体数とした とき，その変化率（増殖率）は，dN

dt は，現在の個体数N(t)に応じて決まるので，M N(t)とかける．も し増殖率が現在の個体数に比例して決まるのであれば，Mは定数で，指数関数の方程式が得られる．し かし，通常，生物種が増えすぎれば，環境の劣化（餌となる植物の減少など）によって，増殖率は抑制

される．そのような場合の最も簡単なモデルとしては，M =a−bN(t)の形で与えられる場合が考えら れる．これを整理してパラメータをつけ変えたものが，上で与えた常微分方程式となるのである．

さて，ロジスティック方程式は，次のようにして求積できてしまう．まず，両辺をN²で割って，

1 N²(t)

dN dt =α

( 1 N(t) −λ

)

となるので，n(t) = ( 1

N(t) −λ )

とおくと，

dn

dt(t) = d dt

( 1 N(t) −λ

)

=− 1 N²(t)

dN

dt =−αn(t) となるから，n(t) = Ce^−αtとなって，初期値N(0) =N₀として，

N(t) = 1

λ+(

N₀⁻¹−λ) e^−αt が一般解として得られる．

この解を図示すると，以下のようになる．

1 λ

0 t

N(t)

常微分方程式として，あるいは力学系として，これを見る場合には，もう少し広い視点で見る必要も 出てくるが，ここでは，初期値N₀が0< N₀ <1/λを満たす場合を考えよう．図示したような解を差分 法によって得ることが出来るかどうかを考えたいのである．

まず，左辺を前進差分で近似してみることにしよう．

N(t+ ∆t)−N(t)

∆t =αN(t)(1−λN(t)) である．これを漸化式の形に変形すると，

N_n+1 = (∆t)αN_n(1−λN_n) +N_n

となる．この差分スキームは前進差分による近似という極めて簡単な操作でできたものであるが，その 振る舞いは極めて興味深い．もちろん，上の図のように振舞う場合もあるが，それ以外に，周期的に振 舞う場合と，カオス的に振舞う場合とが現れるのである．

筆者の数値実験では，α=λ= 1として，N₀ = 0.1とした場合，∆tの値が1.0あたりでは正常にN = 1 に漸近するが，1.5あたりからN = 1の回りで振動しながら漸近するようになり，2.5あたりでは2Step ごとに周期的に振動するようになり，3.0に近づくにつれて，カオス的な振る舞いをするようになる．3.0

以上になると，ある時点で値が負になって，−∞に発散してしまう．

もう少し振る舞いの良い差分スキームを構成したい．そのような方法のひとつとして，例えば，ある 種のゲージ変換不変性を利用するという方法がある．このような方法は，広田良吾によって研究され，

広田の双線形化法などと呼ばれている．

まずN(t) =g(t)/f(t)とおく．すると，もとのロジスティック方程式は，

1 f

dg dt +g d

dt (1

f )

=αg f

(

1−λf g

)

となる．さらに変形して，

1 f

dg dt − g

f² df dt =αg

f (

1−λg f

)

となり，両辺にf²をかけることで，

fdg dt −gdf

dt =αg(f−λg) となる．ここで，f, gの微分項をそれぞれ前進差分で差分化すると，

fg(t+ ∆t)−g(t)

∆t −gf(t+ ∆t)−f(t)

∆t =αg(t)(f(t)−λg(t)) より，

g(t+ ∆t)f(t)−g(t)f(t+ ∆t) =α(∆t)g(t)(f(t)−λg(t))

を得る．さらにここで，f(t)7→f(t)h(t), g(t)7→g(t)h(t)なる変換（ある種のゲージ変換）を施すと，

{g(t+ ∆t)f(t)−g(t)f(t+ ∆t)}h(t+ ∆t)h(t) = α(∆t)g(t)(f(t)−λg(t))(h(t))²

となる．この表式からわかることは，f, gをともに前進差分で置き換えたものは，上で与えたゲージ変換 に対して不変ではないということである．もし右辺のh(t)²がh(t+ ∆t)h(t)であれば，上のゲージ変換 に対して不変であることになるわけである．そこで，そうなるように，差分近似式を書き変えよう．以 下では，3つの方法を紹介する．

まず，最初の方法は，

g(t+ ∆t)f(t)−g(t)f(t+ ∆t) = α(∆t)g(t+ ∆t)(f(t)−λg(t))

とする方法である．（下線部の変更に注意．）これを，両辺f, gで割るなどして，Nに関する漸化式の形に 戻すと，

N_n+1 = Nn

1−(∆t)α(1−λN_n) が得られる．

第2の方法は，

g(t+ ∆t)f(t)−g(t)f(t+ ∆t) = α(∆t)g(t)(f(t+ ∆t)−λg(t+ ∆t))

とする方法である．（やはり下線部の変更に注意．）これを上と同様にN についての漸化式に戻すと，

N_n+1 = (1 + (∆t)α)Nn

1 + (∆t)α(1−λN_n)

となる．

最後に第3の方法として，かなり複雑ではあるが，

g(t+ ∆t)f(t)−g(t)f(t+ ∆t) = 1

2α(∆t) [g(t+ ∆t)f(t) +g(t)f(t+ ∆t)]−α(∆t)λg(t)g(t+ ∆t) とするものがある．（ゲージ変換不変であることの検証は読者に委ねる．）やはり，これもNについての 漸化式に戻すと，

N_n+1 =

(1 + ¹₂(∆t)α) N_n 1− ¹₂α(∆t) + (∆t)αλN_n となる．

このようにして得られた差分スキームは∆tによらず安定であることが知られています．少し複雑な差 分スキームであるが，ある種のゲージ変換不変性を利用すると，良い差分スキームがつくれることもあ るという事実は記憶にとどめておく価値はあるだろう．

もうひとつ，中心差分による方法を試してみよう．N の微分項を中心差分という近似度2の差分近似 式で近似することは，ひとつの試行として有用であり，しかも結果が面白い．

中心差分によってロジスティック方程式を書き換えると，

N(t+ ∆t)−N(t−∆t)

2∆t =αN(t)(1−λN(t)) より，

N_n+1 = 2α(∆t)N_n(1−λN_n) +N_n−1

なる漸化式が得られる．この中心差分による方法を用いるときには，最初の1ステップだけは，別の方 法で計算しておかねばならないことは，すでに注意した．ここでも，やはり前進差分によるものなどを 用いる必要がある．

実際に数値計算してみると，∆t= 0.001として前進差分により初めの1項N1を計算し，その後∆t= 1.0 として中心差分による漸化式で数値計算を行った場合，周期的な振動が見られる．ここで∆t= 0.7くら いにすると，初めN = 1に漸近していっていた解が，N = 1に近づくにつれて激しく振動し，やがて N = 0の近くまでいき，今度は再び緩やかにN = 1に漸近し，また再び激しく振動するという様子が見 られる．

すでに見てきたように，微分方程式を数値的に解く場合には，差分スキームのつくり方によって，本 来（厳密）解ではないはずの数値解が現れてしまうことがある．このような解を幻影解などと言うこと がある．このような解は，例えば熱伝導方程式のところでみた例や，常微分方程式のところでも少しあっ たように，解が最終的に発散してしまうようなものばかりではない．周期的な振動を行う解は，最終的 に+∞や−∞に発散するような解ではないし，そのような周期的な現象は，物理的に見ても比較的妥当 性のある（つまりあっても不思議ではない）現象なのである．しかし，求める微分方程式の解ではない のだから，数値解の考察には慎重さが要求されることがわかるであろう．

ここまで，差分法で常微分方程式を解くことの危うさと「限界」を強調してきたが，次節以降に見る

van der Pol方程式，ロトカ・ヴォルテラ方程式，レスラー方程式，ローレンツ方程式は，特徴的な現象

を差分法によって捉えることが可能であるという意味で，差分法の「宣伝」でもあるし，「可能性」を見 ることでもある．