數(shù)學(xué)建模試驗答案_穩(wěn)定性模型分解

1、實驗08穩(wěn)定性模型(4學(xué)時)(第7章穩(wěn)定性模型)1.(驗證)捕魚業(yè)的持續(xù)收獲一一產(chǎn)量模型p215219產(chǎn)量模型：,_Xx-t)=F(x)=rx1-一一ExINJ其中，x(t)為t時刻漁場中的魚量。r是固有增長率。N是環(huán)境容許的最大魚量。E是捕撈強度，即單位時間捕撈率。要求:運行下面的m文件，并把相應(yīng)結(jié)果填空，即填入“%7.1捕魚業(yè)的持續(xù)收獲一一產(chǎn)量模型%文件名：p215_217.mclear;clc;%無捕撈條件下單位時間的增長量：f(x)=rx(1-x/N)%捕撈條件下單位時間的捕撈量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕撈情況下漁場魚量滿足的方程：x&

2、#39;(t)=F(x)%滿足F(x)=0的點x為方程的平衡點%求方程的平衡點symsrxNE;%定義符號變量Fx=r*x*(1-x/N)-E*x;%創(chuàng)建符號表達(dá)式x=solve(Fx,x)%求解F(x)=0(求根)%得到兩個平衡點，記為：%x0=,x1=,見P216(4)式x0=x(2);x1=x(1);%符號變量x的結(jié)構(gòu)類型成為<2Msym>%求F(x)的微分F'(x)symsx;%定義符號變量x的結(jié)構(gòu)類型為<1Msym>dF=diff(Fx,'x');%求導(dǎo)dF=simple(dF)%簡化符號表達(dá)式%得F'(x)=%求F'(

3、x0)并簡化dFx0=subs(dF,x,x0);%將x=x0代入符號表達(dá)式dFdFx0=simple(dFx0)%得F'(x0)=%求F'(x1)dFx1=subs(dF,x,x1)%得F'(x1)=%若E<r,有F'(x0)<0,F'(x1)>0,故x0點穩(wěn)定，x1點不穩(wěn)定(根據(jù)平衡點穩(wěn)定性的準(zhǔn)則)；%若E>r,則結(jié)果正好相反。%在漁場魚量穩(wěn)定在x0的前提下(E<r),求E使持續(xù)產(chǎn)量h(x0M到最大hm%通過分析(見教材p216圖1),只需求x0*使f(x)達(dá)到最大，且hm=f(x0*)。symsrxNfx=r*x*(1

4、-x/N);%fx=dFdf=diff(fx,'x');x0=solve(df,x)%得x0*=,見P217(6)式hm=subs(fx,x,x0)%得hm=,見P217(7)式%又由x0*=N(1-E/r),可得E*=,見P217(8)式%產(chǎn)量模型的結(jié)論是：%將捕撈率控制在固有增長率的一半(E=r/2)時，能夠獲得最大的持續(xù)產(chǎn)量。符號簡化函數(shù)simple,變量替換函數(shù)sub的用法見提示。給出填空后的M文件(見215217):%7.1捕魚業(yè)的持續(xù)收獲一一產(chǎn)量模型%文件名：p215_217.mclear;clc;%無捕撈條件下單位時間的增長量：f(x)=rx(1-x/N)%捕撈條

5、件下單位時間的捕撈量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕撈情況下漁場魚量滿足的方程：x'(t)=F(x)%滿足F(x)=0的點x為方程的平衡點%求方程的平衡點symsrxNE;%定義符號變量Fx=r*x*(1-x/N)-E*x;%創(chuàng)建符號表達(dá)式x=solve(Fx,x)%求解F(x)=0(求根)%得到兩個平衡點，記為：%x0=-N*(-r+E)/r,x1=_0_,見P216(4)式x0=x(2);x1=x(1);%符號變量x的結(jié)構(gòu)類型成為<2Msym>%求F(x)的微分F'(x)symsx;%定義符號變量x的結(jié)構(gòu)類型為<

6、1Msym>dF=diff(Fx,'x');%求導(dǎo)dF=simple(dF)%簡化符號表達(dá)式%得F'(x)=r-2*r*x/N-E%求F'(x0)并簡化dFx0=subs(dF,x,x0);%將x=x0代入符號表達(dá)式dFdFx0=simple(dFx0)%得F'(x0)=-r+E%求F'(x1)dFx1=subs(dF,x,x1)%得F'(x1)=r-E%若E<r,有F'(x0)<0,F'(x1)>0,故x0點穩(wěn)定，x1點不穩(wěn)定(根據(jù)平衡點穩(wěn)定性的準(zhǔn)則)；%若E>r,則結(jié)果正好相反。%在漁場魚

