先鋒-頂級種群反應(yīng)擴(kuò)散模型的平衡點(diǎn)線性穩(wěn)定性

曹佳敏, 溫俊浩, 翁佩萱

(華南師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院,廣州 510631)

在生物圈中,不同的種群有著不同的特性描述.“先鋒種群”是指在種群密度較低時(shí)生長旺盛,隨著種群密度的不斷增加,競爭加劇而導(dǎo)致生長緩慢的一類種群,其適應(yīng)度函數(shù)在種群密度較低時(shí)取得最大值,隨著種群密度的不斷增加而持續(xù)單調(diào)下降,該種群的適應(yīng)度函數(shù)只有唯一的零點(diǎn).在自然界中這個(gè)種群扮演的是“拓荒者”的角色.例如：松樹、北美鵝掌楸等就是符合這種生長規(guī)律的物種.“頂級種群”是另一類種群,它在種群密度較低或較高時(shí),相應(yīng)的適應(yīng)度函數(shù)的值均小于零,而其適應(yīng)度函數(shù)的最大值一般在適應(yīng)度函數(shù)的2個(gè)零點(diǎn)之間取得.滿足這種生產(chǎn)規(guī)律的物種有:橡樹、楓樹等.

研究多個(gè)物種之間相互作用的動力學(xué)規(guī)律是生物數(shù)學(xué)的重要目標(biāo).一般說來,種群之間相互作用的動力學(xué)模型是以一個(gè)微分系統(tǒng)呈現(xiàn)的. 如果在考慮時(shí)間因素和種群相互作用的同時(shí)也考慮種群在空間的擴(kuò)散, 那么常見的一類模型是反應(yīng)擴(kuò)散方程, 更詳細(xì)的介紹可見文獻(xiàn)[1]、[2]. 本文研究的是如下關(guān)于先鋒種群和頂級種群相互作用的反應(yīng)擴(kuò)散方程:

對系統(tǒng)(1)進(jìn)行無量綱化,可以得到:

根據(jù)前面的描述, 先鋒種群和頂級種群的適應(yīng)度函數(shù)f和g具有以下性質(zhì)：

(i)f′(z)<0且存在唯一的z0>0, 使得f(z0)=0.

相應(yīng)的圖像可以參考圖1. 我們可以舉出一些滿足性質(zhì)(i)和(ii)的適應(yīng)度函數(shù)例子[3-5],例如：f(u)=er(1-u)-a,g(u)=uer(1-u)-a.

圖1 先鋒種群和頂級種群的適應(yīng)度函數(shù)圖像

如果不考慮空間擴(kuò)散因素,則這個(gè)系統(tǒng)對應(yīng)于一個(gè)常微分方程組,其平衡解的穩(wěn)定性與Hopf分支已有較多研究[6-8]. 就我們所知,對于反應(yīng)擴(kuò)散系統(tǒng)(1),文獻(xiàn)[3]、[9]就平衡點(diǎn)處可能出現(xiàn)Hopf分支進(jìn)行了討論,文獻(xiàn)[8]討論了系統(tǒng)(1)的圖靈不穩(wěn)定性. 但是,系統(tǒng)(2)的平衡態(tài)穩(wěn)定性問題尚未有文獻(xiàn)全面討論. 受文獻(xiàn)[10]、[11]的方法啟發(fā),本文主要研究系統(tǒng)(2)在所有平衡點(diǎn)處的線性穩(wěn)定性.

1 系統(tǒng)平衡點(diǎn)分析

其次,可以看出,除了上面的4個(gè)平衡點(diǎn),其他的可能平衡點(diǎn)可以通過求解以下方程組得到:

容易得到方程組(3)、(4)分別有解：

由于生物意義,我們只關(guān)注系統(tǒng)的共存平衡點(diǎn),即滿足u+>0,v+>0的正平衡點(diǎn)(u+,v+). 也就是說,只關(guān)注本身是正平衡點(diǎn)的(u1,v1)和(u2,v2). 通過分析c11u+v=z0,u+c22v=w1,u+c22v=w2這3條直線的位置關(guān)系,可以把系統(tǒng)(2)存在平衡的情況分為如下4類(包括9種子情形. 具體的討論類似于沒有擴(kuò)散的常微分系統(tǒng),可參考文獻(xiàn)[12], 這里略去):

(1)在如下任何一組條件下,系統(tǒng)(2)只有4個(gè)邊界平衡點(diǎn), 沒有共存平衡點(diǎn):

