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

1、實驗08穩(wěn)定性模型(4學時)(第7章穩(wěn)定性模型)1.(驗證)捕魚業(yè)的持續(xù)收獲一一產量模型p215219產量模型：,_Xx-t)=F(x)=rx1-一一ExINJ其中，x(t)為t時刻漁場中的魚量。r是固有增長率。N是環(huán)境容許的最大魚量。E是捕撈強度，即單位時間捕撈率。要求:運行下面的m文件，并把相應結果填空，即填入“%7.1捕魚業(yè)的持續(xù)收獲一一產量模型%文件名：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)建符號表達式x=solve(Fx,x)%求解F(x)=0(求根)%得到兩個平衡點，記為：%x0=,x1=,見P216(4)式x0=x(2);x1=x(1);%符號變量x的結構類型成為<2Msym>%求F(x)的微分F'(x)symsx;%定義符號變量x的結構類型為<1Msym>dF=diff(Fx,'x');%求導dF=simple(dF)%簡化符號表達式%得F'(x)=%求F'(

3、x0)并簡化dFx0=subs(dF,x,x0);%將x=x0代入符號表達式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)定(根據平衡點穩(wěn)定性的準則)；%若E>r,則結果正好相反。%在漁場魚量穩(wěn)定在x0的前提下(E<r),求E使持續(xù)產量h(x0M到最大hm%通過分析(見教材p216圖1),只需求x0*使f(x)達到最大，且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)式%產量模型的結論是：%將捕撈率控制在固有增長率的一半(E=r/2)時，能夠獲得最大的持續(xù)產量。符號簡化函數simple,變量替換函數sub的用法見提示。給出填空后的M文件(見215217):%7.1捕魚業(yè)的持續(xù)收獲一一產量模型%文件名：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)建符號表達式x=solve(Fx,x)%求解F(x)=0(求根)%得到兩個平衡點，記為：%x0=-N*(-r+E)/r,x1=_0_,見P216(4)式x0=x(2);x1=x(1);%符號變量x的結構類型成為<2Msym>%求F(x)的微分F'(x)symsx;%定義符號變量x的結構類型為<

6、1Msym>dF=diff(Fx,'x');%求導dF=simple(dF)%簡化符號表達式%得F'(x)=r-2*r*x/N-E%求F'(x0)并簡化dFx0=subs(dF,x,x0);%將x=x0代入符號表達式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)定(根據平衡點穩(wěn)定性的準則)；%若E>r,則結果正好相反。%在漁場魚

7、量穩(wěn)定在x0的前提下(E<r),求E使持續(xù)產量h(x0M到最大hm%通過分析(見教材p216圖1),只需求x0*使f(x)達到最大，且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)式%產量模型的結論是：%將捕撈率控制在固有增長率的一半(E=r/2)時，能夠獲得最大的持續(xù)產量。2.(驗證、編程)種群的相

8、互競爭P222228模型：xi(t)=riXi(1-ai-)NiN2XiX2X2(t)=2X2(1-02)NiN2其中，xi(t),X2(t)分別是甲乙兩個種群的數量。ri,r2是它們的固有增長率。Ni,N2是它們的最大容量。E：單位數量乙(相對N2)消耗的供養(yǎng)甲的食物量為單位數量甲(相對N1)消耗的供養(yǎng)甲的食物量的E倍。對d2可作相應解釋。2.1 (編程)穩(wěn)定性分析p224225要求：補充如下指出的程序段，然后運行該m文件，對照教材上的相應結果%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)%求方程的平衡點，即解代數方程組(見P224的(5)式)%f(x1,x2)=0%g(x1,x2)=0編寫出該程序段。提示(1)使用符號表達式；(2)用函數solve求解代數方程組，解放入x1,x2;(3)調整解(平衡點)的順序放入P中(見下面注釋所示)，P的結構類型為<4X2sym>,P的第1列對應x1,第2列對應x2。x1x2=x1,x2%顯示結果disp('');P%調整位置后的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%顯示結果p

