ロジスティック成長する個体群に対して 空間非一様性が与える影響
永原 健大郎 （Kentaro Nagahara） 東京工業大学大学院 理工学研究科 数学専攻
e-mail: [email protected]
1 Introduction
環境資源の配置が，ある種の生物の個体数にどのように影響するかという問 題は，保全生物学において重要な問題として残されている[3]．本講演では，
総資源量が一定という制約条件の下で，生物の個体数が最大となる資源配置 は具体的にどのような形で与えられるのか，という問題について論ずること を目的とする．
始めに，本研究の背景に触れるため，Verhulst [6]が提唱した次の方程式 を見てみよう．
d
dtw(t) = w(t){k−w(t)}
これはロジスティック方程式と呼ばれる生物の増殖を表す数理モデルで，人 口の無制限な成長は不可能で，増加率はやがて鈍化し，最終的には定常的に なるという考え方を反映させている．ここで，t ∈ R+は時間，w(t) ∈ Rは 生物種の個体数，k ∈[0,1]⊂Rは内的自然増加率と呼ばれる定数で，その生 物が潜在的に持っている最大の繁殖増加率を表している．
初期値がw(0) >0を満たしているとする．この方程式は変数分離系であり 具体的に解くことができ，
w(t) = k
1 + ( k
w(0) −1 )
e−kt
if (0< k ≤1),
w(t) = 1 t+ 1
w(0)
if (k = 0)
と表される．0< k≤1のとき，解はw(t)→k (t→ ∞)となり，個体数はk まで増殖するが，k = 0のとき，解はw(t)→0 (t → ∞)となり，生物群は絶 滅してしまう．(kは一般にk∈Rとしても同様の議論ができるが，今回は簡 単のために0≤k ≤1として扱った．)
次に，生物が生息域である有界領域内をランダムに移動することを考慮 に入れる．ただし，この生息域は環境が一様ではなく，内的自然増加率が空
間変数xに依存する関数k =m(x)によって表されるものとする．言い換え れば，m(x)は場所xに生物が増殖するのに適した資源がどの程度置かれて いるかを表しており，この関数m(x)を環境容量と呼ぶ．このm(x)は非負値
Lebesgue可測関数であり，次の条件を満たすものとする．
(M1) 0≤m(x)≤1 L-a.e.
(M2) |Ω1|∫
Ωm(x)dx =δ.
ここで，m(x)が0に近い場所は生物にとって住みにくい環境であり，m(x) が1に近いところは生物が住みやすい環境であることを表している．また，δ は総資源量を表す定数で，0< δ < 1を満たす．すなわち，m ≡0 in Ωとい う資源が全くない状況，またはm ≡ 1 in Ωという資源がすべての領域で豊 富にある状況は考えないものとする．以降このδを固定した状況を考える．
以上の状況の下，場所x，時刻tにおける生物種の個体密度分布vm,ε(x, t) は，次の空間非一様性のある反応拡散ロジスティック方程式(あるいはFisher- KPP方程式)と呼ばれる数理モデルを用いて表される[5]．
∂vm,ε
∂t =ε2∆vm,ε+m(x)vm,ε−vm,ε2 ((x, t)∈Ω×R+),
∂vm,ε
∂n = 0 ((x, t)∈∂Ω×R+),
vm,ε(x,0)≥0, vm,ε(x,0)̸≡0 (
x∈Ω) ,
(1)
m∈X :=
{
m∈L∞(Ω)
0≤m(x)≤1, 1
|Ω|
∫
Ω
m(x)dx=δ }
.
Ω⊂Rnは生物の生息域に対応し，滑らかな境界∂Ωを持つ．(1)の第2式に
あるNeumann境界条件は，生物が∂Ωを超えて生息域Ωの外に出ないこと
を意味する．nは∂Ω上の外向き単位法線ベクトルである．また，ε > 0は 正の定数で，拡散係数と呼ばれる生息域内を生物が移動する速さを表すパラ メータである．εが大きければ，生息域内を生物が速いスピードで移動し，ε が小さければ生物は遅いスピードで移動することに対応している．初期条件 vm,ε(x,0)∈C(Ω)は観測を開始した時点での生物の生息域内での分布を表す．
また，X ⊂ L∞(Ω)は，条件(M1),(M2)を満たしている関数m(x)全体を表 している．なお，m(x)が定数ならば，冒頭で紹介したロジスティック方程式 と同様の計算ができる．
ここで，(1)の数理モデルで記述される生物個体群が絶滅(||vm,ε(·, t)||L∞(Ω)→ 0 (t → ∞))せず増殖できるかどうかが問題となるが，これに対しては，
Cantrell-Cosner([1],[2])による結果が知られている．彼らの一連の研究結果 により，m ∈ X のとき，(1)には非自明な正値定常解 um,ε ∈ W2,p(Ω) ∩ C1,α(Ω) (∀p > 1)が一意的に存在し，初期値vm,ε(x,0)に依存せずt → ∞ ですべての解がum,ε(x)に一様収束することが示されている．
この結果を用いて，十分に時間tが経過した後の生物の個体数を調べるた め，(1)の定常解が満たす半線形楕円型偏微分方程式を確認しておこう．
−ε2∆um,ε =m(x)um,ε−u2m,ε (x∈Ω),
∂um,ε
∂n = 0 (x∈∂Ω) (2)
ここで，解um,ε(x)は環境容量m(x)，拡散係数εに依存することに注意 しておく．
これを受け，W.Dingら5名[3]の研究によってNet-Benefitと呼ばれる次 の量が導入された．
J(m, ε) :=
∫
Ω
[um,ε(x)−Bm(x)2] dx ただしBは非負の定数とする．J(m, ε)は，生息数∫
Ωum,ε(x)dxから資源配 置にかかるコスト∫
ΩBm(x)2dxを引いた値を表している．すなわち，J(m, ε) の値が大きい資源配置は，生物の生息数が多くなり，かつコストの少ない資 源配置となっていると言える．
では，資源の総量δ が一定の下で，J(m, ε)が最も大きくなる資源配置 m(x)とは具体的にどのような資源配置なのだろうか．これは，J(m, ε)の大 域的最大化解となるm∗(x)∈U を具体的に求める問題となる．
問題に入る前に，大域的最大化解(global maximizer)と局所的最大化解 (local maximizer)について述べる．環境容量m(x)の集合X内のJ(m, ε)の 上限
Jsup := sup
m∈X
J(m, ε)
に対し，J(m∗, ε) = Jsupを満たすm∗ ∈ Xのことを大域的最大化解という．
ここで，拡散係数εは生物種の生態によって決まる値であり，ここでは固定 されている．また，m(x)˜ ∈Xに対し，あるµ >0が存在し，
∫
Ω
|m(x)−m(x)˜ |dx < µ を満たすすべてのm∈Xに対して，
J(m, ε)≤J( ˜m, ε)
が成り立つとき，m(x)˜ を局所的最大化解という．ここで，大域的最大化解は 局所的最大化解の条件も満たしていることに注意する．W.Dingら[3]によっ
て，すべてのB ≥0，ε >0で大域的最大化解m∗(x)∈Xが存在することが 示されている．
具体的形状についてもW.Dingら[3]によって部分的に解答が得られてお り，B ≥0が十分大きい場合には，m(x)≡δがJ(m, ε)の大域的最大化解で あることが証明されている．しかしBが小さいときは，J(m, ε)の値が(1)の 定常解の空間変数xに対する積分∫
Ωum,ε(x)dxに大きく依存するため解析が 難しく，具体的な形状はわかっておらず，研究の及んでいない部分が多い．
そこでDingらは[3]において数値解析を行い，B = 0のとき，すなわち J(m, ε) =∫
Ωum,ε(x)dxのとき，J(m, ε)の大域的最大化解はbang-bang性と 呼ばれる性質を持つという予想を立てた．ここで，環境容量m(x)がbang- bang性を持つとは，m(x)がLegbesgue可測集合E ⊂Ω上でm(x) = 1を満 たし，Ω\E ⊂Ω上でm(x) = 0を満たすことをいう．
本講演では，B = 0，εを固定したとき，環境容量m(x)を摂動させ，m(x) がJ(m)の停留点となる条件を調べた結果を報告し，そこから明らかになっ たJ(m)の大域的最大化解が持つ条件について述べる．次に，Ω = (0,1)と し，これまでとは逆に環境容量m(x)を固定し，εを動かすことを考える．考 える環境容量m(x)がbang-bang性を持つ場合に，特異極限問題の解と比較 することにより∫
Ωu(x)dxの拡散係数ε = 0での漸近展開が得られたことを 説明し，これらを用いた発展的な議論を加えて発表したいと思う．
2
大域的最大化解の性質
定理 2.1. m∗(x)∈Xを大域的最大化解とする．このとき，
ess sup
x∈Ω
m∗(x) = 1, or
ess inf
x∈Ω m∗(x) = 0.
が成り立つ．
証明. 概略を説明する．まず，X0 ⊂Xを次のように定義する．
X0 :=
{
m ∈L∞(Ω)
0< m(x)<1, 1
|Ω|
∫
Ω
m(x)dx=δ }
.
∀m0(x)∈X0に対し，µ >0，g ∈L∞(Ω); ∫
Ωgdx= 0を，m :=m0+µg ∈X0 となるようにとる．このmに対して定まる(2)の解u(x)は，
∫
Ω
u(x)dx=
∫
Ω
u0(x)dx+µ
∫
Ω
u1(x)dx+o(µ)
の形に漸近展開することができる．ここで，u0(x)及びu1(x)は，それぞれ次 の方程式の解である．
−ε2∆u0 =m(x)u0−u20 (x∈Ω),
∂u0
∂n = 0 (x∈∂Ω) (3)
−ε2∆u1 = (m(x)−2u0)u1+gu0 (x∈Ω),
∂u1
∂n = 0 (x∈∂Ω) (4)
ここで，第一変分∫
Ωu1(x)dxに注目すると，
∫
Ω
u1(x)dx = 0⇔m0 ≡δ を示すことができる．更に，m0 ≡ δは∫
Ωudxの大域的最小化解であること も示され，したがって定理の結論が導かれる．
3 J(m, ε)
の
εに対する漸近展開
方程式(2)から，次の結果が導かれる．
ε2
∫
Ω
|∇u|2 u2 dx=
∫
Ω
mdx−
∫
Ω
udx.
定理2.1よりも詳しい大域的最大化解m∗(x)の形状を調べるためには，この 左辺が最大となるようなm∗ ∈Xを探すことが重要となるが，その解析は容 易ではない．しかし，mがC1級と仮定した場合，
∫
Ω
|∇u|2dx <
∫
Ω
|∇m|2dx
となることが知られている[4]．したがってm(x)の勾配の平均とu(x)の勾配 の平均は重要な関係を持っていることが伺える．そこで，Dingら[3]の予想 をもとに，環境容量m(x)が不連続ではあるが最も変化が激しいbang-bang 性を持つときのJ(m, ε)について詳しく調べようと考えた．以降，次の状況 を考える．
Ω = (0, a)とし，環境容量m1(x)を次のように定義する．
m1(x) = {
0 (x∈(0, x0)), 1 (x∈(x0, a)),
この時，J(m1, ε)のεに関する漸近展開を調べた結果，次の式を得た．
定理 3.1.
J(m1, ε) =δ+C1ε+C2K6(sin π 12)1
x60ε7+o(ε7), (5)
C1 :=
√ 2√
3− (
3−
√
3 + 2√ 3
) , C2 := 12
√ 2√
3 (√
3√ 3 +√
3− (
√4
3 + 9 5
))
ここで，δはm1(x)の定義からδ =a−x0である．
証明. 概略を説明する．m1(x)の定義から，方程式(2)の解uは次を満たす．
−ε2uxx =−u2, 0< x < x0, ux(0) = 0,
−ε2uxx =u−u2, x0 < x < a, ux(a) = 0.
また，原点でC1級に接続する条件により，
u(x0+) =u(x0−), ux(x0+) =ux(x0−),
となる．ここで，変数変換と原点で接続する条件により，方程式を(0, x0)，
(x0, a)で分けて考えると，方程式(2)の解uを解くことは次の連立微分方程 式を解くことに帰着される．
−uAxx =−(uA)2, x∈(
−x0 ε ,0
) , uAx
(−x0 ε
)
= 0, uAx(0) =
√2
3(γ3−α3), uA
(−x0
ε )
=α, uA(0) = γ, uAx(x)≥0.
(6)
−uBxx = (uB)−(uB)2, x∈ (
0,1−x0 ε
) , uBx
(1−x0 ε
)
= 0, uBx(0) =
√2
3(γ3−α3), uB
(1−x0 ε
)
=β, uB(0) =γ, uBx(x)≥0..
(7)
ここで，uA(x)は(0, x0)上での解u(x)，uB(x)は(x0, a)上での解u(x)を表 す．さらに，α:=u(0)，β :=u(a)，γ :=u(x0+) = u(x0−)と約束する．
方程式から，(0, a)上でux >0が示される．さらに，(6)，(7)の第1式の 両辺にそれぞれuAx，uBx を掛けて両辺を積分すると，次の保存量が得られる．
duA dx =
√2
3((uA)3−α3)
( x∈[
−x0 ε ,0
] )
, (8) duB
dx =
√ 2 3
(
(uB)3− 3
2(uB)2−β3+3 2β2
) ( x∈
[
0,1−x0 ε
] )
(9)
また，Liang-Lou[4]の結果から，
α →0 (ε→0), (10)
β →1 (ε→0) (11)
が従う．したがって，ε→0としたとき，保存量は duA
dx =
√2 3
((uA)3)
(x∈(−∞,0]), (12) duB
dx =
√ 2 3
(
(uB)3 −3
2(uB)2 +1 2
)
(x∈[0,∞)) (13) となる．この保存量から得られる解uA, uBが極限問題の解に相当し，具体的 に解くことができる．
次に，極限問題の解と拡散係数εがとても小さいときの解とを比較する ことにより，J(m1, ε)の漸近展開を得ることを考える．保存量(8)，(9)を用 いて，次に述べる関係式を得ることができる．
√
√2α 3
x0
ε = 2K(sin π
12)−sn−1
vu uu ut1−
γ
α −(1 +√ 3) γ
α −1 +√ 3
2
,sin π 12
(14)
tanh
(1−x0
2ε )
<1− 1
4(1−β)<1 (15) ここで，
K(k) :=
∫ 1
0
√ 1
(1−ξ2)(1−k2ξ2)dξ (0≤k <1)
は，第一種完全楕円積分と呼ばれる母数kに関する関数であり，k→1のと き，K(k)→ ∞となる．また，
sn−1(y, k) :=
∫ y
0
√ 1
(1−ξ2)(1−k2ξ2)dξ (−1≤y ≤1)
は，Jacobiの楕円関数のうち，sn関数と呼ばれる関数の逆関数である．
この(14)，(15)，およびα，βの0次オーダー(10)，(11)から，次のようにα の漸近展開，及びβの評価を出すことができる．
α = 2√ 3
x20 K2(sin π
12)ε2+o(ε2), (16) 1−16 e−1−xε0
1 +e−1−xε0
< β < 1. (17) この評価を用いて，スケール変換後の方程式の解の積分との間に，
1 ε
∫ 1 0
u(x)dx =
∫ 0
−x0/ε
uA(x)dx+
∫ (1−x0)/ε 0
uB(x)dx=:IA+IB (18) が成り立つことに注目し，IA，IBについて計算する．
これらの積分を具体的に行うのは通常困難であるが，例えばIAについて，変 数変換
d
dxuA(x) = d dyuA(y)
を施すと，非積分関数が極限方程式の解で，積分区間にεの影響が現れるよ うにできる．従って保存量(12)，uA(y)の具体表示，及びαの漸近展開(16) も合わせて考慮すると，
IA=√ 6ε
∫ ∞
ε(γ3−α3)−1/6
1 t2
( 1 +
( 24√
3
x60 K6(sin π
12) +o(1) )
t6 )−1/3
dt (ε →0) と書くことができる．ここからβの評価(17)を使いつつ，更に計算を進める ことにより，
IA=
√ 2√
3−CA,6 1
x60ε6+o(ε6) (ε →0), (19) CA,6 := 12
√ 2√
3 (9
5−√ 3
)
K6(sin π 12) を得ることができる．IBについても，変数変換
d
dxuB(x) = d dyuB(y)
を施し，保存量(13)，及びβの評価(17)を用いることで，
IB = 1−x0
ε −
( 3−
√
3 + 2√ 3
)
+CB,6 1
x60ε6+o(ε6) (ε→0), (20)
CB,6 := 12√
2(3−√
3)K6(sin π 12)
を示すことができる．よって，各項の展開(19)，(20)を足し合わせ，(18)に 注意すると，(5)を得ることができる．
4
数値解析
本講演では，J(m, ε)を最大にする環境容量m∗を具体的に調べるため，射影勾 配法を実装してその解析を試みたことを発表する．射影勾配法とは，J(m, ε) のmに関する第一変分を調べ，J(m, ε)が増加する方向にmを修正する技法 である．この時，修正したmがXに入るように調整をしなくてはならない．
しかし，残念ながら射影勾配法を実装しても数値解析で分かることは停留点 の形状のみで，最適制御をいきなり調べることは不可能である．理論的な発 展と並行して研究を進めたい．
本講演では，以下に述べるXδのどのような環境容量からスタートしても，
J(m, ε)が増加する方向に環境容量を修正するアルゴリズムを実装すること
を考え，これに成功したことを発表する．X内での処理についても，講演で 触れたい．
注：このセクションでは，拡散係数εは動かさないため，J(m, ε)をJ(m) と表記する．
命題 4.1.
Xδ :=
{
m∈L∞(Ω) 1
|Ω|
∫
Ω
m(x)dx=δ }
とする．任意の初期条件m0 ∈ Xδに対し，修正量mn ∈ Xδを次のように定 めると，J(m)の増大列になる．
mn=mn−1+µln−1(∀n∈N) ただし，
µ:sufficiently small,
ln−1 ∈L∞(Ω) =un−1p1,n−1+p2,n−1 (∀n ∈N) である．ここで，
un−1 =u(mn−1),
p1,n−1は次の方程式の解
−ε2∆p1,n−1−(mn−1−2un−1)p1,n−1 = 1 x∈Ω
∂p1,n−1
∂n = 0 x∈∂Ω (21)
p2,n−1は次で与えられる定数である．
p2,n−1 = −1
|Ω|
∫
Ω
(un−1p1,n−1−2mn−1B)dx.
参考文献
[1] R. S. Cantrell and C. Cosner, Diffusive logistic equations with indefinite weights: population models in a disrupted envronments, Proc. Roy. Soc.
Edinburgh, 112A (1989), 293–318.
[2] R. S. Cantrell and C. Cosner, Spatial Ecology via Reaction-Diffusion Equations, Series in Mathematical and Computational Biology, John Wiley and Sons, Chichester, UK, 2003,
[3] W. Ding, H. Finotti, S. Lenhart, Y. Lou, Q. Ye, Optimal control of growth coefficient on a steady-state population model, Nonlinear Analy- sis: Real World Applications, 11 (2010), 688–704.
[4] S. Liang, Y. Lou,On the dependence of population size upon random dis- persal rate, Special issue on ”PDE Models from Biological Processess,”
Disc. Cont. Dynam. Sys. Series B, 17 (2012), 2771–2788.
[5] J. G. Skellam, Random dispersal in theoretical populations, Biometrika, 38 (1951), 196-218.
[6] P. F. Verhulst, Notice sur la loi que la population suit dans son ac- croissement, Correspondance Math´ematique et Physique Publi´ee par A. Qu´etelet,10 (1838), 113–121.