(2)在如下任何一組條件下,系統(tǒng)(2)有4個(gè)邊界平衡點(diǎn)和1個(gè)共生平衡點(diǎn)(u2,v2):

(3)在以下任何一組條件下,系統(tǒng)(2)有4個(gè)邊界平衡點(diǎn)和1個(gè)共生平衡點(diǎn)(u1,v1):

(4)在以下任何一組條件下,系統(tǒng)(2)有4個(gè)邊界平衡點(diǎn)和2個(gè)共生平衡點(diǎn)(u1,v1)和(u2,v2):

2 平衡點(diǎn)的線性穩(wěn)定性

我們將討論系統(tǒng)(2)各個(gè)平衡點(diǎn)的線性穩(wěn)定性. 首先,假設(shè)(u0,v0)是系統(tǒng)的任意一個(gè)平衡點(diǎn), 把系統(tǒng)(2)在(u0,v0)點(diǎn)處線性化得到線性化系統(tǒng)：

考慮線性方程(6)形如

方程(8)具有非零解當(dāng)且僅當(dāng)以下行列式為零成立:

其中

A(,σ,u0,v0)=+d1σ2-f(a)-c11u0f′(a),

B(,σ,u0,v0)=+d2σ2-g(b)-c22v0g′(b),

整理上面的式子,可以得到特征方程:

A(,σ,u0,v0)B(,σ,u0,v0)-u0v0f′(a)g′(b)=0.

(9)

下面將依據(jù)方程(9)的特征根分布情況來討論系統(tǒng)各個(gè)平衡點(diǎn)的線性穩(wěn)定性.

定理1 (0,0)點(diǎn)為線性不穩(wěn)定的.

證明取(u0,v0)=(0,0), 這時(shí)a=0,b=0, 代入方程(9)得(+d1σ2-f(0))(+d2σ2-g(0))=0.由f的性質(zhì)知f(0)>0,從而可以取σ>0充分小使得=-d1σ2+f(0)>0, 即存在正的特征根, 故(0,0)為線性不穩(wěn)定的. 證畢.

(+d1σ2-z0f′(z0))(

(+d2σ2-w1g′(w1))=0.

(+d2σ2-w2g′(w2))=0.

定理5 在子情形(2b)、(4b)下,正平衡點(diǎn)(u2,v2)為線性漸近穩(wěn)定的;在子情形(2a)、(4a)下，(u2,v2)為線性不穩(wěn)定的.

證明取(u0,v0)=(u2,v2),有a=z0,b=w2, 代入方程(9)得:

u2v2f′(z0)g′(w2)=0.

(10)

令α=α(σ)和β=β(σ)如下：

d1c22v2g′(w2))σ2+(c11c22-1)u2v2f′(z0)g′(w2).

由于f′(z0)<0,g′(w2)<0,可得

Δ=α2-4β=[(d1+d2)σ2-(c11u2f′(z0)+

c22v2g′(w2))]2-4[d1d2σ4-(d2c11u2f′(z0)+

d1c22v2g′(w2))σ2+(c11c22-1)u2v2f′(z0)g′(w2)]=

(d1-d2)2σ4+2(d2-d1)σ2(c11u2f′(z0)-

c22v2g′(w2))+(c11u2f′(z0)-c22v2g′(w2))2+

4u2v2f′(z0)g′(w2)=[(d1-d2)σ2-(c11u2f′(z0)-

c22v2g′(w2))]2+4u2v2f′(z0)g′(w2)>0.

(11)

首先考慮子情形(2b)和(4b), 這時(shí)c11c22>1. 由于f′(z0)<0,g′(w2)<0, 故對任意σ, 總有α>0,β>0.由一元二次方程的求根公式可知2+α+β=0的根均為負(fù)數(shù),從而(u2,v2)點(diǎn)為線性漸近穩(wěn)定的.

其次考慮子情形(2a)和 (4a), 這時(shí)有c11c22<1, 同樣地有α>0.當(dāng)σ=0 時(shí),β=β(0)=(c11c22-1)u2v2f′(z0)g′(w2)<0.此時(shí)方程的根為一正一負(fù), 故總可以找到σ使得方程有正根, 從而(u2,v2)點(diǎn)為線性不穩(wěn)定的. 證畢.

從第1節(jié)可知平衡點(diǎn)(u1,v1)在子情形(3a)、(3b)、(4a)、(4b)下出現(xiàn). 為了討論平衡點(diǎn)(u1,v1)的穩(wěn)定性,取(u0,v0)=(u1,v1),有a=z0,b=w1, 代入式(9)得:

