一類食餌-捕食者模型的二階龍格庫塔方法穩(wěn)定性及分支分析

柯于勝,陳伯山,劉唯一，2

(1.湖北師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 湖北 黃石 4350022.咸寧職業(yè)技術(shù)學(xué)院 工學(xué)院，湖北 咸寧 437100)

一類食餌-捕食者模型的二階龍格庫塔方法穩(wěn)定性及分支分析

柯于勝1,陳伯山1,劉唯一1，2

(1.湖北師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 湖北 黃石 4350022.咸寧職業(yè)技術(shù)學(xué)院 工學(xué)院，湖北 咸寧 437100)

研究了一類食餌-捕食者模型在R6內(nèi)的離散化及其動(dòng)力學(xué)行為.首先,利用二階龍格-庫塔方法將一類食餌-捕食者模型離散化,得到一類新的離散奇異系統(tǒng),然后運(yùn)用微分代數(shù)系統(tǒng)理論與分支理論討論了系統(tǒng)在平衡點(diǎn)處的局部穩(wěn)定性與分支問題,證明了Neimark-Sacker分支的存在性,并且選取捕獲努力研究了Neimark-Sacker分支及其方向,最后通過數(shù)值模擬證明了我們的結(jié)論.

食餌-捕食者模型;微分代數(shù)系統(tǒng);龍格-庫塔方法

在生態(tài)動(dòng)力系統(tǒng)中,傳統(tǒng)的Lotka-Volterra食餌-捕食者模型是一類十分重要的模型,在捕撈業(yè)以及生態(tài)管理中廣泛地應(yīng)用,具有重要的意義.在食餌-捕食者系統(tǒng)中,對(duì)于生命長(zhǎng)、世代重疊并且數(shù)量很大的種群,常近似地用微分方程來描述[1～2].生命短、世代不重疊或者數(shù)量少的種群,均常用差分方程來描述.近年來,這種用差分方程來描述的離散食餌-捕食者系統(tǒng)得到越來越多的關(guān)注.學(xué)者們研究了離散的Lotka-Volterra食餌-捕食者系統(tǒng)各方面的內(nèi)容,其中包括平衡點(diǎn)的穩(wěn)定性,分支情況,周期解的存在性和混沌控制等[3～7].

在差分動(dòng)力系統(tǒng)的研究中,許多用差分系統(tǒng)所描述的離散模型都是通過連續(xù)模型的離散化獲得的.然而,對(duì)于同樣的連續(xù)模型,不同的離散化方法可以得到不同的離散模型,進(jìn)而可以得到不同的定理與結(jié)論.大多數(shù)文獻(xiàn)中的離散化模型都是通過對(duì)連續(xù)模型進(jìn)行歐拉方法的離散化而得到的,并沒有獲得較為精確的離散化系統(tǒng).因此,本文希望在傳統(tǒng)的Lotka-Volterra食餌-捕食者系統(tǒng)的基礎(chǔ)上,用一種新的方法對(duì)連續(xù)模型進(jìn)行離散化,從而得到一類理想的離散模型,并分析該模型的動(dòng)力學(xué)行為.

1 模型建立

本文將考慮如下Lotka-Volterra食餌-捕食者模型:

(1)

其中x=x(t) 和y=y(t) 分別表示食餌種群和捕食者種群在時(shí)刻t的數(shù)量,k表示食餌種群的容納量,r是食餌種群的內(nèi)稟增長(zhǎng)率,b是捕食者的功能性反應(yīng),c是食餌向捕食者轉(zhuǎn)化的速率,d是捕食者種群的死亡率,e是捕獲努力.

首先,利用無量綱變換,將系統(tǒng) (1)化簡(jiǎn)為:

(2)

然后利用二階龍格-庫塔方法對(duì)系統(tǒng) (2)進(jìn)行離散化,得到下面的離散奇異系統(tǒng)

(3)

通過如上分析,建立了一類定義在R6中的二維流形(坐標(biāo)軸分別為x,y,f1,g1,f2,g2).其中,

f1(g1) 表示系統(tǒng) (2)的解x=x(t)(y=y(t)) 在區(qū)間段[ti,ti+1] 左端點(diǎn)的斜率,f2(g2) 表示系統(tǒng) (2)的解x=x(t)(y=y(t))在區(qū)間段[ti,ti+1] 右端點(diǎn)的斜率,δ表示步長(zhǎng).

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