7、量穩(wěn)定在x0的前提下(E<r),求E使持續(xù)產(chǎn)量h(x0M到最大hm%通過分析(見教材p216圖1),只需求x0*使f(x)達(dá)到最大，且hm=f(x0*)。symsrxNfx=r*x*(1-x/N);%fx=dFdf=diff(fx,'x');x0=solve(df,x)%得x0*=1/2*N,見P217(6)式hm=subs(fx,x,x0)%得hm=1/4*r*N,見P217(7)式%又由x0*=N(1-E/r),可得E*=r/2_,見P217(8)式%產(chǎn)量模型的結(jié)論是：%將捕撈率控制在固有增長率的一半(E=r/2)時，能夠獲得最大的持續(xù)產(chǎn)量。2.(驗證、編程)種群的相

8、互競爭P222228模型：xi(t)=riXi(1-ai-)NiN2XiX2X2(t)=2X2(1-02)NiN2其中，xi(t),X2(t)分別是甲乙兩個種群的數(shù)量。ri,r2是它們的固有增長率。Ni,N2是它們的最大容量。E：單位數(shù)量乙(相對N2)消耗的供養(yǎng)甲的食物量為單位數(shù)量甲(相對N1)消耗的供養(yǎng)甲的食物量的E倍。對d2可作相應(yīng)解釋。2.1 (編程)穩(wěn)定性分析p224225要求：補充如下指出的程序段，然后運行該m文件，對照教材上的相應(yīng)結(jié)果%7.3種群的相互競爭一一穩(wěn)定性分析%文件名：p224_225.mclear;clc;%甲乙兩個種群滿足的增長方程：%x1,(t)=f(x1,x2)=

9、r1*x1*(1-x1/N1-k1*x2/N2)%x2，(t)=g(x1,x2)=r2*x2*(1-k2*x1/N1-x2/N2)%求方程的平衡點，即解代數(shù)方程組(見P224的(5)式)%f(x1,x2)=0%g(x1,x2)=0編寫出該程序段。提示(1)使用符號表達(dá)式；(2)用函數(shù)solve求解代數(shù)方程組，解放入x1,x2;(3)調(diào)整解(平衡點)的順序放入P中(見下面注釋所示)，P的結(jié)構(gòu)類型為<4X2sym>,P的第1列對應(yīng)x1,第2列對應(yīng)x2。x1x2=x1,x2%顯示結(jié)果disp('');P%調(diào)整位置后的4個平衡點：%P(1,:)=P1(N1,0)%P(2,:

10、)=P2(0,N2)%P(3,:)=P3(N1*(-1+k1)/(-1+k2*k1),N2*(-1+k2)/(-1+k2*k1)%P(4,:)=P4(0,0)%平衡點位于第一象限才有意義，故要求P3:k1,k2同小于1,或同大于1。%判斷平衡點的穩(wěn)定性參考教材p245的(18),(19)式。symsx1x2;%重新定義fx1=diff(f,'x1');fx2=diff(f,'x2');gx1=diff(g,'x1');gx2=diff(g,'x2');disp('');A=fx1,fx2;gx1,gx2%顯示結(jié)果p

11、=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2);%替換p=simple(p)%簡化符號表達(dá)式pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp('');Ppq%顯示結(jié)果%得到教材p225表1的前3列，經(jīng)測算可得該表的第4列，即穩(wěn)定條件。補充后的程序和運行結(jié)果(見225表1):2341%7.3種群的相互競爭一一穩(wěn)定性分析%文件名：p224_225.mclear;clc;%甲乙兩個種群滿足的增長方程：%x1'(t)=f(x1,x2)=r1*x1*(1-x1/N1-k1*x2/N2)%x2'

12、(t)=g(x1,x2)=r2*x2*(1-k2*x1/N1-x2/N2)%求方程的平衡點，即解代數(shù)方程組(見P224的(5)式)%f(x1,x2)=0%g(x1,x2)=0%編寫出該程序段。symsx1x2r1r2N1N2k1k2;f=r1*x1*(1-x1/N1-k1*x2/N2);g=r2*x2*(1-k2*x1/N1-x2/N2);x1,x2=solve(f,g,x1,x2);P=x1(2,3,4,1),x2(2,3,4,1);x1x2=x1,x2%顯示結(jié)果disp('');P%調(diào)整位置后的4個平衡點：%P(1,:)=P1(N1,0)%P(2,:)=P2(0,N2)%P

13、(3,:)=P3(N1*(-1+k1)/(-1+k2*k1),N2*(-1+k2)/(-1+k2*k1)%P(4,:)=P4(0,0)%平衡點位于第一象限才有意義，故要求P3:k1,k2同小于1,或同大于1%判斷平衡點的穩(wěn)定性參考教材p245的(18),(19)式。symsx1x2;%重新定義fx1=diff(f,'x1');fx2=diff(f,'x2');gx1=diff(g,'x1');gx2=diff(g,'x2');disp('');A=fx1,fx2;gx1,gx2%顯示結(jié)果p=subs(-(fx1+g

14、x2),x1,x2,P(:,1),P(:,2);%替換p=simple(p)%簡化符號表達(dá)式pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp('');Ppq%顯示結(jié)果%得到教材p225表1的前3列，經(jīng)測算可得該表的第4列，即穩(wěn)定條件2.2(驗證、編程)計算與驗證p227微分方程組,、“x1x2、X1=r1x1(1-f-)jNiN2x1x2X2(t)=摟2(1-02)NiN2、=0.5,二2二1.6,r1=2.5,r2=1.8,N1=1.6,N2=1(1)(驗證)當(dāng)x1(0)=x2(0)=0.1時，求微分方程的數(shù)值解，將解的數(shù)

15、值分別畫出x1(t)和x2(t)的曲線，它們同在一個圖形窗口中。程序：%命令文件ts=0:0.2:8;x0=0.1,0.1;t,x=ode45cp227',ts,x0);plot(t,x(:,1),t,x(:,2);%x(:1)是x1(t)的一組函數(shù)值,x(:,2)對應(yīng)x2(t)gridon;axis(0802);text(2.4,1.55,'x1(t)');text(2.2,0.25,'x2(t)');%函數(shù)文件名：p227.mfunctiony=p227(t,x)k1=0.5;k2=1.6;r1=2.5;r2=1.8;N1=1.6;N2=1;y=r1

16、*x(1)*(1-x(1)/N1-k1*x(2)/N2),r2*x(2)*(1-k2*x(1)/N1-x(2)/N2)'(1)運行程序并給出結(jié)果(比較227圖2(a):(2)(驗證)將xi(0)=1,X2(0)=2代入(1)中的程序并運行。給出結(jié)果(比較227圖2(b):(3)(編程)在同一圖形窗口內(nèi)，畫(1)和(2)的相軌線，相軌線是以xi(t)為橫坐標(biāo)，X2(t)為縱坐標(biāo)所得到的一條曲線。具體要求參照下圖。圖中的兩條“點線”直線，一條的兩個端點為(0,1)和(1,0),另一條的兩個端點為(0,2刑(1.6,0)。(3)給出程序及其運行結(jié)果(比較圖!圖2(c):%命令文件ts=0:0

17、.2:8;x0=0.1,0.1;t,x=ode45cp227',ts,x0);plot(x(:,1),x(:,2);%x(:1)是x1(t)的組函數(shù)值，x(:,2)對應(yīng)x2(t)text(0.03,0.3,'(0.1,0.1)');holdon;x0=1,2;t,x=ode45cp227',ts,x0);plot(x(:,1),x(:,2);text(1.02,1.9,'(1,2)');plot(0,1,1,0,':r',0,1.6,2,0,':r');%畫兩條直線gridon;3.(編程)種群的相互依存一一穩(wěn)定性

18、分析P228229模型:X1(t)"喘域)XiX2、X2(t)=32(-102)NiN2其中,X1（t）,X2（t）分別是甲乙兩個種群的數(shù)量。ri,r2是它們的固有增長率。Ni,N2是它們的最大容量。01：單位數(shù)量乙（相對N2）提供的供養(yǎng)甲的食物量為單位數(shù)量甲（相對N1）消耗的供養(yǎng)甲的食物量的E倍。對d2可作相應(yīng)解釋。要求：修改題2.1的程序，求模型的平衡點及穩(wěn)定性。給出程序及其運行結(jié)果（比較2291表1,注：只要最終結(jié)果）：clear;clc;symsx1x2r1r2N1N2k1k2;f=r1*x1*(1-x1/N1+k1*x2/N2);g=r2*x2*(-1+k2*x1/N1-x

19、2/N2);x1,x2=solve(f,g);P=x1(2,4,1,3),x2(2,4,1,3);symsx1x2;%Uffcr義fx1=diff(f,'x1');fx2=diff(f,'x2');gx1=diff(g,'x1');gx2=diff(g,'x2');A=fx1,fx2;gx1,gx2;p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2);%替換p=simple(p)%簡化#pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);Ppq%顯示結(jié)果j匚4|bb

20、j1ndif!PX|匕I.fiai/噸"i#duh*?iJid"rIfvll?十HL屯r)irS-IdTrN11心叱-(kL-Dj舊產(chǎn)2-7r4rl腔-ME*k>r：)/fJi+UE1J*W+r拄皿-L)>/(k+k?-L)10,0,r!-ri*-rl»r2)Q-IE,kl'rl-r：-rL-rl>fT2r(kl-1)1fit»|4.(驗證)食餌-捕食者模型p230232模型的方程：X(t)=rx-axyx(0)=xy(t)7ybxyy(0)=y要求：設(shè)r=1,d=0.5,a=0.1,b=0.02,xo=25,yo=2。輸入p2