u1v1f′(z0)g′(w1)=0.

(12)

令p=p(σ)和q=q(σ)如下：

那么,式(12)被整理為等價(jià)形式

2+p+q=0.

(14)

方程(14)的根的判別式記為Δ(σ),則可以算出類似于式(11)的計(jì)算,只要把其中的g′(w2)換成g′(w1), 即可得到Δ=Δ(σ)的公式：

c22v1g′(w1))]2-4[d1d2σ4-(d2c11u1f′(z0)+

d1c22v1g′(w1))σ2+(c11c22-1)u1v1f′(z0)g′(w1)]=

[(d1-d2)σ2-(c11u1f′(z0)-c22v1g′(w1))]2+

4u1v1f′(z0)g′(w1).

(15)

下面2個(gè)定理的討論都將依據(jù)特征方程(14)的根的實(shí)部符號來闡述平衡點(diǎn)(u1,v1)的線性穩(wěn)定性.

定理6 在子情形(3b)、(4b)下,正平衡點(diǎn)(u1,v1)為線性不穩(wěn)定的.

證明在子情形(3b)、(4b)下,有c11c22>1. 則

Δ(0)=(c11u1f′(z0)-c22v1g′(w1))2+

4u1v1f′(z0)g′(w1)=(c11u1|f′(z0)|+

c22v1g′(w1))2-4u1v1|f′(z0)|g′(w1)=

(c11u1|f′(z0)|)2+(c22v1g′(w1))2+

2c11c22u1v1|f′(z0)|g′(w1)-4u1v1|f′(z0)|g′(w1)≥

4(c11c22-1)u1v1|f′(z0)|g′(w1)> 0.

另一方面, 有

q(0)=(c11c22-1)u1v1f′(z0)g′(w1)<0.

無論p(0)為任何實(shí)數(shù),根據(jù)一元二次方程的求根公式,方程2+p(0)+q(0)=0的根總是一正一負(fù). 也就是說,在子情形(3b)、(4b)下,總可以找到σ使得方程2+p+q=0有正根,從而(u1,v1)為線性不穩(wěn)定的. 證畢.

為了在下面的定理中討論子情形(3a)和(4a), 我們先做一些準(zhǔn)備工作. 注意這時(shí)有c11c22<1.為了方便討論,引進(jìn)一些記號. 令

顯然有ξ≥0,ζ>0,μ<0,ν<0,γ>0.

根據(jù)以上定義給出下面的引理.

引理1 在子情形(3a)、(4a)下,以下結(jié)論成立：

(w1)當(dāng)θ>0 時(shí),η>(d1-d2)ζ;

(w2)當(dāng)θ=0 時(shí),η=(d1-d2)ζ;

(w3)當(dāng)θ<0 時(shí),η<(d1-d2)ζ.

把這些記號代入式(13)和式(15)中得到:

p(ξ)=(d1+d2)ξ-θ,q(ξ)=d1d2ξ2-ηξ+γ,

Δ(ξ)=[(d1-d2)ξ-μ]2+ν=

(d1-d2)2ξ2+2(d2-d1)μξ+μ2+ν.

先分析p(ξ)=0. 由于p′(ξ)=d1+d2>0, 那么對?ξ≥0有p(ξ)單調(diào)遞增,對p(ξ)有:

(p1)當(dāng)θ<0 時(shí),對?ξ≥0有p(ξ)>0;

(p2)當(dāng)θ=0 時(shí),對?ξ≥0有p(ξ)≥0;

最后分析q(ξ)=0,即

d1d2ξ2-ηξ+γ=0.

(16)

我們分以下4種情況來分析：

顯然在[0,ξ1)∪(ξ2,+∞)上有q(ξ)>0, 在(ξ1,ξ2)上有q(ξ)<0.

定理7 在子情形(3a)、(4a)下, 如果滿足下列條件之一：

(1)d1≤d2,θ<0;

則(u1,v1)為線性漸近穩(wěn)定的.

證明由θ<0和(p1)推知:在定理的任何一個(gè)條件下,總有p(ξ)>0對ξ≥0.

先證明在條件(1)下(u1,v1)為線性漸近穩(wěn)定的. 由d1≤d2,θ<0 和(w3)推知η<0, 從而由(q4)推知q(ξ)>0對ξ≥0.由方程(13)和一元二次方程的求根公式知道:當(dāng)Δ>0時(shí),方程(13)的根為兩負(fù)根; 當(dāng)Δ=0時(shí),方程(13)的根為1,2=-p/2;當(dāng)Δ< 0時(shí),方程(13)有一對共軛復(fù)根具有負(fù)實(shí)部, 故(u1,v1)為線性漸近穩(wěn)定的.