容易看出系統(tǒng)(3) 有平衡點(diǎn)E0(x0,y0,f10,g10,f20,g20)其中

x0=(d+e)/c,y0=1-x0

f10=f20=g10=g20=0

現(xiàn)在考慮系統(tǒng)(3) 在平衡點(diǎn)E0(x0,y0,f10,g10,f20,g20) 的穩(wěn)定性.先求出系統(tǒng)(3) 在平衡點(diǎn)

E0(x0,y0,f10,g10,f20,g20) 的雅可比矩陣為

那么J(E0) 相應(yīng)的特征方程為

λ2+Pλ+Q=0

(4)

對(duì)于特征方程(4),F(λ)=λ2+Pλ+Q, 則

其中

B1=c,B2=2-c

根據(jù)系統(tǒng)(3)知道00 .

經(jīng)過分析,得到離散系統(tǒng) (3)在平衡點(diǎn)的局部穩(wěn)定性定理如下:

定理1 系統(tǒng)(3) 的平衡點(diǎn)E0(x0,y0,f10,g10,f20,g20) 為:

1)匯(源)當(dāng)且僅當(dāng)F(1)>0并且Q<1(Q>1 );

2)鞍點(diǎn)當(dāng)且僅當(dāng)F(1)<0;

3)非雙曲的當(dāng)且僅當(dāng)如下條件之一成立:

i)F(1)=0; ii)Q=1 .

本文著重研究系統(tǒng)(3) 的平衡點(diǎn)E0(x0,y0,f10,g10,f20,g20) 在參數(shù)經(jīng)歷一個(gè)閥值時(shí)發(fā)生Neimark-Sacker分支的動(dòng)力學(xué)形為.事實(shí)上,不妨選定系數(shù)c和d,取滿足條件

的適合參數(shù)值c,d,便可以找到相應(yīng)的e值,令其為e0.

3 平衡點(diǎn)E0 的Neimark-Sacker分支分析

選取參數(shù)(c,d,e)∈H,令e=e0+e*,并且選取e*為分支參數(shù),系統(tǒng)(3) 在微小擾動(dòng)下變?yōu)?/p>

(5)

其中 |e*|?1,為一個(gè)很小的擾動(dòng)參數(shù).

運(yùn)用中心流形定理和微分代數(shù)系統(tǒng)理論,我們得到系統(tǒng)(5) 的標(biāo)準(zhǔn)型:

(6)

系統(tǒng)(6)在(z1,z2)=(0,0) 處的線性化系統(tǒng)的特征方程為

λ2+P(e*)λ+Q(e*)=0

通過化簡(jiǎn)我們可以得到系統(tǒng)(6)在e*=0 的范式為

(7)

系統(tǒng)(7)中的G(u,v) 和K(u,v) 為

G(u,v)=a11u2+a12uv+a22v2+a111u3+a112u2v+a122uv2+a222v3+O((|u|+|v|)4)

K(u,v)=b11u2+b12uv+b22v2+b111u3+b112u2v+b122uv2+b222v3+O((|u|+|v|)4)

由G1(z1,z2) 和G2(z1,z2) 可以計(jì)算系統(tǒng)(6)在 (0,0)處的二階,三階導(dǎo)數(shù)分別為

Guu=2a11,Guv=a12,Gvv=2a22,Guuu1=6a111,Guuv=2a112,Guvv=2a122,Gvvv=6a222,

Kuu= 2b11,Kuv=b12,Kvv=2b22,Kuuu=6b111,Kuuv=2b112,Kuvv=2b122,Kvvv=6b222,

如果系統(tǒng)(7)發(fā)生Neimark-Sacker分支,必須要求下面的判定值不為零:

其中

經(jīng)過分析,有如下結(jié)論:

定理2 設(shè)(c,d,δ,e)∈H,若β≠0,則當(dāng)參數(shù)e在e0的小范圍內(nèi)變化時(shí),系統(tǒng)(3)的平衡點(diǎn)E0會(huì)經(jīng)歷一個(gè)Neimark-Sacker分支;若β>0 時(shí),則當(dāng)e>e0時(shí)平衡點(diǎn)E0不穩(wěn)定且有一個(gè)漸進(jìn)穩(wěn)定的周期軌道;當(dāng)ee0時(shí)平衡點(diǎn)E0是漸近穩(wěn)定的.

4 數(shù)值模擬

