西農(nóng) 建模實(shí)驗(yàn)六

1、實(shí)驗(yàn)六 機(jī)理模型與平衡原理實(shí)驗(yàn)?zāi)康?如果對(duì)所研究的問(wèn)題了解的比較深入，知道產(chǎn)生現(xiàn)象的內(nèi)在的機(jī)理，那么依據(jù)機(jī)理建模，則模型具有更好的可靠性和廣泛性。不考慮隨機(jī)因素，假設(shè)每一時(shí)刻是確定的如果對(duì)系統(tǒng)狀態(tài)的觀測(cè)和描述只在離散的時(shí)間點(diǎn)上，則構(gòu)成差分方程模型；如果考慮系統(tǒng)隨時(shí)間連續(xù)變化，則是微分方程模型。本節(jié)主要以這兩類(lèi)方程為例，介紹用MATLAB軟件求解機(jī)理模型的基本方法。差分方程模型一、實(shí)驗(yàn)題目 由一對(duì)兔子開(kāi)始，一年可以繁殖出多少只兔子？如果一對(duì)兔子每個(gè)月可以生一對(duì)小兔子，兔子在出生兩個(gè)月后就具有繁殖能力，由一對(duì)剛出生一個(gè)月的兔子開(kāi)始，一年內(nèi)兔子種群數(shù)量如何變化。求這個(gè)種群的穩(wěn)定分布和固有增長(zhǎng)率。二、