21、31的程序并運行，結(jié)果與教材p232的圖1和圖2比較。給出2個M文件(見231)和程序運行后輸出的圖形(比較232圖1、2):函數(shù)M文件：functionxdot=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=(r-a*x(2).*x(1);(-d+b*x(1).*x(2);命令M文件：ts=0:0.1:15;x0=25,2;t,x=ode45fshier',ts,x0);t,x,plot(t,x),grid,gtext('x(t)'),gtext('y(t)'),%運行中在圖上標(biāo)注pause,plot(x(:,1),x(:

22、,2),grid,x(t),y(t)圖形：相軌線y(x)圖形:5.(驗證)差分形式的阻滯增長模型p236242阻滯增長模型用微分方程描述為：x(t)=”(1一字也可用差分方程描述為:k=°12jH上式可簡化為一階非線性差分方程：Xk.1=b4(1-Xk),k=0,1,2,|考察給定b和x0值后，當(dāng)k一訓(xùn)寸，xk的收斂情況(實際上取k足夠大就可以了)。5.1數(shù)值解法和圖解法p238240(1)取xo=0.2,分別取b=1.7,2.6,3.3,3.45,3.55,3.57對方程Xki=bXk(1-Xk),k=0,1,2,111計算出X1X100的值，顯示X81X100的值。觀察收斂與否。

23、(結(jié)果與教材p238239表1比較)下面是計算程序，在兩處下劃線的位置填入滿足要求的內(nèi)容。%不同b值下方程x(k+1)=b*x(k)*(1-x(k),k=0,1,2,.的計算結(jié)果%文件名：p238table_1.mclc;clearall;formatcompactb=1.7,2.6,3.3,3.45,3.55,3.57;x=zeros(100,length(b);x0=0.2;%初值;%此處填入一條正確語句fork=1:99;%此處填入一條正確語句endK=(81:100)'%將取81100行disp(num2str(NaN,b;K,x(K,:),4);%取4位有效數(shù)字，NaN表示不

24、是數(shù)值(1)對程序正確填空，然后運行。填入的正確語句和程序的運行結(jié)果(對應(yīng)不同的b值)見238表1:x(1,:)=b*x0*(1-x0);x(k+1,:)=b.*x(k,:).*(1-x(k,:);CouandVin.dovE3叵區(qū)1.91.72.63.33.4E3.553.57310.4113D.6L5447940,4475-0.MS0.475gS20.U18CU0.82360.8530.88740.B903出0.411M口54Qt4794U.43Z&.35二U.MSBK40.411K0.8236D.R4SS0,E12T0.nQ6950.4118D.1540.4T網(wǎng)Ol44750.G

25、40St).5死鴕Cl41ISD.61540.82360,8530.83170.羽，S370.1113060,47940.43270.37030.轉(zhuǎn)5口S30.4113U6LE40.82360.82730.E27B39.41IBD.6154Ci.47940.44打CJ.506U.5DF9WH.4118OiLE40.8236凡函mQ.S3740.B9Q2910.41IS61540474Cl4327(J.犯4HCLM二口常0.41130161sd0,825o.三處g0.31270.亞930.41IEDl6L540.47940l447j0.5M50.5607%0.41130,S23CG.陰M0.8S

26、1T。冏柒第0.41ISCk61B40,47940.43270.37030.3TBB960,*11180,61540.8236以總歸0.8278.,84四0.41ISD,6154fi.4740,4471r,EC60.479?第0.4113口540.82560.S&30.距X0.89L的0.4113-G1540.47M0l432?0.354B0.34661000.1113(X61B4Q,鴕蒸0.812T0.3055(2)運行以下程序，觀察顯示的圖形，與題(1)的數(shù)據(jù)對照，注意收斂的倍周期。%圖解法，圖1和圖2%文件名：p238fig_1_2.mclear;clc;closeall;f=(

27、x,b)b.*x.*(1-x);%定義f是函數(shù)的句柄，函數(shù)b*x*(1-x)含兩個變量x,bb=1.7,2.6,3.3,3.45,3.55,3.57;x=ones(101,length(b);x(1,:)=0.2;fork=1:100x(k+1,:)=f(x(k,:),b);endforn=1:length(b)figure(n);%指定圖形窗口figurensubplot(1,2,1)%開始迭代的圖形fplot(x)f(x,b(n),x,0101);%x是自變量，畫曲線y=bx(1-x)和直線y=xaxissquare;holdon;x0=x(1,n);y0=0;%畫迭代軌跡線fork=2:

28、5x1=x(k,n);y1=x(k,n);plot(x0+i*y0,x0+i*y1,x1+i*y1,'r');%實部為橫坐標(biāo)，虛部為縱坐標(biāo)x0=x1;y0=y1;endtitle(圖解法：開始迭代的圖形(b='num2str(b(n)')');holdoff;subplot(1,2,2);%最后迭代的圖形fplot(x)f(x,b(n),x,0101);axissquare;holdon;xy(1:2:41)=x(81:101,n)+i*x(81:101,n);%盡量不用循環(huán)xy(2:2:40)=x(81:100,n)+i*x(82:101,n);plo

29、t(xy,'r');title(圖解法：最后迭代的圖形(b='num2str(b(n)')');holdoff;end(2)運行程序并給出結(jié)果(對應(yīng)不同的b值)見238圖1、2:圖解法：最后迭代的圖形(b=1.7)10.80.60.40.200.5120倍0圖解法：開始迭代的圖形（b=2.6）圖解法：最后迭代的圖形（b=2.6）20倍圖解法：開始迭代的圖形（b=3.3）21倍圖解法：開始迭代的圖形（b=3.45）圖解法：最后迭代的圖形（b=3.45）圖解法：開始迭代的圖形（b=3.55）23倍圖解法：最后迭代的圖形(b=3.55)1.0.80.60.40

30、.2000.510.80.60.40.200.5圖解法：開始迭代的圖形（b=3.57）1圖解法：最后迭代的圖形（b=3.57）10.80.60.40.2000.51混沌5.2b值下的收斂圖形p238240下面程序是在不同b值下的收斂圖形。%b值下的收斂圖形%文件名：p212tab1figure.m%方程(6)xk+1=b*xk*(1-xk),k=0,1,2,.clear;clc;closeall;b=3.33.453.553.57;x=0.2*ones(size(b);%初值x0=0.2fork=1:100x(k+1,:)=b.*x(k,:).*(1-x(k,:);endplot(b,x(82

31、:101,:),'.r','markersize',5);%可改變值5,調(diào)整數(shù)據(jù)點圖標(biāo)的大小xlabel('b');ylabel('x(k)');gridon要求:運行以上程序。在運行結(jié)果的圖形中，從對應(yīng)的b值上的“點數(shù)”判斷倍周期收斂。提示:放大圖形。程序的運行結(jié)果（見238表1）:Fjfijife1Lt|fx1 k工dxt.VlLakm-t"TdibIxDaulctapadj：hHtilpijuil£|-a口”><*4=能國1卜,-：，1h-l-tJs«二憶.士口信*a士1*ht*t0

32、9Q8Q7N，K0.60.5n4133363IS3.43S65.3收斂、分岔和混沌p240-242畫出教材p241圖4模型的收斂、分岔和混沌程序的m文件如下:%圖4模型(6)的收斂、分岔和混沌%文件名：p241fig4.m%模型：xk+1=b*xk*(1-xk),k=0,1,2,.clear;clc;closeall;b=2.5:0.0001:4;%參數(shù)b的間隔取值x=0.2*ones(size(b);xx=;n=100000;fork=1:nx=b.*x.*(1-x);ifk>=n-50xx=xx;x;endendplot(b,xx,'.b',markersize,5)

33、;%可改變值5,調(diào)整數(shù)據(jù)點圖標(biāo)的大小title('圖4模型的收斂、分岔和混沌,);xlabel(，參數(shù)b,);ylabel(,x(k)(k足夠大)，)；gridon;運行以上m文件。運行結(jié)果(比較241圖4):事工1-Zt£aiVawZxixnrlTa:I/D-aijcfcez,IhxndrviHalp白金席%、叫國安£2口國I國朗士悔理的性能，子史和混沌1r，口省.上g.一.一.*15.42n倍周期圖形p240-242要求：在上題的程序中，修改b值，使運行后顯示以下圖形:(1)單周期(1<b<3):EileEditVie*lasertloolsIesJ

34、rtottindo*Help5日d)k盟總zi0"圉日模型的收被、分岔和混他U,>111J；：!:；:“；-":"I-.一a1>e£1(2)2倍周期(3<b<3.449):/FLgUZEt口叵僅.91.1e*Insert工。工與Rudrt<xptindo*Help，k鼠一2立130因口模型的收斂“分岔和您隨3.053.13.153.23.253.33.353.43.4536參抽b576E6了0.0&Doo(1<卷嗖)等(3)22倍周期(3.449<b<3.544):>FlgUEE1fZT|XE

35、ileIditInsertIoo>ls取drt叮tindwrHelp日白U)珞-：A-QF®配3ns"0圉日模型的收斂、分岔和混沌(4)23倍周期(3.544<b<3.564):F£tuie1L叵|區(qū)EditVie*Insert工口o1工上Rudrta;ptindo*HelpLRJU，玲.鼠-T盟總S圖b模型的收散、分茗和混沌5.5（編程）求2n倍周期的b值區(qū)間p240241求2An倍周期U斂的b的上限的程序如下：functionbmax=fun(bn_1,n)%求2入門倍周期收斂的b的上限。%bn_1是2A(n-1)倍周期收斂的b的上限(或2=

36、倍周期收斂的b的下限)c=bn_1;d=3.57;%2倍周期時收斂b>bn_1,b<3.57y=zeros(1,2*2An);%行向量，用于存放xk最后的2X2n個值whiled-c>1e-12%使區(qū)間(c,d)足夠小b=(d+c)/2;x=0.2;%b取區(qū)間的中點fori=1:100000x=b*x*(1-x);endy(1)=x;%取最后2>2八門個xkfork=1:2A(n+1)-1y(k+1)=b*y(k)*(1-y(k);endifnorm(y(1:2An)-y(2An+1:2A(n+1)<1e-12%范數(shù)，比較c=b;%滿足2An倍周期收斂的b給區(qū)間(

37、c,d)的下限celsed=b;%不滿足2八門倍周期收斂的b給區(qū)間(c,d)的上限dendendbmax=c;%2八門倍周期收斂的b的上限近似要求:編寫程序，調(diào)用以上函數(shù)文件，求單倍周期、2倍周期、25倍周期收斂的b的上限近似值，輸出要求有10位有效數(shù)字。運行結(jié)果輸出形式如下:提示：可使用num2str函數(shù)。用下面的程序結(jié)構(gòu)，會使運行速度加快。functionmain()自己編寫的程序；將上面的函數(shù)文件復(fù)制到此處。運行的程序及輸出結(jié)果：functionmain()clc;n=5;%求單(賦0)、2倍(賦1)、2倍周期收斂的b的上、下限B=ones(n+2,1);fori=0:nB(i+2)=f

38、un(B(i+1),i);end%單周期收斂時，B(1)<b<B(2);2倍時，B(2)<b<B(3);disp(num2str(-1:n)',B,10);%10位有效數(shù)字functionbmax=fun(bn_1,n)%求2入門倍周期收斂的b的上限。%bn_1是2A(n-1)倍周期收斂的b的上限(或2=倍周期收斂的b的下限)c=bn_1;d=3.57;%2倍周期時收斂b>bn_1,b<3.57y=zeros(1,2*2An);%行向量，用于存放xk最后的2X2n個值whiled-c>1e-12%使區(qū)間(c,d)足夠小b=(d+c)/2;x=0

39、.2;%b取區(qū)間的中點fori=1:100000x=b*x*(1-x);endy(1)=x;%取最后2>2八門個xkfork=1:2A(n+1)-1y(k+1)=b*y(k)*(1-y(k);endifnorm(y(1:2An)-y(2An+1:2A(n+1)<1e-12%范數(shù)，比較c=b;%滿足2An倍周期收斂的b給區(qū)間(c,d)的下限celsed=b;%不滿足2八門倍周期收斂的b給區(qū)間(c,d)的上限dendendbmax=c;%2八門倍周期收斂的b的上限近似conaiid一.rz-iFn'irx£d七口也。耳3-1102.99976909713.449355

40、96723.54402147133.56435787943.56873999353.569679465A»|注：vpa的用法附1:實驗提示第1題提示：符號簡化函數(shù)simple,變量替換函數(shù)sub符號簡化函數(shù)simple的格式：simple(S)對符號表達(dá)式S嘗試多種不同的算法簡化，以顯示S表達(dá)式的長度最短的簡化形式變量替換函數(shù)sub的格式:subs(S,OLD,NEW)將符號表達(dá)式S中的OLD變量替換為NEW變量附2:第7章穩(wěn)定性模型2157.1捕魚業(yè)的持續(xù)收獲7第/邕穩(wěn)定性模型雖然動態(tài)過程的變化規(guī)律一般要用微分方程建立的動態(tài)模型來描述，但是對于某些實際問題，建模的主要目的并不是要尋

41、求動態(tài)過程每個瞬時的性態(tài)，而是研究某種意義下穩(wěn)定狀態(tài)的特征，特別是當(dāng)時間充分長以后動態(tài)過程的變化趨勢.比如，在什么條件下描述過程的變量會越來越接近某些確定的數(shù)值，在什么情況下又會越來越遠(yuǎn)離這些數(shù)值而導(dǎo)致過程不穩(wěn)定.為了分析這種穩(wěn)定與不穩(wěn)定的規(guī)律，常常不需要求解微分方程（并且我們將會看到，即使對于不太復(fù)雜的方程，解析解也不是總能得到的），而可以利用微分方程穩(wěn)定性理論，直接研究平衡狀態(tài)的穩(wěn)定性就行了.在下面建模過程中，大多數(shù)情況下我們將直接引用微分方程穩(wěn)定性理論的結(jié)果.不熟悉這方面內(nèi)容的讀者可以參閱17節(jié).差分方程的穩(wěn)定性與微分方程有很多相似之處.也有明顯的差別，工5和工7節(jié)將介紹這方面的內(nèi)容.7

42、.1捕魚業(yè)的持續(xù)收獲可持續(xù)發(fā)展是一項基本國策，對于像漁業(yè)、林業(yè)這樣的再生資源，一定要注意適度開發(fā)，不能為了一時的高產(chǎn)去“竭澤而漁”，應(yīng)該在持續(xù)穩(wěn)產(chǎn)的前提下追求產(chǎn)量或效益的最優(yōu)化.考察一個漁場，其中的魚量在天然環(huán)境下按一定規(guī)律增及，如果捕撈量恰好等于增長最，那么漁場魚量將保持不變，這個捕撈量就可以持續(xù)下去本節(jié)要建立在捕播情況下漁場魚鼠遵從的方程，分析魚址穩(wěn)定的條件，并且在穩(wěn)定的前提下討論如何控制捕撈使持續(xù)產(chǎn)量或經(jīng)濟效益達(dá)到最大，最后研究所謂捕撈過度的問題產(chǎn)模型記時刻t漁場中魚量為X。，關(guān)于的自然增長和人工捕撈作如下假設(shè)：L在無捕撈條件下的增氏服從1。機比規(guī)律（見5.6節(jié)的阻滯增長模型），即*(0

43、=f(x)=rx|-N/(1)是固有增長率，川是環(huán)境容許的最大魚最，用人工）表示單位時間的增長量.2.單位時間的捕撈量（即產(chǎn)量）與漁場魚量工成正比，比例常數(shù)E表示單位時間捕撈率,又稱為捕撈強度，可以用比如捕魚網(wǎng)眼的大小或出海漁船數(shù)量來控制其大小.于是單位時間的捕撈量為(2)h（x）-Ex根據(jù)以上假設(shè)并記尸（#）=/（）-A（x）得到捕措情況下漁場魚量滿足的方程x（t）=F（富）=1_木）-Ex我們并不需要解方程（3）以得到工（。的動態(tài)變化過程.同希望知道漁場的穩(wěn)定魚量和保持穩(wěn)定的條件，即時間上足夠長以后漁場魚最工。）的趨向，并由此確定最大持續(xù)產(chǎn)量.為此可以直接求方程（3）的平衡點并分析其穩(wěn)定性

44、.令Ex-0得到兩個平衡點2=N(1-g,J?!=0不難算出F'(不)=E-r.k(盯)-r-E所以若£<r（5）有尸（如）<0,尸（盯）>0,故通點穩(wěn)定百點不穩(wěn)定（判斷平衡點穩(wěn)定性的準(zhǔn)則見7.7節(jié)）彳若E>r,則結(jié)果正好相反.y-f(x)圖】最大持續(xù)產(chǎn)拼的圖解法xf=M2xDNE是捕撈率，r是最大的增長率，上述分析表明，只要捕撈適度（£<1）.就可使?jié)O場魚量穩(wěn)定在質(zhì)，從而獲得持續(xù)產(chǎn)最M不）二取°而當(dāng)捕撈過度時（E>力，漁場魚量將趨向盯=0,當(dāng)然談不上獲得持續(xù)產(chǎn)量了.進(jìn)一步討論漁場魚量穩(wěn)定在與的前提下，如何控制捕撈強度E

45、使持續(xù)產(chǎn)量最大的問甌用圖解法可以非常簡單地得到216；結(jié)果.根據(jù)（1）,（2）式作拋物線耳）和直線=乂*）=&,如圖1.注意到y(tǒng)=/«）在原點的切線為？二%所以在條件（5）下,二株必與有交點P,P的橫坐標(biāo)就是穩(wěn)定平衡點小.根據(jù)假設(shè)2,尸點的縱坐標(biāo)為為穩(wěn)定條件下單位時間的持續(xù)產(chǎn)量.由圖1立刻知道，當(dāng)y=&與在拋物線頂點P*相交時可獲得最大的持續(xù)產(chǎn)量,此時的穩(wěn)定平衡點為(6)且單位時間的最大持續(xù)產(chǎn)量為而由（4）式不難算出保持漁場魚量穩(wěn)定在工；的捕撈率為E=2綜上所述，產(chǎn)量模型的結(jié)論是將捕撈率控制在固有增長率r的一半，更簡單一些，可以說使?jié)O場魚量保持在最大魚量N的一半時，能

46、夠獲得最大的持續(xù)產(chǎn)量.效益模型從經(jīng)濟角度看不應(yīng)追求產(chǎn)量最大，而應(yīng)考慮效益最佳.如果經(jīng)濟效益用從捕撈所得的收入中扣除開支后的利潤來衡量，并且簡單地假設(shè)：魚的銷售單價為常數(shù)P,單位捕撈率（如每條出海漁船）的費用為常數(shù)c,那么單位時間的收入T和支出S分別為T-ph(x)=pEx,S=cE(9)單位時間的利潤為RT-S-pEx-cE在穩(wěn)定條件*=褊下，以（4）代入（10）式，得(10)/?(£)=T(E)-S(£)=pNE(T-(11)用微分法容易求出使利潤K（E）達(dá)到最大的捕撈強度為(12)將4代入（4）式，可得最大利潤下的漁場穩(wěn)定魚量5及單位時間的持續(xù)產(chǎn)量K為(13)(14)2

47、17將（12）14）式與產(chǎn)量模型中的（6）（8）式相比較可以看出,在最大效益原則下捕撈率和持續(xù)產(chǎn)量均有所減少，而渣場應(yīng)保持的穩(wěn)定魚量有所增加，并且減少或增加的部分隨著捕撈成本e的增長而變大，循著銷售價格p的增長而變小.請讀者分析這些結(jié)果是符合實際情況的.捕措過度上面的效益模型是以計劃捕撈（或稱封閉式捕撈）為基礎(chǔ)的，即漁場由單獨的經(jīng)營者有計劃地捕撈，可以追求最大利潤.如果漁場向眾多盲目的經(jīng)營者開放，比如在公海上無規(guī)則地捕撈,那么即使只有微薄的利潤，經(jīng)營者也會蟒擁而去，這種情況稱為盲目捕撈（或開放式捕撈）.這種捕撈方式將導(dǎo)致捕撈過度，下面討論這個模型.（11）式給出了利潤與捕撈強度的關(guān)系令R（E）

48、=0的解為跖，可得(15)當(dāng)E風(fēng)時，利潤靛（E）0,盲目的經(jīng)營者們會加大捕措強度；若EaE-利潤R（E）<0*他們當(dāng)然要減小強度,所以乙是盲目捕撈下的臨界強度M也可由圖解法確定.在圖2中以£為橫坐標(biāo)，按（11）式畫出T（E）和S（E）,它們交點的橫坐標(biāo)即為名（圖2中的號1或后加）,由（15）式或圖2容易知道，酌存在的必要條件（即4>。）是圖2盲目捕撈強度的圖解法即售價大于（相對于總量而言）成本.并且由（15）式可知，成本越低，售價越高，則瓦越大,將（15）代人（4）式，得到盲目捕撈下的漁場穩(wěn)定魚量為（17）%完全由成本-價格比決定，隨著價格的上升和成本的下降%將迅速減少，

49、出現(xiàn)捕撈過度.比校32）和（15）式可知=2號,即盲目捕撈強度比最大效益下捕撈噩度大一倍.從（15）式和圖2還可以得到，當(dāng)竟<p<2時，（昂<）/</，如圖2中218凰，稱經(jīng)濟學(xué)捕撈過度；當(dāng)P>2右時，曷>£，如圖2中”，稱生態(tài)學(xué)捕撈過度.:*本節(jié)完*評注為了研究漁業(yè)的產(chǎn)量、效益及捕撈過度問題，首先在對魚的自然增長和捕撈情況的合理假設(shè)F,建立漁場魚量的基本方程（3）,并利用平衡點穩(wěn)定性分析確定了保持漁場魚鼠楫定的條件.產(chǎn)量、效益和捕撈過度3個模型在穩(wěn)定的前提下步步深入，數(shù)學(xué)推導(dǎo)過程十分簡單，卻得到了在定性關(guān)系上與實際情況完全符合的結(jié)果.如果改變對魚

50、的自然增長和人工捕撈的假設(shè)，模型及結(jié)果將隨之變化（習(xí)題1,2）.72軍備競賽兩個國家或兩個國家集團之間由于相互不信任和各種矛盾的存在、發(fā)展而不斷增加自己的軍事力量，防御對方可能發(fā)動的戰(zhàn)爭.能否用一個數(shù)學(xué)模型描述這種軍備竟賽的過程，從定性和定髭的角度對競賽的結(jié)果作出解釋或預(yù)測.本節(jié)介紹LF.Richardsgl939年提出的一個模型.當(dāng)然影響軍備競賽的因素是錯綜復(fù)雜的，無法用數(shù)學(xué)工具給以恰當(dāng)?shù)膱A滿的描述.這個模型只不過告訴我們，一個復(fù)雜的實際過程可以被合理地簡化到什么程度，得到的結(jié)果又怎樣用來解釋實際現(xiàn)象，并且如果允許軍備比賽沿著機械的本能和固定的模式發(fā)展的話.會導(dǎo)致什么樣的結(jié)局廊.模型假設(shè)與構(gòu)

51、成為了方便起見，用軍備這個詞表示軍事力量的總和，如兵力、裝備、軍事頸算等.甲乙雙方在時刻t的軍備分別記作工3）和江。，假定它們的變化只取決于下面3個因素工1,由于相互不信任及矛盾的發(fā)展，一方軍備越大，另一方軍備增加得越快.2,由于各方本身經(jīng)濟實力的限制.任一方軍備越大，對軍備增長的制約作用越大.3由于相互敵視或領(lǐng)土爭端，每一方都存在膂增加軍備的固有潛力.進(jìn)一步假定前兩個因素的影響是線性的，第3個因素的影響是常數(shù)，那么雙上）和貝。的變化過程可用微分方程組：£（£）=-5+Gy+g（I）y（0三自=By+%表示，其中的系數(shù)均大于或等于0.是對方軍備剌激程度的度最e,©

52、是己方經(jīng)濟實力制約程度的度量灑力是巴方軍備競賽的固有潛力-如果我們感興趣的是軍備競賽的結(jié)局由什么因素決定，而不關(guān)心競賽的過程.那么只需用微分方程穩(wěn)定性理論討論時間充分長以后工（。，¥”）的變化趨勢，即方程（1）的平衡點的穩(wěn)定情況令（1）式右端等于0,容易算出平衡點（與，九）為fl印-H'九一4一奴方程（1）的系數(shù)矩陣為F-a*1A=>，p于是按照判斷平衡點穩(wěn)定性的方法計算（見7.7節(jié)（9）-（13）式）p=-（an+a12）=a>0（3）qdetA=ap-kl（4）由穩(wěn)定性準(zhǔn)則（見7.7節(jié)（15）式八當(dāng)印奴（5）時，平衡點（與，為）是穩(wěn)定的：反之.是不穩(wěn)定的.這就

53、是說，在（5）式的條件下，時間足夠長以后雙方的軍備將分別趨向一個有限值，軍備競賽是穩(wěn)定的、模型的定性解釋根據(jù)方程（1）和平衡點穩(wěn)定性的分析，可以解釋幾個簡單而又重要的現(xiàn)象.1.條件（5）表明，當(dāng)雙方的經(jīng)濟制約程度邛大于雙方的軍備刺激程度kl時，軍備競賽才會趨向穩(wěn)定.反之聲U）4（f）將趨向無窮，競賽無限地進(jìn)行下去，可能導(dǎo)致戰(zhàn)爭.2,由（2）式，如果g=A=0,則0=0,兀=0是方程（1）的平衡點，并且在條件（5）下它是穩(wěn)定的，于是如果在某個時候品有=y（/0）=。，%了就永遠(yuǎn)保持為E這種情況可以解釋為雙方不存在任何敵視和爭端，通過裁軍可以達(dá)到持久和平.兩個友好的鄰國正是這樣.3 .如果即使由于

54、某種原因（如裁軍協(xié)定）在某個時候雙方軍備大減，不妨設(shè)M%）=y（r0）=。，那么因為x=h,也將使雙方重整軍備.這說明未經(jīng)和解的裁軍（即不消除敵視或領(lǐng)土爭端）是不會持久的.4 .如果由于某種原因（如戰(zhàn)敗或協(xié)議），在某個時候一方的軍備大減，不妨設(shè)立色）=0,那么因為五二行+心也將使該方重整軍備.這說明存在不信任（無*0）或固有爭端3*0）的單方面裁軍也不會持久.模型參數(shù)的估計為了利用（5）式判斷軍備競賽是否會趨于穩(wěn)定，需要知道白班,AN的數(shù)直估計這些參數(shù)無疑是很困難的，下面是RichardQi提出的一種方法.1. 的估計設(shè)工（0）=0,當(dāng)二較小時，忽略g和-ox的作用,并近似地假定丁二為不變，由方程（1）得i二%220如果當(dāng)t=丁時，則由（6）式得到（7）這說明A"是甲方軍備從0到趕上乙方軍備外所需的時間.例如，譙國從1933年開始重整軍備，只用了約3年的時間就趕上了它的鄰國.假設(shè)它增加軍備的固有潛力g被制約效應(yīng)ax所抵消，那么可以認(rèn)為德國的A-=3年，即上=0.3./可以類似地估計，或者合理地假定它與國家的經(jīng)濟實力成正比.這樣若乙國的經(jīng)濟實力是德國的2倍，則可以估計20.6.嘉a,尸的估計設(shè)g=0,y=0,由方程（】）可得x（t）=.（0”7以代人，算出這表示儀7是在乙方無軍備時甲方軍備減少到原來的,所需的時間.Richards