證畢.

定理8 在子情形(3a)、(4a)下, 如果滿足下列條件之一：

(1)θ>0;

則(u1,v1)為線性不穩(wěn)定的.

注1 注意定理7中的條件(1)和(w2)、(w3)推出η≤(d1-d2)ζ≤0, 這時(shí)定理8中的條件(1)、(2)自然都不成立.也就是說,定理7中的任何一組條件都是同時(shí)否定定理8中的條件(1)、(2)的.

3 數(shù)值模擬

考慮如下模型：

顯然,系統(tǒng)(17)中f(z)=z0-z,g(w)=-(w1-w)×(w2-w).通過計(jì)算得f′(z0)=-1,g′(w1)=w2-w1>0,g′(w2)=w1-w2<0,g′(w*)=0, 這里w2>w1>0,w*=(w1+w2)/2.

首先考慮子情形(3a),下面總?cè)?/p>

將式(18)代入系統(tǒng)(17)有

經(jīng)過計(jì)算得(u1,v1)=(4/3,1/3).

系統(tǒng)(19)在(u1,v1)=(4/3,1/3)處的線性化方程為

取d1=1,d2=1/2, 此時(shí)θ=-5/12,η=-1/12,γ=1/2,顯然滿足定理7中的條件(2),這時(shí)有(u1,v1)=(4/3,1/3)是線性漸近穩(wěn)定的,模擬結(jié)果見圖2.

對于原系統(tǒng)(19)，仍然取d1=1,d2=1/2,此時(shí)根據(jù)模擬的結(jié)果(圖3),我們猜測(u1,v1)=(4/3,1/3)為局部漸近穩(wěn)定的.

這一節(jié)只給出了定理7中子情形(3a)下條件(2)的數(shù)值模擬,對于定理7中子情形(3a)下的條件(1),我們只需隨意取定一組d1、d2,滿足0≤d1≤d2,其余的參數(shù)如式(18)所取,對于定理7中子情形(4a)的討論只需取z0=1,w1=3/2,w2=2,c11=1/4,c22=1/2.而d1,d2同子情形(3a)下的取法,均可以進(jìn)行數(shù)值模擬,將不再贅述.

參考文獻(xiàn)：

[1] Shigesada N, Kawasaki K. Biology invasions: Theory and practice[M]. Oxford: Oxford University Press, 1997.

[2] Murray J D. Mathematical biology: I, II[M]. New York: Spriner-Verlag, 2002.

[3] Selgrade J F, Namkong G. Stable periodic behavior in a pioneer-climax models[J]. Nature Resource Modeling, 1990, 4:215-227.

[4] Ricker W E. Stock and recruitment[J]. Journal of the Fisheries Research Board of Canada, 1954, 11:559-623.

[5] Franke J E, Yakubu A A. Pioneer exclusion in a one-hump discrete pioneer-climax competing system[J]. Journal of Mathematical Biology, 1994,32: 771-787.

[6] Summer S. Hopf bifurcation in pioneer-climax species models[J].Mathematical Biosciences, 1996, 137: 1-24.

[7] Summer S. Stable periodic behavior in pioneer-climax competing species models with constant rate forcing[J]. Nature Resource Modeling, 1998,11: 155-171.

[8] Buchanan J R. Turning instability in pioneer-climax species interactions[J]. Mathematical Biosciences, 2005, 194: 126-199.

[9] Liu J X, Wei J J. Bifurcation analysis of a diffusive model of pioneer and climax species interaction[J]. Advance in Difference Equation, 2001, 24: 11-52.

[10] 馬知恩, 周義倉. 常微分方程定性與穩(wěn)定性方法[M]. 北京：科學(xué)出版社, 2004.

[11] 吳敏, 翁佩萱. 具有階段結(jié)構(gòu)的多時(shí)滯SIR擴(kuò)散模型的穩(wěn)定性[J]. 華南師范大學(xué)學(xué)報(bào):自然科學(xué)版, 2013, 45(2): 20-23.

Wu M, Weng P X. Stability of a stage-structured diffusive SIR model with delays[J]. Journal of South China Normal University: Nature Science Edition, 2013, 45(2): 20-23.

[12] Buchanan J R. Asymptotic behavior of two interacting pioneer-climax species[J].Fields Institute Communication, 1999, 21: 51-63.