11、=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);disp('');Ppq%顯示結果%得到教材p225表1的前3列，經測算可得該表的第4列，即穩(wěn)定條件。補充后的程序和運行結果(見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)%求方程的平衡點，即解代數方程組(見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%顯示結果disp('');P%調整位置后的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%顯示結果p=subs(-(fx1+g

14、x2),x1,x2,P(:,1),P(:,2);%替換p=simple(p)%簡化符號表達式pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp('');Ppq%顯示結果%得到教材p225表1的前3列，經測算可得該表的第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)(驗證)當x1(0)=x2(0)=0.1時，求微分方程的數值解，將解的數

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)的一組函數值,x(:,2)對應x2(t)gridon;axis(0802);text(2.4,1.55,'x1(t)');text(2.2,0.25,'x2(t)');%函數文件名：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)運行程序并給出結果(比較227圖2(a):10(2)(驗證)將xi(0)=1,X2(0)=2代入(1)中的程序并運行。給出結果(比較227圖2(b):(3)(編程)在同一圖形窗口內，畫(1)和(2)的相軌線，相軌線是以xi(t)為橫坐標，X2(t)為縱坐標所得到的一條曲線。具體要求參照下圖。圖中的兩條“點線”直線，一條的兩個端點為(0,1)和(1,0),另一條的兩個端點為(0,2刑(1.6,0)。11(3)給出程序及其運行結果(比較圖!圖2(c):%命令文件ts

17、=0:0.2:8;x0=0.1,0.1;t,x=ode45cp227',ts,x0);plot(x(:,1),x(:,2);%x(:1)是x1(t)的組函數值，x(:,2)對應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.(編程)種群的相互依存一

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

19、1/N1-x2/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%顯示結果

20、j匚4|bbj1ndif!PX|匕I.fiai/噸"i#duh*?iJid"rIfvll?十HL屯r)irS-IdTrN11心叱-(kL-Dj舊產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»|134.(驗證)食餌-捕食者模型p230232模型的方程：X(t)=rx-axyx(0)=xy(t)7ybxyy(0)=y要求：設r=1,d=0.5,a=0.1,b=0.02,xo=25,y

21、o=2。輸入p231的程序并運行，結果與教材p232的圖1和圖2比較。給出2個M文件(見231)和程序運行后輸出的圖形(比較232圖1、2):函數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)'),%運行中在圖上標注pause,plot(x(

22、:,1),x(:,2),grid,x(t),y(t)圖形：14相軌線y(x)圖形:5.(驗證)差分形式的阻滯增長模型p236242阻滯增長模型用微分方程描述為：x(t)=”(1一字也可用差分方程描述為:k=°12jH上式可簡化為一階非線性差分方程：15Xk.1=b4(1-Xk),k=0,1,2,|考察給定b和x0值后，當k一訓寸，xk的收斂情況(實際上取k足夠大就可以了)。5.1數值解法和圖解法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的值，顯示X81X1

23、00的值。觀察收斂與否。(結果與教材p238239表1比較)下面是計算程序，在兩處下劃線的位置填入滿足要求的內容。%不同b值下方程x(k+1)=b*x(k)*(1-x(k),k=0,1,2,.的計算結果%文件名：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

24、位有效數字，NaN表示不是數值(1)對程序正確填空，然后運行。填入的正確語句和程序的運行結果(對應不同的b值)見238表1:x(1,:)=b*x0*(1-x0);x(k+1,:)=b.*x(k,:).*(1-x(k,:);16CouandVin.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.154

25、0.4T網Ol44750.G40St).5死鴕Cl41ISD.61540.82360,8530.83170.羽，S370.11130,47940.43270.37030.轉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,

26、S23CG.陰M0.8S1T。冏柒第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)的數據對照，注意收斂的倍周期。%圖解法，圖1和圖2%文件名：p238fig_1_2.mclear;clc;

27、closeall;f=(x,b)b.*x.*(1-x);%定義f是函數的句柄，函數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);%指定圖形窗口figuren17subplot(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;

28、%畫迭代軌跡線fork=2:5x1=x(k,n);y1=x(k,n);plot(x0+i*y0,x0+i*y1,x1+i*y1,'r');%實部為橫坐標，虛部為縱坐標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

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

30、=3.55)1.0.80.60.40.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

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

32、171;二憶.士口信*a士1*ht*t09Q8Q7N，K0.60.5n4133363IS3.43S65.3收斂、分岔和混沌p240-242畫出教材p241圖4模型的收斂、分岔和混沌21程序的m文件如下:%圖4模型(6)的收斂、分岔和混沌%文件名：p241fig4.m%模型：xk+1=b*xk*(1-xk),k=0,1,2,.clear;clc;closeall;b=2.5:0.0001:4;%參數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,'

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

34、eEditVie*lasertloolsIesJrtottindo*Help5日d)k盟總zi0"圉日模型的收被、分岔和混他U,>111J；：!:；:“；-":"I-.一a1>e£123(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&l

35、t;3.544):>FlgUEE1fZT|XEileIditInsertIoo>ls取drt叮tindwrHelp日白U)珞-：A-QF®配3ns"0圉日模型的收斂、分岔和混沌24(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的上限

