劉 波

(湖北民族學(xué)院理學(xué)院，湖北恩施445000)

對種群動力學(xué)模型而言，計算機數(shù)值仿真是模型檢驗和預(yù)測的一種重要手段［1－3］.計算機數(shù)值仿真可以將抽象復(fù)雜的理論研究結(jié)果，轉(zhuǎn)化成易于人們理解的數(shù)值序列或圖形演示，使復(fù)雜的動力學(xué)分析和數(shù)量巨大的參數(shù)篩選及動態(tài)過程的模擬成為可能［4］.基于Matlab的兩種群動力學(xué)模型數(shù)值仿真手段主要有三種:時間序列圖、斜率場圖、相平面軌跡圖［5－7］，時間序列圖能夠直觀反映種群隨著時間演化的趨勢，例如持續(xù)生存性、絕滅性和周期性.斜率場圖和相平面軌跡圖都屬于相平面圖形，尤其適合反映種群的相關(guān)變化趨勢.斜率場圖的最大優(yōu)勢在于從整體反映斜率場的分布，且不必求解微分方程.相平面軌跡圖的最大優(yōu)勢在于對應(yīng)不同初值條件反映軌線的走向，有助于分析系統(tǒng)的穩(wěn)定性.綜合運用以上三種數(shù)值仿真手段，可以進一步揭示兩種群動力學(xué)模型更為復(fù)雜的漸近行為.

1 兩種群動力學(xué)模型的Matlab仿真

考慮一類經(jīng)典的Lotka－volterra兩種群動力學(xué)模型:

r1和r2分別表示兩種群的內(nèi)稟增長率，a和d分別表示兩種群的種內(nèi)作用系數(shù)，b和c分別表示兩種群的種間作用系數(shù)，為了揭示模型(1)的漸近行為，下面分別采用時間序列圖、斜率場圖、相平面軌跡圖三種不同的仿真手段，編寫Matlab程序?qū)δＰ?1)進行數(shù)值仿真.根據(jù)實際生物意義，以下僅在區(qū)域R2+={(x，y)|x≥0，y≥0}中繪制圖形.

1.1 時間序列圖

在模型(1)中，取 r1=1，a=0，b=－0.1，r2=－0.5，c=0.02，d=0 得到如下系統(tǒng):

對系統(tǒng)(2)，只要給定初值條件(x0，y0)，即可利用函數(shù)ode45求得數(shù)值解，由此繪制的圖形x(t)和y(t)稱為時間序列圖.下面給出繪制系統(tǒng)(2)時間序列圖的Matlab程序及圖例:

function xdot=TimeSeries(t，x)

r1=1;a=0;b=－0.1;r2=－0.5;c=0.02;d=0;

xdot=［x(1)*(r1+a*x(1)+b*x(2));x(2)*(r2+c*x(1)+d*x(2))］;

functionLoadTimeSeries(x0，y0)%LoadTimeSeries(25，2)

tspan=0:0.01:25;cond0=［x0，y0］;

［t，x］=ode45(＇timeseries＇，tspan，cond0);

plot(t，x(:，1)，＇r－＇，t，x(:，2)，＇b－＇);legend(＇x(t)＇，＇y(t)＇);xlabel(＇t＇);ylabel(＇x(t)and y(t)＇);

通過觀察圖1所示時間序列圖，可以分析得知系統(tǒng)(2)在初值條件cond0=［25，2］下，種群x和種群y持續(xù)生存，且呈現(xiàn)出近似周期性.

1.2 斜率場圖

在系統(tǒng)(2)中，消去dt可以得到形如dy/dx的斜率表達式，在R2+中取定一矩形區(qū)域并稱之為場，自場中每一點(xi，yi)引出斜率為Ki的短線段，由此繪制的圖形就是斜率場圖.下面給出繪制系統(tǒng)(2)斜率場圖的Matlab程序及圖例:

function SlopesField(r1，a，b，r2，c，d)%SlopesField(1，0，－0.1，－0.5，0.02，0)

syms x y;fun=y*(r2+c*x+d*y)/(x*(r1+a*x+b*y));hx=120/30;hy=30/30;x0=0;y0=0;hold on

for i=1:30，x=x0+(i－1)*hx;for j=1:30，y=y0+(j－1)*hy;k=eval(fun);

if abs(hx*k)＞hy，plot(［x，x+hy/k*2/3］，［y，y+hy*2/3］);else plot(［x，x+hx*2/3］，［y，y+hx*k*2/3］);end;end;end;

xlabel(＇x＇);ylabel(＇y＇);box on;warning off;

通過觀察圖2所示斜率場圖，可以分析得知系統(tǒng)(2)在不同初值條件下，斜率場圍繞在正平衡點(25，10)周圍，種群x和種群y持續(xù)生存，從場中一點出發(fā)即可判斷出大致走向.

1.3 相平面軌跡圖

在系統(tǒng)(2)中，考慮將數(shù)值解法求得的時間序列數(shù)據(jù)x(t)和y(t)繪制在相平面XOY上，由此得到的圖形稱為相平面軌跡圖.下面給出繪制系統(tǒng)(2)相平面軌跡圖的Matlab程序及圖例:

通過觀察圖3所示相平面軌跡圖，可以分析得知系統(tǒng)(2)在不同初值條件，軌線以正平衡點(25，10)為中心環(huán)狀分布，種群x和種群y持續(xù)生存，進一步驗證可知系統(tǒng)是全局一致漸近穩(wěn)定的［7］.

圖1 時間序列圖Fig.1 Time series plot

圖2 斜率場圖Fig 2 Field of slopes

圖3 相平面軌跡圖Fig 3 Phase plane plot

2 總結(jié)和討論

本文分析和對比了三種不同的數(shù)值仿真手段在兩種群動力學(xué)模型仿真中的應(yīng)用，通過編寫Matlab程序?qū)σ活惤?jīng)典的Lotka－volterra兩種群動力學(xué)模型進行數(shù)值仿真，結(jié)果表明目前使用較為廣泛的時間序列圖在確定周期性上并不準確，斜率場圖最值得關(guān)注的特點是不必求解微分方程，相平面軌跡圖則彌補了斜率場圖在方向性上的不足，并且對系統(tǒng)穩(wěn)定性分析有直接幫助.綜合應(yīng)用本文介紹的三種數(shù)值仿真手段，可以幫助我們更加深刻地認識和理解種群動力學(xué)模型背后復(fù)雜的漸近行為.

［1］ 劉波.單種群動力學(xué)模型的MATLAB仿真［J］.湖北民族學(xué)院學(xué)報:自然科學(xué)版，2011，29(4):368－370.

［2］ 梁仁君，林振山，韓洪凌，等.物種對資源競爭的動力機制及數(shù)值模擬［J］.生態(tài)學(xué)報，2007，27(12):5390－5397.

［3］ 劉波，劉志軍.一類食餌具有S－型變收獲率的非自治離散捕食系統(tǒng)漸近行為［J］.生物數(shù)學(xué)學(xué)報，2012，27(3):499－506.

［4］ 黃樨.數(shù)學(xué)與計算機仿真在生態(tài)學(xué)研究中的應(yīng)用［J］.南京林業(yè)大學(xué)學(xué)報，2009，25(5):63－66.

［5］ 桂占吉.生物動力學(xué)模型與計算機仿真［M］.北京:科學(xué)出版社，2005.

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

［7］ 馬知恩，周義倉.常微分方程定性與穩(wěn)定性方法［M］.北京:科學(xué)出版社，2001.