單種群動(dòng)力學(xué)模型的Matlab仿真

劉 波

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

種群動(dòng)力學(xué)模型[1]用于描述、預(yù)測(cè)以及調(diào)控種群的發(fā)展過(guò)程與發(fā)展趨勢(shì)，是幫助人們制訂和修正種群資源開(kāi)發(fā)策略的重要工具.單種群模型作為整個(gè)種群生態(tài)學(xué)的基石，由于其形式簡(jiǎn)單、參數(shù)容易確定、對(duì)大量實(shí)際數(shù)據(jù)都能完美地?cái)M合，因而受到生物學(xué)家和生物數(shù)學(xué)工作者的青睞[2].對(duì)單種群動(dòng)力學(xué)模型的模型建立和理論分析，已有較為廣泛的研究[2-4]，而作為一種模型檢驗(yàn)和預(yù)測(cè)的重要手段，計(jì)算機(jī)數(shù)值仿真往往沒(méi)有受到應(yīng)有的重視.計(jì)算機(jī)數(shù)值仿真可以將抽象復(fù)雜的理論研究結(jié)果，轉(zhuǎn)化成易于人們理解的數(shù)值序列或圖形演示，使復(fù)雜的動(dòng)力學(xué)分析和數(shù)量巨大的參數(shù)篩選及動(dòng)態(tài)過(guò)程的模擬成為可能[5].基于Matlab的單種群離散模型數(shù)值仿真手段主要有三種：線性聯(lián)結(jié)圖、蛛網(wǎng)圖、費(fèi)根鮑姆圖[6]，目前使用最為廣泛的是線性聯(lián)結(jié)圖[7]，而蛛網(wǎng)圖和費(fèi)根鮑姆圖則較少應(yīng)用.線性聯(lián)結(jié)圖是種群數(shù)量或密度隨著時(shí)間變化的序列圖像，能夠直觀反映種群隨著時(shí)間演化的持續(xù)生存性、絕滅性和周期性.線性聯(lián)結(jié)圖繪制簡(jiǎn)單易于理解，但在確定系統(tǒng)的周期性上存在缺陷，一般只能確定近似周期，如果能將蛛網(wǎng)圖和費(fèi)根鮑姆圖結(jié)合起來(lái)，不但能比較準(zhǔn)確地確定系統(tǒng)的周期性，而且還可以進(jìn)一步揭示系統(tǒng)更為復(fù)雜的漸近行為.

1 單種群離散模型的Matlab仿真

考慮一類簡(jiǎn)單的Logistic單種群離散模型：ut+1=rut(1-ut),r>0.

ut和ut+1分別是種群在t時(shí)刻和t+1時(shí)刻的種群密度，為保證解的正性，通常假設(shè)初值滿足0

1.1 線性聯(lián)結(jié)圖

對(duì)迭代函數(shù)xn+1=f(xn),n=0,1,2,…，在以迭代次數(shù)n為橫軸、迭代值xn為縱軸的平面直角坐標(biāo)系中，不斷以直線段聯(lián)結(jié)相鄰迭代點(diǎn)列形成的折線圖，稱為線性聯(lián)結(jié)圖.線性聯(lián)結(jié)圖可以很好地揭示迭代序列在無(wú)窮遠(yuǎn)處的漸近行為，通常用來(lái)觀察種群動(dòng)力學(xué)模型的持續(xù)生存性、絕滅性和近似周期性.

下面給出繪制系統(tǒng)(1)線性聯(lián)結(jié)圖的Matlab程序及圖例：

function linlink(x,a,nspan) %linlink(0.79,3.8284,[451,1000]);

funstr=strrep('a*x*(1-x)','a',num2str(a)); cmdstr=['x(n)=',strrep(funstr,'x','x(n-1)'),';'];

n_start=nspan(1); n_final=nspan(end);

for n=2:n_final, eval(cmdstr); end;

figure; plot(n_start:n_final,x(n_start:n_final)); xlabel(['n from ',num2str(n_start),' to ',num2str(n_final)]);

