一類捕食與被捕食模型的動(dòng)力學(xué)穩(wěn)定性分析

孫冬梅，黃東衛(wèi)

(天津工業(yè)大學(xué) 理學(xué)院，天津 300387)

?

一類捕食與被捕食模型的動(dòng)力學(xué)穩(wěn)定性分析

孫冬梅，黃東衛(wèi)

(天津工業(yè)大學(xué) 理學(xué)院，天津 300387)

針對(duì)一種被捕食者與兩種捕食者的密度之間相互關(guān)系,建立了描述者間動(dòng)力學(xué)演化行為的數(shù)學(xué)模型,利用Routh-Hurwitz判別法，分析了模型中各種群演化過程性態(tài)的穩(wěn)定性. 并選取模型中被捕食者密度的增長率為關(guān)鍵參數(shù)，得到了該動(dòng)力系統(tǒng)發(fā)生Hopf分岔的條件,對(duì)系統(tǒng)進(jìn)行數(shù)值模擬，計(jì)算出其Lyapunov指數(shù)譜，以及發(fā)生Hopf分岔時(shí)的參數(shù)閾值, 結(jié)果表明被捕食者密度參數(shù)對(duì)系統(tǒng)性態(tài)變化起到了關(guān)鍵的控制作用，有助于認(rèn)識(shí)此類系統(tǒng)產(chǎn)生分岔及混沌現(xiàn)象的演化機(jī)理.

動(dòng)力學(xué)模型；Routh-Hurwitz判別法；Hopf分岔；Lyapunov指數(shù)譜

環(huán)境問題是近些年來很重要的問題，影響人類以及生物種群的存在和發(fā)展人與生物種群之間的關(guān)系問題也是被廣泛關(guān)注的，所以為了人與自然和諧相處，研究種群成為許多學(xué)者感興趣的課題之一[1].迄今為止，許多科學(xué)家做了大量調(diào)查并發(fā)表了許多論文[2-4].捕食食餌模型是很重要的模型之一，捕食模型主要對(duì)生態(tài)系統(tǒng)中的捕食者與食餌的關(guān)系進(jìn)行了分析研究，保持各生物種群的穩(wěn)定，以免使得物種被滅絕，從而使得生態(tài)系統(tǒng)達(dá)到持久穩(wěn)定狀態(tài)，保持生物種群的多樣性.所以，研究食餌捕食者之間的相互作用關(guān)系具有很重要的意義[5-6].

1 多種群捕食與被捕食模型建立

通過相關(guān)的文獻(xiàn)資料[7-10]，分析多種群捕食者與被捕食者的生態(tài)動(dòng)力學(xué)模型，根據(jù)各種群的食物鏈的反饋機(jī)制(如圖1所示)捕食者——兩種被捕食者，被捕食者被兩種捕食者捕食，存在外物補(bǔ)充，兩種捕食者存在種內(nèi)競爭和種間競爭，并且靠捕食被捕食者繁殖.建立反應(yīng)各生態(tài)要素的多參數(shù)模型,選取主要影響參數(shù)作為控制變量，對(duì)捕食者與被捕食者之間的機(jī)理進(jìn)行分析和討論.

圖1 營養(yǎng)關(guān)系圖

根據(jù)捕食者與被捕食者之間的關(guān)系形成一個(gè)簡單生態(tài)系統(tǒng)，通過利用多種群生態(tài)學(xué)原理[11]，建立捕食者與兩種被捕食者的動(dòng)力學(xué)模型的表達(dá)式為

其中:xi=xi(t),x1表示被捕食者隨時(shí)間變化的密度，x2,x3分別表示兩種捕食者隨時(shí)間變化的密度；a10為被捕食者內(nèi)稟增長率，a20和a30分別為兩種捕食者內(nèi)稟增長率；b1和b2共同決定被捕食者密度增長的半飽和參數(shù)；a12、a13分別為被兩種捕食者捕食的相對(duì)減少比例；a21、a31分別為兩種捕食者捕食后的轉(zhuǎn)換的增長比例；c1、c2分別為兩種捕食者種內(nèi)競爭減少的比例，d1、d2分別為兩種捕食者相互競爭減少的比例.

