ロジスティック曲線

ロジスティック曲線
KENZOU
2005 年 2 月 13 日（2 月 14 日:米国人口統計追記）
最近、日本の出生比率が極端に下がってきたと新聞を賑わしています。日本人類存亡の危機に瀕しているとは
少し大袈裟ですが、全世界で比類ない独特の日本文化を継承する子孫が少なくなってくるというのは寂しいこと
です。さて、以下は生物の自然増加現象等を説明する際によく引き合いにだされるロジスティック曲線について
のレポートです。これは非線形微分方程式の解となります。
ロジスティック曲線
1
1.1
マルサスのモデルとロジステック方程式
生物の増える速さは現在の生物の個体数に比例する、これをマルサスのモデルという。
dx
= kx,
dt
k=
1 dx
(k > 0 は常数で個体数増加率)
x dt
(1)
マルサスのモデルは生物個体が指数関数的に増大することを表すが、これは現実とは合わない。実際
には、増加率そのものが個体数の大きさによって変化する。個体数が増加すれば個体数増加を抑制する
力（えさが不足するとか）が加わると考えた方が現実的である。即ち、ある程度までは急激に増えて、
増えすぎると抑制しようとする力が働いてきて、最後にある値で安定してくる。このことを考慮して
オランダの数理生物学者 Verhulst（ヴェアファルスト）が 1987 年に人口増加の方程式を提案した。こ
れが有名なロジステック方程式である。人口の上限を Q とし、その時点での人口を P とすると
µ
¶
dP
Q−P
= kP
(2)
dt
Q
この微分方程式は P について１次ではないので、非線形微分方程式である。P = Q となると右辺はゼ
ロ、つまり人口増加率はゼロとなるから人口の増加が止まる。
さて、これからこのロジスティク方程式を解いていくが、(7) の方程式を次のように簡略化して書く
ことにする1 。
dx
= R(x)x,
R(x) = a − bx (a > 0, b > 0)
(3)
dt
a は増殖率またはマルサス係数であり、環境的影響を受けない成長率である。 b は混雑定数 (crowdness
constant) と呼ばれる。
1
式 (7) との関係は Q = a/b, k = a で結ばれる。Q は環境容量とも呼ばれ、すぐ後でその意味が明らかになる。
1
1.2
ロジスティック方程式を解く
ロジスティック方程式を解く。
dx
= (a − bx)x (a > 0, b > 0)
dt
(4)
ロジスティック方程式 (4) を Mathematica で解いてみる。初期条件は t = 0 で x[0] = x0 とする。
In[1] := sol1 = x[t]/.DSolve[{x′ [t] == (a − bx[t]) x[t], x[0] == x0}, x[t], t}][[1]]
Out[1] =
aeat x0
a − bx0 + beat x0
解 sol を整理してみよう。分母分子を beat x0 で割ると分母は
In[2] := a − bx0 + beat x0)/(beat x0)//Expand
−at
Out[2] = 1 − e−at + aebx0
となる。従って解は次の形に整理される2 。
x[t] =
³a´
b
1
,
1 + c e−at
但し c =
a
−1 bx0
(5)
ここで解 sol1 の t → ∞ の極限をとってみよう。(5) より目視で
lim x[t] = a/b
t→∞
at
ae x0
となるのは明らかだが、ここでは Out[1] = a−bx0+be
at x0 の解から求めてみる。 Mathematica でこのよ
うな式の Limit を計算する場合、定数 a, b, x0 の符号を明確にしてやらないと計算してくれない。そこ
で組み込み関数 Assuming を使って符合を明確にして、Limit 計算をする。
In[3] := Assuming[a > 0 && b > 0 && x0 > 0, Limit[sol, t → ∞]]
a
Out[3] =
b
以上のことから、x(個体数) はどんなに増えても定数 a/b を超えることがない。a/b は環境容量とも呼
ばれる。
1.3
ロジスティク曲線を描く
f (t) =
1
のそれぞれ定数 a, c を変化させた場合のグラフのアウトラインを調べてみよう。
1 + c e−at
(1) c を一定にして a を変えていった場合：f (t) =
1
1 + 10 e−at
−→ 曲線の立ち上がりが a の増加につれて急峻になってくる。
1
・
・図 (1)
Plot[Table[ 1+10e
−at , {a, 20, 100, 20}]//Evaluate, {t, −.15, .4}]・
2
この形は Fermi 分布関数（f (ϵ) =
1
1+e(ϵ−µ)/kB T
) によく似ていることは興味深い。
2
(2) a を一定にして c を変えていった場合：f (t) =
1
1 + c e−20t
−→ 曲線の y 切片が c の増加につれて下方に向かう。
・
・図 (2)
Plot[Table[ 1+ce1−20t , {c, 1, 10, 1}]//Evaluate, {t, −.15, .4}]・
(3) c がマイナス、即ち x0 > a/b の場合：f (t) =
1
1 − 10 e−20t
−→ グラフの様相は一変する。
・
・
・図 (3) 1
Plot[ 1−10e
−20t , {t, −/15, .4}]
このケースは初期値 (x0) が環境容量 (a/b) より大きいということで、物理的には意味がない。
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
10
5
-0.1
0.1
0.2
0.3
0.4
-5
-10
-0.1
0.1
0.2
0.3
-0.1
0.4
図 1: c 一定で a を変えた場合
0.1
0.2
0.3
0.4
図 2: a 一定で c を変えた場合
図 3: c < 0 の場合
さて、ロジスティック方程式の解 sol1 のグラフを描くこととする。ここで a = 20, x0 = 10 とし、b
を 0.1 から 0.15 まで 0.01 刻みに変えた。
In[4] := sol2 = sol1/.a → 20/.x0 → 10;
Plot[Table[sol2, {b, 0.1, 0.15, 0.01}]//Evaluate, {t, −.2, .5}, GridLines → Automatic, Frame → True]
200
150
100
50
0
-0.2
-0.1
0
0.1
0.2
0.3
図 4: ロジスティク曲線
3
0.4
0.5
1.4
生物が一定の割合で捕獲される場合のロジスティク曲線
捕獲率を α(0 < α < 1) とすると、ロジスティック方程式は
dx
= (a − bx)x − αx (a > 0, b > 0, 0 < α < 1)
dt
(6)
となる。この微分方程式の解を sol3 としておく。
In[5] := sol3 = x[t]/.DSolve[x′ [t] == (a − bx[t])x[t] − cx[t], x[0[== x0, x[t], t][[1]]
a = 20, b = 0.1, c = 0.9, x0 = 10 として sol2, sol4 の曲線を描いてみよう。
In[6] := sol2 = sol1/.a → 20/.x0 → 10/.b → 0.1;
sol4 = sol3/.a → 20/.x0 → 10/.b → 0.1;
Plot[sol2, sol4, t, −0.2, 0.5, GrifLines → Automatic, Frame → True]
200
150
100
50
0
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
図 5: 一定の捕獲がある場合との比較
sol3 の t → ∞ の極限値を求めてみると
In[7] := Assuming[(a − c) > 0&&mathsf b > 0&&x0 > 0&&c > 0, Limit[sol3, t → ∞]]
a−c
Out[7] =
b
となって、一定の捕獲率がない場合の環境容量 a/b より小さい値となることが分かる。
ロジスティック曲線と米国の人口統計
2
最後の締めくくりとして Verhulst モデルの人口増加の方程式（ロジスティック方程式）の威力がど
れほどのものか、1790 年から 1930 年までの米国の人口実測値とロジスティック方程式からの予測値を
示しておく3 。ここで物理的意味が分かりやすい (7) を使うことにする。
µ
¶
dP
Q−P
= kP
(7)
dt
Q
3
http://www.fit.ac.jp/elec/lab/kawabatalab
4
6 t = 0 での初期人口数）
パラメータは P 0 = 3.9 × 10（
, k = 0.3134, Q = 197 × 106 （環境容量：人口増
加の最大値）とする。
In[8] := Psol1 = P[t]/.DSolve[P′ [t] == kP[t](1 − P[t]/Q), P[0] == P0, P[t], t][[1]];
Psol2 = Psol1/.k → 0.3134/.Q → 197 ∗ 106 /.P0 → 3.9 ∗ 106 ;
data = Table[Pso2, t, 0, 14, 1]
Out[8] = {3.9 × 106 , 5.29687 × 106 , 7.17547 × 106 , 9.68667 × 106 , 1.30165 × 107 , 1.73846 × 107 ,
2.30353 × 107 , 3.02133 × 107 , 3.91253 × 107 , 4.98799 × 107 , 6.24218 × 107
7.64778 × 107 , 9.15461 × 107 , 1.06949 × 108 , 1.21946 × 108
Verhulst モデルは 140 年余りにわたって正確に米国の人口を予測していることは注目に値する。しか
しこのモデルは 1930 年以降も米国の総人口を予測できるかとなると、これは怪しい。というのは、米
国の人口はこのモデルの人口飽和値 197 × 106 （2090 年頃）を既に上回っている。
1.2·10
8
2·10
8
8
1·10
8
1.5·10
7
8·10
8
7
1·10
6·10
7
4·10
7
5·10
7
2·10
0
0
0
2
4
6
8
10
12
0
14
図 6: 1790 年から 1930 年
5
10
15
20
図 7: 1790 年∼2090 年
表 1: Verhulst モデルと米国の人口統計
年度
米国の人口 (×106 )
予測値 (×106 )
1790
1800
1810
1820
1830
1840
1850
1860
1870
1880
1890
1900
1910
1920
1930
3.9
5.3
7.2
9.6
12.9
17.1
23.2
31.4
38.6
50.2
62.9
76.0
92.0
106.5
123.2
3.9
5.3
7.2
9.7
13.0
17.4
23.0
30.2
38.1
49.9
62.4
76.5
91.6
107.0
122.0
5
25
30