我們令c=12,d=8.5,δ=1,根據(jù)前面分析知道方程Q=1有且只有一個(gè)正根e0≈0.405,此時(shí)系統(tǒng)有一個(gè)平衡點(diǎn)E0(0.7421,0.2579,0,0,0)當(dāng)參數(shù)e在e0的范圍內(nèi)變化時(shí)會(huì)經(jīng)歷一個(gè)Neimark-Sacker分支,經(jīng)過計(jì)算得

因此,若e>e0,則平衡點(diǎn)E0漸近穩(wěn)定,若e

在圖1中,e=0.5,則平衡點(diǎn)E0(0.7421,0.2579,0.,0.,0,0) 漸近穩(wěn)定.

在圖2中,e=0.3,則平衡點(diǎn)E0(0.7421,0.2579,0,0,0,0) 不穩(wěn)定.

在圖3中,e=0.405,則在平衡點(diǎn)E0(0.7421,0.2579,0,0,0,0) 處產(chǎn)生一個(gè)漸近穩(wěn)定的周期軌道.

圖1 當(dāng)c=12,d=8.5,δ=1,e=0.5 時(shí)的Neimark-Sacker分支圖

圖2 當(dāng)c=12,d=8.5,δ=1,e=0.3 時(shí)的Neimark-Sacker分支圖

圖3 當(dāng)c=12,d=8.5,δ=1,e=0.405 時(shí)的Neimark-Sacker分支圖

[1]陳伯山, 劉永清. 非線性微分代數(shù)系統(tǒng)的穩(wěn)定性[J]. 控制理論與應(yīng)用, 2000, 17(1), 40～44.

[2]陳伯山, 廖曉昕. 微分代數(shù)系統(tǒng)的標(biāo)準(zhǔn)型與分支[J]. 應(yīng)用數(shù)學(xué)學(xué)報(bào), 2000, 23(3), 429～443.

[3]He Zhiming, Lai Xin. Bifurcation and chaotic behavior of a discrete-time predator-prey system[J]. Nonlinear Aalysis: Real World Application, 2011, 12: 403～417.

[4]Jing Zhujun, Yang Jianping. Bifurcation and chaos in discrete-time predator-prey system[J]. Chaos, Solitons and Fractals, 2006, 27: 259～277.

[5]Liu Xiaoli, Xiao Dongmei. Complex dynamic behaviors of a discrete-time predator-prey system[J]. Chaos, Solitons and Fractals, 2007, 32: 80～94.

[6]Hu Zengyun, Teng Zhidong, Zhang Long. Stability and Bifurcation analysis of a discrete predator-prey model with nonmonotonic functional response. Nonlinear Analysis: Real World Application, 2011, 12: 2356～2377.

[7]Zhang Guodong, Zhu Lulu, Chen Boshan. Hopf bifurcation and stability for a differential-algebraic biological economic system[J]. Applied Mathematics and Computation, 2010, 217: 330～338.

[8]劉唯一，傅朝金，陳 靜，等.一類具有非線性收獲率的捕食者-食餌生態(tài)經(jīng)濟(jì)系統(tǒng)的分支分析[J].湖北師范學(xué)院學(xué)報(bào)(自然科學(xué)版)，2014，34(2):86～89.

StabilityandbifurcationanalysisofasecondorderRunge-Kuttamethodforaprey-predatormodel

KE Yu-sheng1, CHEN Bo-shan1,LIU Wei-yi1,2

(1. College of Mathematics and Statistic,Hubei Normal University,Huangshi 435002,China;2.Institute of Technology,Xianning Vocational Technical College,Xianning 437100,China)

This paper investigate the discretization and the dynamical behavior of a prey-predator model in R6, First, applying the second order Runge-Kutta method, we obtain a new discrete system. Second, by using the differential-algebraic theory and bifurcation theory, the local stability and bifurcation of the proposed model around an interior equilibrium is discussed, and prove the Neimark-Sacker bifurcation exist. Then we analysis the Neimark-Sacker bifurcation and its direction by choosing the harvesting effort as the bifurcation parameter. Finally, a numerical simulation is presented to illustrate the theoretical results.

prey-predator model; local stability; Runge-Kutta method

2014—01—02

柯于勝(1987— )，男，湖北黃石人，碩士研究生，主要從事微分方程與控制論.

O193

A

1009-2714(2014)03- 0068- 06

10.3969/j.issn.1009-2714.2014.03.016