2 分析系統(tǒng)模型的運(yùn)動(dòng)穩(wěn)定性

根據(jù)方程的一些基本特征,考慮食物鏈模型中各元素的物理意義以及在實(shí)際情況中的相互作用和影響.考慮Lyapunov運(yùn)動(dòng)穩(wěn)定性的判定條件及其相關(guān)穩(wěn)定理論去判斷上述方程的穩(wěn)定性.

先求方程的平衡點(diǎn)，令上述方程組的左端為零，即

解上述方程組得平衡點(diǎn)Q(x1*,x2*,x3*).通過坐標(biāo)平移，可令

則方程可化為

則可得式(1)的Jacobi矩陣K為

對(duì)于平衡態(tài)Q1(x1,0,0)，此時(shí)僅有被捕食者存在，此時(shí)可得Jacobi矩陣為

對(duì)于平衡態(tài)Q1(x1,x2,0)，這時(shí)只存在被捕食者和其中一種捕食者，此時(shí)的Jacobi矩陣為

此時(shí)如果運(yùn)用特征值判斷系統(tǒng)穩(wěn)定性，由于不確定量太多，不容易進(jìn)行判斷，所以直接運(yùn)用Routh-Hurwitz判別法去判定系統(tǒng)的穩(wěn)定性.由Routh-Hurwitz判別法可得三階方程穩(wěn)定的條件為:若已知系統(tǒng)特征方程為a3l3+a2l2+a1l+a0=0，則當(dāng)所有系數(shù)均為正數(shù)時(shí)系統(tǒng)穩(wěn)定，且要求

對(duì)于平衡態(tài)Q1(x1,x2,0)可求得此時(shí)系統(tǒng)的特征方程式為

λ3+α2λ2+α1λ1+α0=0

其中:

根據(jù)Routh-Hurwitz判別法可知，不等式條件(R)成立時(shí)，系統(tǒng)將在平衡態(tài)點(diǎn)Q2(x1,x2,0)穩(wěn)定.

隨著參數(shù)取值的變化，平衡態(tài)Q2(x1,x2,0)出現(xiàn)跨臨界分叉和超臨界分岔現(xiàn)象.

3 數(shù)值模擬

對(duì)于平衡態(tài)Q3(x1,x2,x3)，此時(shí)被捕食者和兩種捕食者都存在，此時(shí)Jacobi矩陣和特征方程太過繁瑣，將采取將參數(shù)賦值的方法進(jìn)行分析說明，這里我們將a10作為控制變量對(duì)方程平衡態(tài)的穩(wěn)定性進(jìn)行研究，取參數(shù)a20=0.013 25，a30=0.018，a21=0.4，a12=0.5，a31=0.3，a13=0.4，b1=2，b2=3，c1=0.08，c2=0.02，d1=0.02，d2=0.01.此時(shí)各變量的初始值如下：

x(0)=0.8 ，y(0)=0.3，z(0)=0.2；

將各參數(shù)代入系統(tǒng)方程，并進(jìn)行數(shù)值模擬，可以發(fā)現(xiàn)：當(dāng)a10在0.01～0.264 6之間變化時(shí)，系統(tǒng)的相圖將由極限環(huán)逐步變化為點(diǎn)之后再呈現(xiàn)極限環(huán)狀態(tài)，此時(shí)主要出現(xiàn)兩種形態(tài).當(dāng)在0.264 7～0.264 8之間取值時(shí)系統(tǒng)開始出現(xiàn)不穩(wěn)定狀態(tài)，此時(shí)為系統(tǒng)的分界點(diǎn).當(dāng)a10取值大于0.264 8時(shí)，系統(tǒng)將呈現(xiàn)發(fā)散狀態(tài)，此時(shí)將不存在混沌吸引子，需要人為控制或改變其他控制變量，從而使系統(tǒng)再次達(dá)到平衡狀態(tài)[13].下面主要就三種情況進(jìn)行詳細(xì)分析.