2、實(shí)驗(yàn)內(nèi)容解 假設(shè)（a） 兔子每經(jīng)過(guò)一個(gè)月底就增加一個(gè)月齡；（b） 月齡大于等于2的兔子都具有繁殖能力；（c） 具有繁殖能力的兔子每一個(gè)月一定生一對(duì)兔子；（d） 兔子不離開(kāi)群體（不考慮死亡） 記第n個(gè)月初的幼兔（一月齡兔）數(shù)量為a0(n）,成兔（月齡大于等于2）數(shù)量為a1(n），則兔子總數(shù)為a(n)= a0(n）+a1(n）,平衡關(guān)系為： 建立模型：這個(gè)一階差分方程的矩陣表達(dá)式為其中， 利用迭代方法求數(shù)值解，也就是按時(shí)間步長(zhǎng)法仿真種群增長(zhǎng)的動(dòng)態(tài)過(guò)程，模擬幼兔和成兔占整體比例隨時(shí)間的變化。a=01;11;x=10;fork=2:12y=a*x(:,k-1);x=xy;endzz=repmat(su

3、m(x),21);z=x./zz;t=1:12;plot(t,x(1,:),r,t,x(2,:),b),grid;plot(t,z(1,:),r,t,z(2,:),b),grid;由數(shù)值模擬結(jié)果可見(jiàn)，兔子數(shù)量遞增，但是幼兔和成兔在種群中所占比例很快會(huì)趨于一個(gè)極限。由線性差分方程的定性理論可知，這個(gè)極限就是對(duì)應(yīng)于差分方程系數(shù)矩陣A主特征值得歸一化特征向量。因?yàn)锳是非負(fù)矩陣，由矩陣?yán)碚撝?，A主特征值是正實(shí)數(shù)，是最大的特征值。v,d=eig(a)v=-0.85070.52570.52570.8507d=-0.6180001.6180max(max(d)ans=1.6180v(:,2)/sum(v(:

4、,2)ans=0.38200.6180由數(shù)值計(jì)算可知，系數(shù)矩陣模A最大的特征值r=1.618，生物上稱(chēng)之為種群的內(nèi)稟增長(zhǎng)率，是個(gè)大于一的實(shí)數(shù)。因此種群數(shù)量隨時(shí)間遞增。相應(yīng)的歸一化的特征向量的兩個(gè)分量0.382和0.618正是幼兔和成兔在種群中所占比例趨近的穩(wěn)定值，生物上稱(chēng)之為種群的穩(wěn)定分布。 從這個(gè)例子的討論可見(jiàn)，數(shù)值模擬能幫我們對(duì)系統(tǒng)的變化有更直觀的認(rèn)識(shí)。但是只有通過(guò)計(jì)算方程組系數(shù)矩陣的特征值和特征向量，運(yùn)用差分方程定性分析才能對(duì)解得漸進(jìn)性質(zhì)給出確定的結(jié)論。微分方程模型一、實(shí)驗(yàn)題目 藍(lán)鯨的內(nèi)稟增長(zhǎng)率每年估計(jì)為5%，估計(jì)藍(lán)鯨的最大環(huán)境載量為150000條，磷蝦是藍(lán)鯨喜歡的一種食物。磷蝦的最大飽

5、和種群為500噸/英畝*。當(dāng)缺少捕食者，環(huán)境不擁擠時(shí)，磷蝦種群以每年25%的速率增長(zhǎng)。磷蝦500噸/英畝可以提高藍(lán)鯨2%的年增長(zhǎng)率，同時(shí)150000條藍(lán)鯨將減少磷蝦10%的年增長(zhǎng)率。（1） 組建一個(gè)藍(lán)鯨和磷蝦的動(dòng)態(tài)模型，模擬兩個(gè)種群隨時(shí)間的變化。假設(shè)初始狀態(tài)為藍(lán)鯨5000條，磷蝦750噸/英畝；（2） 確定藍(lán)鯨與磷蝦是否可以長(zhǎng)期共存；（3） 假設(shè)捕撈使得藍(lán)鯨只剩下它的平衡態(tài)的5%，而磷蝦保持平衡態(tài)的數(shù)量。描述一旦停止捕撈將發(fā)生什么情況。藍(lán)鯨恢復(fù)需要多長(zhǎng)時(shí)間？磷蝦群體將發(fā)生什么變化？進(jìn)一步，給出藍(lán)鯨種群恢復(fù)時(shí)間對(duì)它所受傷害程度的依賴(lài)關(guān)系。二、實(shí)驗(yàn)內(nèi)容 解 ：（1）假設(shè)（a） 藍(lán)鯨和磷蝦單獨(dú)存在時(shí)

6、，兩種群都依靠有限的自然資源增長(zhǎng)，即遵循Logistic模型；（b） 藍(lán)鯨捕食磷蝦，藍(lán)鯨平均增長(zhǎng)率的增加量正比于單位區(qū)域內(nèi)磷蝦數(shù)量，磷蝦平均增長(zhǎng)率的減少量正比于藍(lán)鯨數(shù)量 記N1(t)為藍(lán)鯨在t時(shí)刻的種群數(shù)量（條），N2(t)為磷蝦在t時(shí)刻的種群質(zhì)量（噸/英畝），于是，依據(jù)假設(shè)，建立藍(lán)鯨和磷蝦兩個(gè)種群的動(dòng)態(tài)模型 其中r1=0.05和r2=0.25分別表示藍(lán)鯨和磷蝦種群的內(nèi)稟增長(zhǎng)率，k1=150000(條）和k2=500（噸/英畝）分別表示藍(lán)鯨和磷蝦種群環(huán)境載量，a12=0.02/500表示每英畝每噸磷蝦對(duì)藍(lán)鯨平均增長(zhǎng)率的改變，a21=0.10/150000表示每條藍(lán)鯨對(duì)磷蝦平均增長(zhǎng)率的改變。 下

7、一步，為了便于數(shù)值模擬，保持?jǐn)?shù)量級(jí)的協(xié)調(diào)，把鯨魚(yú)的單位改為以千條為單位。當(dāng)初始狀態(tài)為藍(lán)鯨5千條、磷蝦750噸/英畝時(shí)，動(dòng)態(tài)模型的數(shù)值模擬MATLAB指令為vG=(t,y)0.05*y(1)*(1-y(1)/150)+(0.02/500)*y(2)*y(1);0.25*y(2)*(1-y(2)/500)-(0.1/150)*y(1)*y(2); t,y=ode45(vG,0,150,5,750);plot(t,y(:,1),-.),grid on,hold on plot(t,y(:,2),-),grid on,xlabel(時(shí)間/年),ylabel(種群數(shù)量/千條); gtext(虛線表示藍(lán)鯨

8、，實(shí)線表示磷蝦);由圖可見(jiàn)，藍(lán)鯨的數(shù)量隨著時(shí)間的增加而逐漸增加，磷蝦的數(shù)量隨時(shí)間的增加而逐漸減少，最后都分別趨近一個(gè)穩(wěn)定值。（2）由常微分方程理論知，方程組的常數(shù)值解為方程的平衡點(diǎn)，由平衡點(diǎn)的穩(wěn)定性確定了方程解的動(dòng)態(tài)性質(zhì)，即解的漸近行為。上一問(wèn)的數(shù)值結(jié)果表明，方程組具有一個(gè)穩(wěn)定的正平衡點(diǎn)，也就是說(shuō)藍(lán)鯨與磷蝦可以長(zhǎng)期共存。首先求方程組的平衡點(diǎn)，Matlab指令為： syms y1 y2 f1=0.05*y1*(1-y1/150)+(0.02/500)*y2*y1; f2=0.25*y2*(1-y2/500)-(0.1/150)*y1*y2; y1steady,y2steady=solve(f1,

9、f2); disp(平衡點(diǎn)是：);平衡點(diǎn)是： disp(vpa(y1steady,6) vpa(y2steady,6); 0, 0 150.0, 0 0, 500.0 181.034, 258.621接著，考慮方程組在平衡點(diǎn)（N1,N2）=(xi,yi),i=1,2,3,4附近的局部線性方程零解的穩(wěn)定性這些線性方程組零解的穩(wěn)定性由其系數(shù)矩陣的特征值確定。利用Matlab指令求系數(shù)矩陣的特征值x=0,150,0,181.034;y=0,0,500,258.621;fori=1:4A=0.05-0.1*x(i)/150+0.02*y(i)/500,0.02*x(i)/500;-0.1*y(i)/1

10、500.25-0.5*y(i)/500-0.1*x(i)/150,;b=eig(A);disp(b(1)b(2)end0.05000.2500-0.05000.1500-0.25000.0700-0.0948+0.0077i-0.0948-0.0077i得到的結(jié)果表明，在正平衡點(diǎn)(181.034，258.621)，相應(yīng)的兩個(gè)特征值的實(shí)部都是負(fù)的。按照微分方程定性理論可知，方程組正平衡點(diǎn)(181.034，258.621)是漸近穩(wěn)定的，即從任意初值點(diǎn)出發(fā)，解軌線都會(huì)趨近該點(diǎn)。因此，可以斷言，只要停止捕撈，藍(lán)鯨與磷蝦可以長(zhǎng)期共存。(3) 為了確定當(dāng)藍(lán)鯨數(shù)量為平衡態(tài)數(shù)量的r%，磷蝦數(shù)量為平衡態(tài)時(shí)，停止

11、捕撈，藍(lán)鯨恢復(fù)到平衡態(tài)需要的時(shí)間，只有通過(guò)數(shù)值模擬。對(duì)不同的初值N1(0)=r/100181.034,N2(0)=258.621,在一個(gè)充分長(zhǎng)的時(shí)間區(qū)間上求方程組的數(shù)值解，注意到藍(lán)鯨數(shù)量會(huì)遞增趨于平衡態(tài)，可以N1(T)181.034-0.001確定的時(shí)間T近似表示種群恢復(fù)所需的時(shí)間，得到對(duì)應(yīng)的函數(shù)關(guān)系T=T(r). Matlab指令如下：G=(t,N)0.05*N(1)*(1-N(1)/150)+(0.02/500)*N(2)*N(1);0.25*N(2)*(1-N(2)/500)-(0.1/150)*N(1)*N(2);T=;forr=0.05:0.05:0.9t,N=ode45(G,0,2

12、00,181.034*r,258.621);subplot(1,3,1),plot(t,N(:,1),-,t,N(:,2),-.),xlabel(時(shí)間),ylabel(種群數(shù)量),gridon,holdond=find(N(:,1)181.034-0.001,1);T=Tt(d);endsubplot(1,3,2),plot(0.05:0.05:0.9,T),xlabel(損傷程度r),ylabel(恢復(fù)時(shí)間T),gridgtext(實(shí)線表示藍(lán)鯨，虛線表示磷蝦)由圖可知，得到的函數(shù)關(guān)系T=T(r)是不光滑的。注意到計(jì)算機(jī)只會(huì)計(jì)算離散值，盡管ode45采用了四階、五階龍格-庫(kù)塔法，可能是離散的時(shí)

13、間步長(zhǎng)太大，計(jì)算得結(jié)果并未反映連續(xù)系統(tǒng)的真實(shí)規(guī)律。重新規(guī)定步長(zhǎng)，計(jì)算量大些，但能夠保證離散結(jié)果逼近連續(xù)情形的精度。tspan=linspace(0,150,1000);T=;forr=0.05:0.05:0.9t,N=ode45(G,tspan,181.034*r,258.621);d=find(N(:,1)181.034-0.001,1);T=Tt(d);endsubplot(1,3,3),plot(0.05:0.05:0.9,T),xlabel(損傷程度r),ylabel(恢復(fù)時(shí)間T),grid果然得到理想的效果。在利用計(jì)算機(jī)進(jìn)行模型的數(shù)值求解時(shí)，對(duì)模型解的性質(zhì)先有一個(gè)定性分析是非常重要的

14、，以確保離散網(wǎng)格點(diǎn)上的解確實(shí)保持原連續(xù)系統(tǒng)的性質(zhì)。離散與連續(xù)一、實(shí)驗(yàn)題目 上節(jié)的例子說(shuō)明連續(xù)系統(tǒng)仿真必須限制離散步長(zhǎng)，否則可能會(huì)產(chǎn)生不同的結(jié)果，甚至導(dǎo)致復(fù)雜的動(dòng)力學(xué)行為，這些行為并不是連續(xù)系統(tǒng)所具有的。我們?cè)诮柚鷶?shù)值模擬研究連續(xù)系統(tǒng)時(shí)，計(jì)算機(jī)只會(huì)在離散點(diǎn)上計(jì)算，不同的離散格式同樣會(huì)導(dǎo)致不同的動(dòng)力學(xué)行為，有些行為并不是連續(xù)系統(tǒng)所具有的。我們通過(guò)下面的例子說(shuō)明這一點(diǎn)。 利用數(shù)值模擬，討論離散的和連續(xù)的Logistic模型的動(dòng)力學(xué)性質(zhì)。二、實(shí)驗(yàn)內(nèi)容解：先考察連續(xù)Logistic模型的動(dòng)力學(xué)行為。運(yùn)用微分方程定性理論分析，得知對(duì)任意給定的不同參數(shù)r0,該方程從任意初值出發(fā)的軌跡線都將趨近于平衡點(diǎn)N=1

15、.運(yùn)用數(shù)值模擬方法,可以得到圖:Matlab指令為：r=13.6;fori=1:2funlog=(t,x)r(i)*x*(1-x);tt,xx=ode45(funlog,0,10,0.1);subplot(2,2,1),axis(0,10,0.1,1.5),plot(tt,xx,k),holdonendn=16;t=linspace(0,10,n);fori=1:2rr=r(i);x1=dislog1(rr,0.1,n);subplot(2,2,2),axis(0,10,0.1,1.5),plot(t,x1,-k),holdonx2=dislog2(rr,0.1,n);subplot(2,2,

16、3),axis(0,10,0.1,1.5),plot(t,x2,k),holdonx3=dislog3(rr,0.1,n);subplot(2,2,4),axis(0,10,0.1,1.5),plot(t,x3,k),holdonend 求解子程序的函數(shù)文件為：文件一：function N=dislog1(r,N0,n)N=N0 zeros(1,n-1);x=(exp(r*10/n)-1)*N0 zeros(1,n-1);for i=2:n x(i)=exp(r*10/n)*x(i-1)/(1+x(i-1); N(i)=x(i)/(exp(r*10/n)-1);end文件二：function

17、N=dislog2(r,N0,n)N=N0 zeros(1,n-1);x=r*10/n*N0 zeros(1,n-1);for i=2:n x(i)=exp(r*10/n)*x(i-1)*exp(-x(i-1); N(i)=x(i)/(r*10/n);end文件三：function N=dislog3(r,N0,n)N=N0 zeros(1,n-1);x=r*10/n*N0/(r*10/n+1) zeros(1,n-1);for i=2:n x(i)=(r*10/n+1)*x(i-1)*(1-x(i-1); N(i)=x(i)*(r*10/n+1)/(r*10/n);end由左上圖可以看到，對(duì)

18、不同的參數(shù)值r=1，r=3.6，微分方程解軌線的變化特征是一樣的，都會(huì)從原點(diǎn)附近離開(kāi)，單調(diào)增加趨于1.接著，以三種不同格式將微分方程離散化，為了比較離散方程的動(dòng)力學(xué)行為，統(tǒng)一取離散步長(zhǎng)為：在圖中，每個(gè)子圖橫軸變量為t，縱軸變量為N=N（t）.（1） 如果在區(qū)間t，t+t上，對(duì) 積分，分離變量得由此得到記則得到差分方程因?yàn)樵谵D(zhuǎn)化的過(guò)程中沒(méi)有附加任何限制，這個(gè)離散模型僅僅描述了連續(xù)的Logistic模型在給定的離散時(shí)間點(diǎn)上的動(dòng)態(tài)，所以它應(yīng)該具有與連續(xù)模型同樣的動(dòng)態(tài)行為。數(shù)值實(shí)驗(yàn)驗(yàn)證了這一點(diǎn)，即右上圖。（2） 假設(shè)在t，t+t上，單位變化率保持不變，即在t，t+t上積分，得到記得到微分方程當(dāng)t=1

19、時(shí)，稱(chēng)為Ricker模型。該模型常被用于描述魚(yú)群的變化。由于這個(gè)離散模型在小區(qū)間上部分改變了連續(xù)模型，當(dāng)參數(shù)r變化時(shí)，它會(huì)表現(xiàn)出不同于連續(xù)系統(tǒng)的性質(zhì)，如左下圖。（3） 假設(shè)在t，t+t上，變化率dN/dt保持不變，即按最簡(jiǎn)單的歐拉差分格式離散化，記得到差分方程當(dāng)t=1時(shí)，稱(chēng)為離散的Logistic模型。這個(gè)離散模型在小區(qū)間上完全改變了連續(xù)模型。當(dāng)參數(shù)r變化時(shí)，它也會(huì)表現(xiàn)出不同于連續(xù)系統(tǒng)的性質(zhì)，如上面截圖右下圖。在上面的程序中，我們?nèi)∧M步數(shù)n=16,相當(dāng)于將模擬時(shí)間區(qū)間0,10分成步長(zhǎng)為t=10/16的時(shí)間間隔，結(jié)果對(duì)英語(yǔ)較小的參數(shù)值r=1，不同形式的離散格式?jīng)]有導(dǎo)致不同的變化；而對(duì)于較大的參數(shù)值r=3.6，不同形式的離散格式導(dǎo)致了不同的變化。如果取模擬步數(shù)n=100，時(shí)間步長(zhǎng)縮短為t=0.1，則3種離散格式對(duì)于不同的參數(shù)值r=1和r=3.6都能保持原來(lái)系統(tǒng)的動(dòng)態(tài)變化性質(zhì)。可見(jiàn)只要限制離散步長(zhǎng)，采用不同格式的離散化的系統(tǒng)仿真都能保持連續(xù)系統(tǒng)的特征。最