ylabel([cmdstr,' with x(0)=',num2str(x(1))]); title('Linear Connection Diagram'); zoom xon;

(a)linlink(0.20,3.5000,[800,900]);(b) linlink(0.16,3.5000,[120,300]);(c)linlink(0.79,3.8284,[640,820])圖1 線性聯(lián)結(jié)圖Fig.1 Linear conneition graph

通過(guò)觀察圖1(a)～(c)所示線性聯(lián)結(jié)圖，可以分析得知種群在不同初值條件和參數(shù)條件下，系統(tǒng)的持續(xù)生存性容易確定，周期性尚不明確，例如圖2(a)振蕩發(fā)散，存在近似周期；圖2(b)和圖2(c)振蕩發(fā)散，周期性不易確定.

1.2 蛛網(wǎng)圖

對(duì)迭代函數(shù)，在平面直角坐標(biāo)系畫出直線和，自上某初始點(diǎn)引垂直線段與相交，再以交點(diǎn)引水平線段與相交，如此反復(fù)形成的折線圖形似階梯或蛛網(wǎng)，稱為蛛網(wǎng)圖.蛛網(wǎng)圖可以很好地揭示迭代序列在無(wú)窮遠(yuǎn)處的漸近行為，通常用來(lái)觀察種群動(dòng)力學(xué)模型的持續(xù)生存性、絕滅性和準(zhǔn)確周期性.

下面給出繪制系統(tǒng)(1)蛛網(wǎng)圖的Matlab程序及圖例：

function spilink(x,a,nspan) %spilink(0.2,3.5644,[1,500])

funstr=strrep('a*x*(1-x)','a',num2str(a)); cmdstr=['x(n)=',strrep(funstr,'x','x(n-1)'),';'];

n_start=nspan(1); n_final=nspan(end);

for n=2:n_final+1, eval(cmdstr); end;

figure; ezplot('x',[-0.2,1.2]); hold on; ezplot(funstr,[-0.2,1.2]);

ylabel([cmdstr,' with x(0)=',num2str(x(1))]); title('Spider Connection Diagram')

for n=n_start:n_final

plot([x(n),x(n)],[x(n),x(n+1)],'r-'); plot([x(n),x(n+1)],[x(n+1),x(n+1)],'r-');

xlabel(['n of x(n) is arriving at ',num2str(n)]); drawnow;

end

(a)spilink(0.1,2.5300,[1,100]);(b)spilink(0.2,3.5644,[900,1200]);(c)spilink(0.1,3.1691,[1300,1400])圖2 蛛網(wǎng)圖Fig.2 Spider graph

通過(guò)觀察圖2(a)～(c)所示蛛網(wǎng)圖，可以分析得知種群在不同初值條件和參數(shù)條件下，系統(tǒng)的持續(xù)生存性、周期性容易明確，例如圖2(a)收斂，存在穩(wěn)定點(diǎn)；圖2(b)振蕩發(fā)散，存在近似8-周期；圖2(c)振蕩發(fā)散，存在準(zhǔn)確2-周期.

1.3 費(fèi)根鮑姆圖

對(duì)迭代函數(shù)，在以參數(shù)為橫軸、迭代序列為縱軸的平面直角坐標(biāo)系中，對(duì)應(yīng)不同的，在過(guò)點(diǎn)垂直橫軸的方向，繪制所有點(diǎn)，如此反復(fù)形成的散點(diǎn)圖形，稱為費(fèi)根鮑姆圖.費(fèi)根鮑姆圖可以很好地揭示迭代序列在不同參數(shù)條件下的漸近行為，通常用來(lái)觀察種群動(dòng)力學(xué)模型的周期性、倍周期分岔和混沌等復(fù)雜現(xiàn)象.

下面給出繪制系統(tǒng)(1)費(fèi)根鮑姆圖的Matlab程序及圖例：

function fgblink(x,aspan,nspan) %fgblink(0.2,[2,4],[900,1300])