以a10為控制變量計(jì)算系統(tǒng)的Lyapunov指數(shù)(如圖2)，由圖2可以看出系統(tǒng)的最大Lyapunov指數(shù)隨著a10的改變而改變，所以，系統(tǒng)的穩(wěn)定性也

隨著a10的改變而改變；當(dāng)a10<0.28時(shí)，系統(tǒng)的最大Lyapunov指數(shù)小于零，所以此時(shí)系統(tǒng)是穩(wěn)定的；當(dāng)a10≥0.28時(shí)，系統(tǒng)的最大Lyapunov指數(shù)大于零，此時(shí)系統(tǒng)處于不穩(wěn)定狀態(tài)，系統(tǒng)的Lyapunov指數(shù)大幅度增加，系統(tǒng)出現(xiàn)不穩(wěn)定狀態(tài).

圖2 系統(tǒng)關(guān)于a10變化時(shí)的Lyapunov指數(shù)圖

當(dāng)a10=0.22時(shí)，通過數(shù)值模擬可得系統(tǒng)的相圖和各因素的時(shí)間歷程圖見圖3.

圖3 a10=0.22時(shí)系統(tǒng)的相圖及時(shí)間歷程圖

圖3中(A)為系統(tǒng)的相圖，(B)為X的時(shí)間歷程圖，(C)為Y的時(shí)間歷程圖，(D)為Z的時(shí)間歷程圖；通過數(shù)值模擬結(jié)果，主要表明隨著時(shí)間的發(fā)展被捕食者與兩種捕食者隨時(shí)間發(fā)生周期性變化，由圖像可以看出，此時(shí)一種捕食者競爭不過另一種捕食者而滅絕，而沒有滅絕的捕食者與被捕食者呈現(xiàn)周期增長，隨時(shí)間增長系統(tǒng)會(huì)達(dá)一種平衡狀態(tài)，呈現(xiàn)出極限環(huán).此時(shí)系統(tǒng)是相對(duì)穩(wěn)定的.

當(dāng)a10=0.145 1時(shí)，通過數(shù)值模擬得到系統(tǒng)的相圖和各因素的時(shí)間歷程圖見圖4.

同樣由模擬結(jié)果可以表明:其中一種捕食者競爭不過另一種捕食者而滅絕，但隨著時(shí)間的變化，另一種捕食者與被捕食者最后達(dá)到了一種穩(wěn)定的狀態(tài)，種群的數(shù)量變化波動(dòng)很小，這種狀態(tài)是相比較穩(wěn)定的一種狀態(tài)，此時(shí)如果改變其他參數(shù)也可以使得滅絕的捕食者被滅絕的時(shí)間延長甚至使得三種群的數(shù)量最后也都達(dá)到一種穩(wěn)定狀態(tài).

圖4 a10=0.145 1時(shí)系統(tǒng)的相圖及時(shí)間歷程圖

當(dāng)a10>0.264 7時(shí)，此時(shí)得到的結(jié)果沒有意義，系統(tǒng)不再出現(xiàn)極限環(huán)或平衡態(tài)，此時(shí)可以改變其他參變量進(jìn)行人為干預(yù)，也可達(dá)到一種平衡態(tài)，且在較長時(shí)間內(nèi)各物種的變化有一定波動(dòng)，不會(huì)出現(xiàn)滅絕.不妨取a10=0.27，重新定義參數(shù)：a21=0.5，a31=0.575，c1=0.0182，a10=0.27，a20=0.02.其余參數(shù)不變.通過數(shù)值模擬得到系統(tǒng)重新定義參數(shù)之前的相圖和各因素的時(shí)間歷程圖分別如圖5.

圖5 a10=0.27時(shí)系統(tǒng)的相圖及時(shí)間歷程圖

