西農 建模實驗六

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

2、實驗內容解 假設（a） 兔子每經(jīng)過一個月底就增加一個月齡；（b） 月齡大于等于2的兔子都具有繁殖能力；（c） 具有繁殖能力的兔子每一個月一定生一對兔子；（d） 兔子不離開群體（不考慮死亡） 記第n個月初的幼兔（一月齡兔）數(shù)量為a0(n）,成兔（月齡大于等于2）數(shù)量為a1(n），則兔子總數(shù)為a(n)= a0(n）+a1(n）,平衡關系為： 建立模型：這個一階差分方程的矩陣表達式為其中， 利用迭代方法求數(shù)值解，也就是按時間步長法仿真種群增長的動態(tài)過程，模擬幼兔和成兔占整體比例隨時間的變化。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ù)值模擬結果可見，兔子數(shù)量遞增，但是幼兔和成兔在種群中所占比例很快會趨于一個極限。由線性差分方程的定性理論可知，這個極限就是對應于差分方程系數(shù)矩陣A主特征值得歸一化特征向量。因為A是非負矩陣，由矩陣理論知，A主特征值是正實數(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ù)值計算可知，系數(shù)矩陣模A最大的特征值r=1.618，生物上稱之為種群的內稟增長率，是個大于一的實數(shù)。因此種群數(shù)量隨時間遞增。相應的歸一化的特征向量的兩個分量0.382和0.618正是幼兔和成兔在種群中所占比例趨近的穩(wěn)定值，生物上稱之為種群的穩(wěn)定分布。 從這個例子的討論可見，數(shù)值模擬能幫我們對系統(tǒng)的變化有更直觀的認識。但是只有通過計算方程組系數(shù)矩陣的特征值和特征向量，運用差分方程定性分析才能對解得漸進性質給出確定的結論。微分方程模型一、實驗題目 藍鯨的內稟增長率每年估計為5%，估計藍鯨的最大環(huán)境載量為150000條，磷蝦是藍鯨喜歡的一種食物。磷蝦的最大飽

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

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

7、一步，為了便于數(shù)值模擬，保持數(shù)量級的協(xié)調，把鯨魚的單位改為以千條為單位。當初始狀態(tài)為藍鯨5千條、磷蝦750噸/英畝時，動態(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(時間/年),ylabel(種群數(shù)量/千條); gtext(虛線表示藍鯨

8、，實線表示磷蝦);由圖可見，藍鯨的數(shù)量隨著時間的增加而逐漸增加，磷蝦的數(shù)量隨時間的增加而逐漸減少，最后都分別趨近一個穩(wěn)定值。（2）由常微分方程理論知，方程組的常數(shù)值解為方程的平衡點，由平衡點的穩(wěn)定性確定了方程解的動態(tài)性質，即解的漸近行為。上一問的數(shù)值結果表明，方程組具有一個穩(wě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(平衡點是：);平衡點是： disp(vpa(y1steady,6) vpa(y2steady,6); 0, 0 150.0, 0 0, 500.0 181.034, 258.621接著，考慮方程組在平衡點（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得到的結果表明，在正平衡點(181.034，258.621)，相應的兩個特征值的實部都是負的。按照微分方程定性理論可知，方程組正平衡點(181.034，258.621)是漸近穩(wěn)定的，即從任意初值點出發(fā)，解軌線都會趨近該點。因此，可以斷言，只要停止捕撈，藍鯨與磷蝦可以長期共存。(3) 為了確定當藍鯨數(shù)量為平衡態(tài)數(shù)量的r%，磷蝦數(shù)量為平衡態(tài)時，停止

11、捕撈，藍鯨恢復到平衡態(tài)需要的時間，只有通過數(shù)值模擬。對不同的初值N1(0)=r/100181.034,N2(0)=258.621,在一個充分長的時間區(qū)間上求方程組的數(shù)值解，注意到藍鯨數(shù)量會遞增趨于平衡態(tài)，可以N1(T)181.034-0.001確定的時間T近似表示種群恢復所需的時間，得到對應的函數(shù)關系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(時間),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(恢復時間T),gridgtext(實線表示藍鯨，虛線表示磷蝦)由圖可知，得到的函數(shù)關系T=T(r)是不光滑的。注意到計算機只會計算離散值，盡管ode45采用了四階、五階龍格-庫塔法，可能是離散的時

13、間步長太大，計算得結果并未反映連續(xù)系統(tǒng)的真實規(guī)律。重新規(guī)定步長，計算量大些，但能夠保證離散結果逼近連續(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(恢復時間T),grid果然得到理想的效果。在利用計算機進行模型的數(shù)值求解時，對模型解的性質先有一個定性分析是非常重要的

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

15、.運用數(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由左上圖可以看到，對

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

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