warning off; a_start=aspan(1); a_final=aspan(end); n_start=nspan(1); n_final=nspan(end); figure;

for a=a_start:0.01:a_final

funstr=strrep('a*x*(1-x)','a',num2str(a));cmdstr=['x(n)=',strrep(funstr,'x','x(n-1)'),';'];

for n=2:n_final+1, eval(cmdstr); end;

plot(a*ones(1,max(nspan)-min(nspan)+1),x(n_start:n_final),'r+'); drawnow; hold on;

title(['a is arriving at ',num2str(a)]);

end

xlabel('a'); ylabel(['a*x*(1-x) with x(0)=',num2str(x(1))]); zoom xon;

(a) fgblink(0.2,[0,4],[1,500]);(b) fgblink(0.2,[2,4],[900,1300]);(c)fgblink(0.2,[3.70,3.95],[10000,12000])圖3 費(fèi)極鮑姆圖Fig.3 Feigenbanm graph

通過(guò)觀察圖3(a)～(c)所示費(fèi)根鮑姆圖，可以分析得知種群在不同參數(shù)條件下，系統(tǒng)出現(xiàn)倍周期分岔和混沌等復(fù)雜現(xiàn)象，例如圖3(a)在取1.25時(shí)，存在穩(wěn)定點(diǎn)；圖3(b)在間于2.8和3.7時(shí)，出現(xiàn)倍周期分岔和混沌；圖3(c)在間于3.7和3.9時(shí)，出現(xiàn)準(zhǔn)確3-周期或5-周期.

2 總結(jié)和討論

本文分析和對(duì)比了三種不同的數(shù)值仿真手段在種群動(dòng)力學(xué)模型仿真中的應(yīng)用，通過(guò)編寫Matlab程序?qū)σ活惤?jīng)典的Logistic單種群離散模型進(jìn)行數(shù)值仿真，結(jié)果表明目前使用最為廣泛的線性聯(lián)結(jié)圖只能用于確定近似周期，蛛網(wǎng)圖則可以很好地彌補(bǔ)這一缺陷，而費(fèi)根鮑姆圖在探究動(dòng)力學(xué)模型更為復(fù)雜的漸近行為時(shí)，具有獨(dú)特的優(yōu)勢(shì).

眾多學(xué)者指出對(duì)于種群數(shù)目小，世代不重疊種群來(lái)說(shuō)，差分系統(tǒng)表示的離散模型比連續(xù)模型更為合適[8]，而且離散模型可以提供比連續(xù)模型更為有效的計(jì)算機(jī)數(shù)值仿真.因此，某些連續(xù)模型在一定條件下也可以通過(guò)差分化轉(zhuǎn)化成相應(yīng)的離散模型，綜合應(yīng)用本文介紹的三種數(shù)值仿真手段，可以幫助更加深刻地理解數(shù)學(xué)模型背后的生態(tài)學(xué)意義.

[1] 陳蘭蓀,王東達(dá).數(shù)學(xué)、物理學(xué)與生態(tài)學(xué)的結(jié)合─一種群動(dòng)力學(xué)模型[J].物理,1994,23(7):408-413.

[2] 唐三一,肖燕妮.單種群生物動(dòng)力系統(tǒng)[M].北京:科學(xué)出版社,2008.

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

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

[5] 黃樨.數(shù)學(xué)與計(jì)算機(jī)仿真在生態(tài)學(xué)研究中的應(yīng)用[J].南京林業(yè)大學(xué)學(xué)報(bào),2009,25(5):63-66.

[6] 桂占吉.生物動(dòng)力學(xué)模型與計(jì)算機(jī)仿真[M].北京:科學(xué)出版社,2005.

[7] 杜斌,張彥民.MATLAB在種群動(dòng)力學(xué)模型中的仿真研究[J].商洛學(xué)院學(xué)報(bào),2010,24(2):28-31.

[8] Murry J D. Mathematical biology[M].Berlin:Springer-Verlag,1989.