通過圖5可以看到，系統(tǒng)是不穩(wěn)定的，此時(shí)被捕食者的在初始時(shí)刻的種群密度較小，兩種捕食者的種群密度也極不穩(wěn)定，但一定時(shí)間后，被捕食者的種群密度劇增，其中一種捕食者滅絕，另一種捕食者劇增，這種情況極其不穩(wěn)定，所以此時(shí)需要進(jìn)行人為干預(yù)才可以使得系統(tǒng)重新達(dá)到一種理想的穩(wěn)定狀態(tài).

通過人為干預(yù)重新改變參數(shù)后可得到如下的相圖和各因素的時(shí)間歷程圖，見圖6.

圖6 a10=0.27時(shí)系統(tǒng)的相圖及時(shí)間歷程圖

通過圖6可以看出，參變量發(fā)生變化后，此時(shí)系統(tǒng)重新達(dá)到了一種平衡態(tài)，在較長時(shí)間內(nèi)，被捕食者和捕食者均波動(dòng)生長，被捕食者種群密度的增大或減小，直接影響兩種捕食者種群密度的增大或減小，各物種不會(huì)出現(xiàn)劇增或滅絕狀態(tài)，這種情況是相對(duì)穩(wěn)定的.

從上面的分析中可以看出，系統(tǒng)的穩(wěn)定性隨a10的不斷變化而變化.當(dāng)a10<0.264 7時(shí)，系統(tǒng)是相對(duì)穩(wěn)定的；當(dāng)a10>0.264 7時(shí)，系統(tǒng)是不穩(wěn)定的，此時(shí)要達(dá)到穩(wěn)定狀態(tài)則需要人為干預(yù)，從而使系統(tǒng)再次達(dá)到一種平衡態(tài).所以，控制被捕食者的種群密度是維持物種種群平衡的重要因素.因此我們可以得到預(yù)防和控制某物種爆發(fā)或急劇減少的一種方法:及時(shí)檢測被捕食者的種群密度，當(dāng)被捕食者的增長率超過0.27之后，應(yīng)該及時(shí)采取人為干預(yù)，改變其他因素，從而改變其他參變量，使得各物種重新處于穩(wěn)定狀態(tài).

4 結(jié) 語

本文利用現(xiàn)代動(dòng)力學(xué)理論和系統(tǒng)運(yùn)動(dòng)穩(wěn)定性理論對(duì)一類捕食者與被捕食者模型的穩(wěn)定性作了深入的分析和討論.結(jié)果表明，對(duì)于此類多種群的生態(tài)動(dòng)力學(xué)系統(tǒng)，當(dāng)參數(shù)變化達(dá)到某些閾值時(shí)，系統(tǒng)會(huì)發(fā)生Hopf分岔現(xiàn)象.其他很多文章也對(duì)生態(tài)系統(tǒng)做了一些分析[14-16]，結(jié)合本文可得出系統(tǒng)對(duì)應(yīng)的各種穩(wěn)定狀態(tài)是有條件的，參數(shù)和實(shí)際發(fā)生條件是比較一致的.研究結(jié)果有助于更好的理解各種捕食者與被捕食者種群動(dòng)力學(xué)行為以及建立更有效的捕食與被捕食模型，具有一定的借鑒意義.但本文也有一些不足，缺乏一定的精確的數(shù)據(jù)來源.現(xiàn)實(shí)中在研究種群競爭時(shí)，可以在準(zhǔn)確數(shù)據(jù)來源的基礎(chǔ)上更加深入的研究，通過數(shù)學(xué)的方法使得種群維持在一種穩(wěn)定狀態(tài)，并及時(shí)給出理論上的分析與建議供相關(guān)機(jī)構(gòu)參考.

[1] CUI J, SUN Y. Permanence of predator-prey system with infinite delay [J]. Electronic Journal of Differential Equations I, 2004(81): 1-12.

[2] 黃運(yùn)金, 張朝陽. 無窮時(shí)滯Lotka-Volterra競爭系統(tǒng)的持久性[J]. 福建農(nóng)林大學(xué)學(xué)報(bào):自然科學(xué)版, 2010, 39(2): 222-224．