36、。%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.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-1

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

38、下限B=ones(n+2,1);fori=0:nB(i+2)=fun(B(i+1),i);end%單周期收斂時，B(1)<b<B(2);2倍時，B(2)<b<B(3);disp(num2str(-1:n)',B,10);%10位有效數字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最后的2X2個值whiled-c>

39、;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浸防個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%范數，比較c=b;%滿足2An倍周期收斂的b給區(qū)間(c,d)的下限celsed=b;%不滿足2八門倍周期收斂的b給區(qū)間(c,d)的上限dendendbmax=c;%2八門倍周期收斂的b的上限近似27conaiid一.rz-iFn'irx£d七口

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

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

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

43、件下的增氏服從1。機比規(guī)律（見5.6節(jié)的阻滯增長模型），即30*(0=f(x)=rx|-N/(1)是固有增長率，川是環(huán)境容許的最大魚最，用人工）表示單位時間的增長量.2.單位時間的捕撈量（即產量）與漁場魚量工成正比，比例常數E表示單位時間捕撈率,又稱為捕撈強度，可以用比如捕魚網眼的大小或出海漁船數量來控制其大小.于是單位時間的捕撈量為(2)h（x）-Ex根據以上假設并記尸（#）=/（）-A（x）得到捕措情況下漁場魚量滿足的方程x（t）=F（富）=1_木）-Ex我們并不需要解方程（3）以得到工（。的動態(tài)變化過程.同希望知道漁場的穩(wěn)定魚量和保持穩(wěn)定的條件，即時間上足夠長以后漁場魚最工。）的趨向，并

44、由此確定最大持續(xù)產量.為此可以直接求方程（3）的平衡點并分析其穩(wěn)定性.令Ex-0得到兩個平衡點2=N(1-g,J?!=0不難算出F'(不)=E-r.k(盯)-r-E所以若£<r（5）有尸（如）<0,尸（盯）>0,故通點穩(wěn)定百點不穩(wěn)定（判斷平衡點穩(wěn)定性的準則見7.7節(jié)）彳若E>r,則結果正好相反.y-f(x)圖】最大持續(xù)產拼的圖解法xf=M2xDNE是捕撈率，r是最大的增長率，上述分析表明，只要捕撈適度（£<1）.就可使?jié)O場魚量穩(wěn)定在質，從而獲得持續(xù)產最M不）二取°而當捕撈過度時（E>力，漁場魚量將趨向盯=0,當然談不上獲

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

46、增長率r的一半，更簡單一些，可以說使?jié)O場魚量保持在最大魚量N的一半時，能夠獲得最大的持續(xù)產量.效益模型從經濟角度看不應追求產量最大，而應考慮效益最佳.如果經濟效益用從捕撈所得的收入中扣除開支后的利潤來衡量，并且簡單地假設：魚的銷售單價為常數P,單位捕撈率（如每條出海漁船）的費用為常數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）達到最大的捕撈強度為(12)將4代入（4）式，

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

48、致捕撈過度，下面討論這個模型.（11）式給出了利潤與捕撈強度的關系令R（E）=0的解為跖，可得(15)當E風時，利潤靛（E）0,盲目的經營者們會加大捕措強度；若EaE-利潤R（E）<0*他們當然要減小強度,所以乙是盲目捕撈下的臨界強度圖2盲目捕撈強度的圖解法M也可由圖解法確定.在圖2中以£為橫坐標，按（11）式畫出T（E）和S（E）,它們交點的橫坐標即為名（圖2中的號1或后加）,由（15）式或圖2容易知道，酌存在的必要條件（即4>。）是即售價大于（相對于總量而言）成本.并且由（15）式可知，成本越低，售價越高，則瓦越大,將（15）代人（4）式，得到盲目捕撈下的漁場穩(wěn)定魚量

49、為（17）%完全由成本-價格比決定，隨著價格的上升和成本的下降%將迅速減少，出現捕撈過度.比校32）和（15）式可知=2號,即盲目捕撈強度比最大效益下捕撈噩度大一倍.從（15）式和圖2還可以得到，當竟<p<2時，（昂<）/</，如圖2中凰，稱經濟學捕撈過度；當P>2右時，曷>£，如圖2中”，稱生態(tài)學捕撈過度.21833:*本節(jié)完*評注為了研究漁業(yè)的產量、效益及捕撈過度問題，首先在對魚的自然增長和捕撈情況的合理假設F,建立漁場魚量的基本方程（3）,并利用平衡點穩(wěn)定性分析確定了保持漁場魚鼠楫定的條件.產量、效益和捕撈過度3個模型在穩(wěn)定的前提下步步深入，

50、數學推導過程十分簡單，卻得到了在定性關系上與實際情況完全符合的結果.如果改變對魚的自然增長和人工捕撈的假設，模型及結果將隨之變化（習題1,2）.72軍備競賽兩個國家或兩個國家集團之間由于相互不信任和各種矛盾的存在、發(fā)展而不斷增加自己的軍事力量，防御對方可能發(fā)動的戰(zhàn)爭.能否用一個數學模型描述這種軍備竟賽的過程，從定性和定髭的角度對競賽的結果作出解釋或預測.本節(jié)介紹LF.Richardsgl939年提出的一個模型.當然影響軍備競賽的因素是錯綜復雜的，無法用數學工具給以恰當的圓滿的描述.這個模型只不過告訴我們，一個復雜的實際過程可以被合理地簡化到什么程度，得到的結果又怎樣用來解釋實際現象，并且如果允

51、許軍備比賽沿著機械的本能和固定的模式發(fā)展的話.會導致什么樣的結局廊.模型假設與構成為了方便起見，用軍備這個詞表示軍事力量的總和，如兵力、裝備、軍事頸算等.甲乙雙方在時刻t的軍備分別記作工3）和江。，假定它們的變化只取決于下面3個因素工1,由于相互不信任及矛盾的發(fā)展，一方軍備越大，另一方軍備增加得越快.2,由于各方本身經濟實力的限制.任一方軍備越大，對軍備增長的制約作用越大.3由于相互敵視或領土爭端，每一方都存在膂增加軍備的固有潛力.進一步假定前兩個因素的影響是線性的，第3個因素的影響是常數，那么雙上）和貝。的變化過程可用微分方程組：£（£）=-5+Gy+g（I）y（0三自=

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

53、7.7節(jié)（15）式八當印奴（5）時，平衡點（與，為）是穩(wěn)定的：反之.是不穩(wěn)定的.這就是說，在（5）式的條件下，時間足夠長以后雙方的軍備將分別趨向一個有限值，軍備競賽是穩(wěn)定的、模型的定性解釋根據方程（1）和平衡點穩(wěn)定性的分析，可以解釋幾個簡單而又重要的現象.1.條件（5）表明，當雙方的經濟制約程度邛大于雙方的軍備刺激程度kl時，軍備競賽才會趨向穩(wěn)定.反之聲U）4（f）將趨向無窮，競賽無限地進行下去，可能導致戰(zhàn)爭.2,由（2）式，如果g=A=0,則0=0,兀=0是方程（1）的平衡點，并且在條件（5）下它是穩(wěn)定的，于是如果在某個時候品有=y（/0）=。，%了就永遠保持為E這種情況可以解釋為雙方不存在

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