[3] 陳蘭蓀. 數(shù)學(xué)生態(tài)學(xué)模型與研究方法[M]. 北京: 科學(xué)出版社, 1988. 129-132．

[4] 陳蘭蓀, 陳 健.非線性生物動(dòng)力系統(tǒng)[M]. 北京: 科學(xué)出版社, 1993.165-170.

[5] 王洪禮, 馮劍豐, 沈 菲, 等. 赤潮藻類非線性動(dòng)力學(xué)模型的分岔及穩(wěn)定性研究[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2005, 26(6): 671-676.

[6] 馬知恩. 種群生態(tài)學(xué)的數(shù)學(xué)建模與研究[M]. 合肥: 安徽教育出版社, 1996. 16-159.

[7] 岳宗敏, 胡志興. 一類具功能反應(yīng)的食餌-捕食兩種群模型的極限環(huán)的唯一性[J]. 生物數(shù)學(xué)學(xué)報(bào), 2005, 20(2): 169-172.

[8] YANG W, LI X, BAI Z. Permanence of periodic Holling type-IV predator-prey system with stage structure for prey[J]. Mathematical and Computer Modelling, 2008, 48: 677-684.

[9] LU C, DING X, LIU M. The numerical simulation of periodic solutions for a predator-prey system[J]. Computers and Math with Applications, 2010, 59: 868-879.

[10] YAN X P, ZHANG C C. Hopf bifurcation in a delayed Lokta-Volterra predator-prey system[J]. Nonlinear Analysis Real World Applications, 2008, 9(1): 114-127.

[11] LV J, WANG K. Asymptotic properties of a stochastic predator-prey system with Holiing II functional response[J]. Commun Nonlinear SciNumer Somulat, 2011, 16: 4037-4048.

[12] WANG L L, LI W T. Periodic solutions and permanence for a delayed nonauton-omous ratio-dependent predator-prey model with Holling type functional response[J]. Comput Appl Math, 2004, 162: 341-357.

[13] FENG W. Dynamics in 3-species predator-prey models with time delays[J]. Discrete and Continuous Dynamical Systems, 2007, supplement: 364-372.

[14] TOAHA S, HASSAN M A. Stability analysis of predator-prey population model with time delay and constant rate of harvesting[J]. Journal of Mathematics, 2008, 40: 37-48.

[15] GE Z H, HE Y N. Diffusion effect and stability analysis of a predator-prey system described by a delayed reaction-diffusion equations [J]. J Math Anal Appl, 2008, 339: 1432-1450.

[16] CUI R H, SHI J P, WU B Y. StrogAllee effect in a diffusive predator-prey system with a protection zone[J]. J Diff Equns, 2014, 256: 108-129.

Stability analysis of dynamical model of two predators and one prey

SUN Dong-mei, HUANG Dong-wei

(School of Science, Tianjin Polytechnic University, Tianjin 300387, China)

In this paper, the interrelation between the density of two predators and one prey was investigated, and a mathematical model was established to describe dynamical behaviors between the three. The stability state of system was analyzed and judged by Routh-Hurwitz criterion. And chose the increase rate of the density of prey growth as the key parameter, could get the conditions for Hopf bifurcation occurrence in the system. And calculated the Lyapunov exponent spectrum, and presented simulation results of the system, and obtained the threshold parameter of Hopf bifurcation occurrence. The results showed that the density of prey played a vital role in the change of system state，which would be helpful to understand and identify the evolutionary mechanism of bifurcation and chaos phenomena of the system.

dynamical model; Routh-Hurwitz criterion; Hopf bifurcation; Lyapunov exponent spectrum

2016-05-25.

國家自然科學(xué)基金項(xiàng)目(11501410).

孫冬梅(1991-)，女，碩士，研究方向：非線性動(dòng)力系統(tǒng)理論及應(yīng)用.

黃東衛(wèi)(1966-)，男，博士，教授，碩士生導(dǎo)師，研究方向：非線性動(dòng)力系統(tǒng)理論及應(yīng)用.

O175

A

1672-0946(2016)06-0697-06