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

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

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

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

4、=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é)論是：%將捕撈率控制在固有增長(zhǎng)率的一半(E=r/2)時(shí)，能夠獲得最大的持續(xù)產(chǎn)量。符號(hào)簡(jiǎn)化函數(shù)simple,變量替換函數(shù)sub的用法見提示。給出填空后的M文件(見215217):%7.1捕魚業(yè)的持續(xù)收獲一一產(chǎn)量模型%文件名：p215_217.mclear;clc;%無(wú)捕撈條件下單位時(shí)間的增長(zhǎng)量：f(x)=rx(1-x/N)%捕撈條件下單位時(shí)間的捕撈

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

6、dF=diff(Fx,'x');%求導(dǎo)dF=simple(dF)%簡(jiǎn)化符號(hào)表達(dá)式%得F'(x)=r-2*r*x/N-E%求F'(x0)并簡(jiǎn)化dFx0=subs(dF,x,x0);%將x=x0代入符號(hào)表達(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點(diǎn)穩(wěn)定，x1點(diǎn)不穩(wěn)定(根據(jù)平衡點(diǎn)穩(wěn)定性的準(zhǔn)則)；%若E>r,則結(jié)果正好相反。%在漁場(chǎng)魚量穩(wěn)定在x0的前提

7、下(E<r),求E使持續(xù)產(chǎn)量h(x0M到最大hm%通過(guò)分析(見教材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é)論是：%將捕撈率控制在固有增長(zhǎng)率的一半(E=r/2)時(shí)，能夠獲得最大的持續(xù)產(chǎn)量。2.(驗(yàn)證、編程)種群的相互競(jìng)爭(zhēng)P22222

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

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

10、%P(3,:)=P3(N1*(-1+k1)/(-1+k2*k1),N2*(-1+k2)/(-1+k2*k1)%P(4,:)=P4(0,0)%平衡點(diǎn)位于第一象限才有意義，故要求P3:k1,k2同小于1,或同大于1。%判斷平衡點(diǎn)的穩(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(-(fx

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

12、2)=r2*x2*(1-k2*x1/N1-x2/N2)%求方程的平衡點(diǎn)，即解代數(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個(gè)平衡點(diǎn)：%P(1,:)=P1(N1,0)%P(2,:)=P2(0,N2)%P(3,:)=P3(N

13、1*(-1+k1)/(-1+k2*k1),N2*(-1+k2)/(-1+k2*k1)%P(4,:)=P4(0,0)%平衡點(diǎn)位于第一象限才有意義，故要求P3:k1,k2同小于1,或同大于1%判斷平衡點(diǎn)的穩(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+gx2),x1,x2,

14、P(:,1),P(:,2);%替換p=simple(p)%簡(jiǎn)化符號(hào)表達(dá)式pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp('');Ppq%顯示結(jié)果%得到教材p225表1的前3列，經(jīng)測(cè)算可得該表的第4列，即穩(wěn)定條件2.2(驗(yàn)證、編程)計(jì)算與驗(yàn)證p227微分方程組xi(t)=r1x1(1-X2NiN2X2(t)=許(1-二X1X22N1N201二0.5,二1.6,1=2.5,R=1.8,M=1.6,N2=1(1)(驗(yàn)證)當(dāng)x1(0)=x2(0)=0.1時(shí)，求微分方程的數(shù)值解，將解的數(shù)值分別畫出x1(t)和x2(t)的曲線，它們同

15、在一個(gè)圖形窗口中。程序：%命令文件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)對(duì)應(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*x(1)*(1-x(1)/N1-k1*x(2

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

17、e45cp227',ts,x0);plot(x(:,1),x(:,2);%x(:1)是x1(t)的組函數(shù)值，x(:,2)對(duì)應(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)定性分析P228229模型：xi(t)=1為(1

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

19、ve(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)%簡(jiǎn)化#pq=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);Ppq%顯示結(jié)果j匚4|bbj1ndif!PX|匕I.fi

20、ai/噸"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.(驗(yàn)證)食餌-捕食者模型p230232模型的方程：X(t)=rx-axyx(0)=%y(t)=-dybxyy(0)=Vo要求：設(shè)r=1,d=0.5,a=0.1,b=0.02,xo=25,yo=2。輸入p231的程序并運(yùn)行，結(jié)果與

21、教材p232的圖1和圖2比較。給出2個(gè)M文件(見231)和程序運(yùn)行后輸出的圖形(比較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)'),%運(yùn)行中在圖上標(biāo)注pause,plot(x(:,1),x(:,2),grid,x(t

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

23、察收斂與否。(結(jié)果與教材p238239表1比較)下面是計(jì)算程序，在兩處下劃線的位置填入滿足要求的內(nèi)容。%不同b值下方程x(k+1)=b*x(k)*(1-x(k),k=0,1,2,.的計(jì)算結(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;%初值;%此處填入一條正確語(yǔ)句fork=1:99;%此處填入一條正確語(yǔ)句endK=(81:100)'%將取81100行disp(num2str(NaN,b;K,x(K,:),4);%取4位有效數(shù)字，

24、NaN表示不是數(shù)值(1)對(duì)程序正確填空，然后運(yùn)行。填入的正確語(yǔ)句和程序的運(yùn)行結(jié)果(對(duì)應(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)Ol4

25、4750.G40St).5死鴕Cl41ISD.61540.82360,8530.83170.羽，S370.11130,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.陰M

26、0.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)運(yùn)行以下程序，觀察顯示的圖形，與題(1)的數(shù)據(jù)對(duì)照，注意收斂的倍周期。%圖解法，圖1和圖2%文件名：p238fig_1_2.mclear;clc;closeall

27、;f=(x,b)b.*x.*(1-x);%定義f是函數(shù)的句柄，函數(shù)b*x*(1-x)含兩個(gè)變量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;%畫迭代軌跡線for

28、k=2:5x1=x(k,n);y1=x(k,n);plot(x0+i*y0,x0+i*y1,x1+i*y1,'r');%實(shí)部為橫坐標(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)

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

30、（b=3.45）圖解法：開始迭代的圖形（b=3.55）23倍圖解法：最后迭代的圖形(b=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);%初

31、值x0=0.2fork=1:100x(k+1,:)=b.*x(k,:).*(1-x(k,:);endplot(b,x(82:101,:),'.r','markersize',5);%可改變值5,調(diào)整數(shù)據(jù)點(diǎn)圖標(biāo)的大小xlabel('b');ylabel('x(k)');gridon要求:運(yùn)行以上程序。在運(yùn)行結(jié)果的圖形中，從對(duì)應(yīng)的b值上的“點(diǎn)數(shù)”判斷倍周期收斂。提示:放大圖形。程序的運(yùn)行結(jié)果（見238表1）:5.3收斂、分岔和混沌p240242畫出教材p241圖4模型的收斂、分岔和混沌程序的m文件如下:%圖4模型(6)的收斂、分岔和混

32、沌%文件名：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);%可改變值5,調(diào)整數(shù)據(jù)點(diǎn)圖標(biāo)的大小title('圖4模型的收斂、分岔和混沌,);xlabel(，參數(shù)b,);ylabel(,x(k)(k足夠大)，)；gridon;運(yùn)行以上m文件。

33、運(yùn)行結(jié)果(比較241圖4):fjehebiL'r|x1=1.XdaiVaw.ZxixhtIT0:IvD-aiJirtcrp工ndnr世1竄?白N、工"？*公口口國(guó)I國(guó)5.42n倍周期圖形p240242要求：在上題的程序中，修改b值，使運(yùn)行后顯示以下圖形:(1)單周期(1<b<3):>FLgUZEt匚巴摩EileEditVie*lasertloolsIesJrtottindo*HelpLfdA玲K1忑為立1X二總0"圉日模型的收被、分岔和混他U>；111r；：!:；,.,1',_：-：_r：；：:，Q.D-11Jitlllll；!II!

34、I1IBIIIIIIB|Je|IIIIIIJI(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|XEileIditInsertIoo>ls取drt叮tindwrHelp日白U)珞-：A-QF®配3ns"0圉日模型的

35、收斂、分岔和混沌(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=倍周期收斂的b的下限)c=bn_1;d=3.57;%2倍周期時(shí)收斂b>bn_1,b<3.57y=zeros(1,2*2An);%

36、行向量，用于存放xk最后的2X2n個(gè)值whiled-c>1e-12%使區(qū)間(c,d)足夠小b=(d+c)/2;x=0.2;%b取區(qū)間的中點(diǎn)fori=1:100000x=b*x*(1-x);endy(1)=x;%取最后2>2八門個(gè)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的上限近

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

38、um2str(-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倍周期時(shí)收斂b>bn_1,b<3.57y=zeros(1,2*2An);%行向量，用于存放xk最后的2X2n個(gè)值whiled-c>1e-12%使區(qū)間(c,d)足夠小b=(d+c)/2;x=0.2;%b取區(qū)間的中點(diǎn)fori=1:100000x=b*x*(1-x);endy(1)=x;%取最后2>2八門個(gè)xkfork=1:2A

39、(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.44935596723.54402147133.56435787943.56873999353.569679465A»|注：vpa的用法附1:

40、實(shí)驗(yàn)提示第1題提示：符號(hào)簡(jiǎn)化函數(shù)simple,變量替換函數(shù)sub符號(hào)簡(jiǎn)化函數(shù)simple的格式：simple(S)對(duì)符號(hào)表達(dá)式S嘗試多種不同的算法簡(jiǎn)化，以顯示S表達(dá)式的長(zhǎng)度最短的簡(jiǎn)化形式變量替換函數(shù)sub的格式:subs(S,OLD,NEW)將符號(hào)表達(dá)式S中的OLD變量替換為NEW變量附2:第7章穩(wěn)定性模型2157.1捕魚業(yè)的持續(xù)收獲7第/邕穩(wěn)定性模型雖然動(dòng)態(tài)過(guò)程的變化規(guī)律一般要用微分方程建立的動(dòng)態(tài)模型來(lái)描述，但是對(duì)于某些實(shí)際問(wèn)題，建模的主要目的并不是要尋求動(dòng)態(tài)過(guò)程每個(gè)瞬時(shí)的性態(tài)，而是研究某種意義下穩(wěn)定狀態(tài)的特征，特別是當(dāng)時(shí)間充分長(zhǎng)以后動(dòng)態(tài)過(guò)程的變化趨勢(shì).比如，在什么條件下描述過(guò)程的變量會(huì)越來(lái)

41、越接近某些確定的數(shù)值，在什么情況下又會(huì)越來(lái)越遠(yuǎn)離這些數(shù)值而導(dǎo)致過(guò)程不穩(wěn)定.為了分析這種穩(wěn)定與不穩(wěn)定的規(guī)律，常常不需要求解微分方程（并且我們將會(huì)看到，即使對(duì)于不太復(fù)雜的方程，解析解也不是總能得到的），而可以利用微分方程穩(wěn)定性理論，直接研究平衡狀態(tài)的穩(wěn)定性就行了.在下面建模過(guò)程中，大多數(shù)情況下我們將直接引用微分方程穩(wěn)定性理論的結(jié)果.不熟悉這方面內(nèi)容的讀者可以參閱17節(jié).差分方程的穩(wěn)定性與微分方程有很多相似之處.也有明顯的差別，工5和工7節(jié)將介紹這方面的內(nèi)容.7.1捕魚業(yè)的持續(xù)收獲可持續(xù)發(fā)展是一項(xiàng)基本國(guó)策，對(duì)于像漁業(yè)、林業(yè)這樣的再生資源，一定要注意適度開發(fā)，不能為了一時(shí)的高產(chǎn)去“竭澤而漁”，應(yīng)該在持

42、續(xù)穩(wěn)產(chǎn)的前提下追求產(chǎn)量或效益的最優(yōu)化.考察一個(gè)漁場(chǎng)，其中的魚量在天然環(huán)境下按一定規(guī)律增及，如果捕撈量恰好等于增長(zhǎng)最，那么漁場(chǎng)魚量將保持不變，這個(gè)捕撈量就可以持續(xù)下去本節(jié)要建立在捕播情況下漁場(chǎng)魚鼠遵從的方程，分析魚址穩(wěn)定的條件，并且在穩(wěn)定的前提下討論如何控制捕撈使持續(xù)產(chǎn)量或經(jīng)濟(jì)效益達(dá)到最大，最后研究所謂捕撈過(guò)度的問(wèn)題產(chǎn)模型記時(shí)刻t漁場(chǎng)中魚量為X。，關(guān)于的自然增長(zhǎng)和人工捕撈作如下假設(shè)：L在無(wú)捕撈條件下的增氏服從1。機(jī)比規(guī)律（見5.6節(jié)的阻滯增長(zhǎng)模型），即*(0=f(x)=rx|-N/(1)是固有增長(zhǎng)率，川是環(huán)境容許的最大魚最，用人工）表示單位時(shí)間的增長(zhǎng)量.2.單位時(shí)間的捕撈量（即產(chǎn)量）與漁場(chǎng)魚量工

43、成正比，比例常數(shù)E表示單位時(shí)間捕撈率,又稱為捕撈強(qiáng)度，可以用比如捕魚網(wǎng)眼的大小或出海漁船數(shù)量來(lái)控制其大小.于是單位時(shí)間的捕撈量為(2)h（x）-Ex根據(jù)以上假設(shè)并記尸（#）=/（）-A（x）得到捕措情況下漁場(chǎng)魚量滿足的方程x（t）=F（富）=1_木）-Ex我們并不需要解方程（3）以得到工（。的動(dòng)態(tài)變化過(guò)程.同希望知道漁場(chǎng)的穩(wěn)定魚量和保持穩(wěn)定的條件，即時(shí)間上足夠長(zhǎng)以后漁場(chǎng)魚最工。）的趨向，并由此確定最大持續(xù)產(chǎn)量.為此可以直接求方程（3）的平衡點(diǎn)并分析其穩(wěn)定性.令Ex-0得到兩個(gè)平衡點(diǎn)2=N(1-g,J?!=0不難算出F'(不)=E-r.k(盯)-r-E所以若£<r（5）有

44、尸（如）<0,尸（盯）>0,故通點(diǎn)穩(wěn)定百點(diǎn)不穩(wěn)定（判斷平衡點(diǎn)穩(wěn)定性的準(zhǔn)則見7.7節(jié)）彳若E>r,則結(jié)果正好相反.y-f(x)圖】最大持續(xù)產(chǎn)拼的圖解法xf=M2xDNE是捕撈率，r是最大的增長(zhǎng)率，上述分析表明，只要捕撈適度（£<1）.就可使?jié)O場(chǎng)魚量穩(wěn)定在質(zhì)，從而獲得持續(xù)產(chǎn)最M不）二取°而當(dāng)捕撈過(guò)度時(shí)（E>力，漁場(chǎng)魚量將趨向盯=0,當(dāng)然談不上獲得持續(xù)產(chǎn)量了.進(jìn)一步討論漁場(chǎng)魚量穩(wěn)定在與的前提下，如何控制捕撈強(qiáng)度E使持續(xù)產(chǎn)量最大的問(wèn)甌用圖解法可以非常簡(jiǎn)單地得到216；結(jié)果.根據(jù)（1）,（2）式作拋物線耳）和直線=乂*）=&,如圖1.注意到y(tǒng)=/

45、«）在原點(diǎn)的切線為？二%所以在條件（5）下,二株必與有交點(diǎn)P,P的橫坐標(biāo)就是穩(wěn)定平衡點(diǎn)小.根據(jù)假設(shè)2,尸點(diǎn)的縱坐標(biāo)為為穩(wěn)定條件下單位時(shí)間的持續(xù)產(chǎn)量.由圖1立刻知道，當(dāng)y=&與在拋物線頂點(diǎn)P*相交時(shí)可獲得最大的持續(xù)產(chǎn)量,此時(shí)的穩(wěn)定平衡點(diǎn)為(6)且單位時(shí)間的最大持續(xù)產(chǎn)量為而由（4）式不難算出保持漁場(chǎng)魚量穩(wěn)定在工；的捕撈率為E=2綜上所述，產(chǎn)量模型的結(jié)論是將捕撈率控制在固有增長(zhǎng)率r的一半，更簡(jiǎn)單一些，可以說(shuō)使?jié)O場(chǎng)魚量保持在最大魚量N的一半時(shí)，能夠獲得最大的持續(xù)產(chǎn)量.效益模型從經(jīng)濟(jì)角度看不應(yīng)追求產(chǎn)量最大，而應(yīng)考慮效益最佳.如果經(jīng)濟(jì)效益用從捕撈所得的收入中扣除開支后的利潤(rùn)來(lái)衡量，并且簡(jiǎn)

46、單地假設(shè)：魚的銷售單價(jià)為常數(shù)P,單位捕撈率（如每條出海漁船）的費(fèi)用為常數(shù)c,那么單位時(shí)間的收入T和支出S分別為T-ph(x)=pEx,S=cE(9)單位時(shí)間的利潤(rùn)為RT-S-pEx-cE在穩(wěn)定條件*=褊下，以（4）代入（10）式，得(10)/?(£)=T(E)-S(£)=pNE(T-(11)用微分法容易求出使利潤(rùn)K（E）達(dá)到最大的捕撈強(qiáng)度為(12)將4代入（4）式，可得最大利潤(rùn)下的漁場(chǎng)穩(wěn)定魚量5及單位時(shí)間的持續(xù)產(chǎn)量K為(13)(14)217將（12）14）式與產(chǎn)量模型中的（6）（8）式相比較可以看出,在最大效益原則下捕撈率和持續(xù)產(chǎn)量均有所減少，而渣場(chǎng)應(yīng)保持的穩(wěn)定魚量有所增加

47、，并且減少或增加的部分隨著捕撈成本e的增長(zhǎng)而變大，循著銷售價(jià)格p的增長(zhǎng)而變小.請(qǐng)讀者分析這些結(jié)果是符合實(shí)際情況的.捕措過(guò)度上面的效益模型是以計(jì)劃捕撈（或稱封閉式捕撈）為基礎(chǔ)的，即漁場(chǎng)由單獨(dú)的經(jīng)營(yíng)者有計(jì)劃地捕撈，可以追求最大利潤(rùn).如果漁場(chǎng)向眾多盲目的經(jīng)營(yíng)者開放，比如在公海上無(wú)規(guī)則地捕撈,那么即使只有微薄的利潤(rùn)，經(jīng)營(yíng)者也會(huì)蟒擁而去，這種情況稱為盲目捕撈（或開放式捕撈）.這種捕撈方式將導(dǎo)致捕撈過(guò)度，下面討論這個(gè)模型.（11）式給出了利潤(rùn)與捕撈強(qiáng)度的關(guān)系令R（E）=0的解為跖，可得(15)當(dāng)E風(fēng)時(shí)，利潤(rùn)靛（E）0,盲目的經(jīng)營(yíng)者們會(huì)加大捕措強(qiáng)度；若EaE-利潤(rùn)R（E）<0*他們當(dāng)然要減小強(qiáng)度,所以

48、乙是盲目捕撈下的臨界強(qiáng)度M也可由圖解法確定.在圖2中以£為橫坐標(biāo)，按（11）式畫出T（E）和S（E）,它們交點(diǎn)的橫坐標(biāo)即為名（圖2中的號(hào)1或后加）,由（15）式或圖2容易知道，酌存在的必要條件（即4>。）是圖2盲目捕撈強(qiáng)度的圖解法即售價(jià)大于（相對(duì)于總量而言）成本.并且由（15）式可知，成本越低，售價(jià)越高，則瓦越大,將（15）代人（4）式，得到盲目捕撈下的漁場(chǎng)穩(wěn)定魚量為（17）%完全由成本-價(jià)格比決定，隨著價(jià)格的上升和成本的下降%將迅速減少，出現(xiàn)捕撈過(guò)度.比校32）和（15）式可知=2號(hào),即盲目捕撈強(qiáng)度比最大效益下捕撈噩度大一倍.從（15）式和圖2還可以得到，當(dāng)竟<p<

49、;2時(shí)，（昂<）/</，如圖2中218凰，稱經(jīng)濟(jì)學(xué)捕撈過(guò)度；當(dāng)P>2右時(shí)，曷>£，如圖2中”，稱生態(tài)學(xué)捕撈過(guò)度.:*本節(jié)完*評(píng)注為了研究漁業(yè)的產(chǎn)量、效益及捕撈過(guò)度問(wèn)題，首先在對(duì)魚的自然增長(zhǎng)和捕撈情況的合理假設(shè)F,建立漁場(chǎng)魚量的基本方程（3）,并利用平衡點(diǎn)穩(wěn)定性分析確定了保持漁場(chǎng)魚鼠楫定的條件.產(chǎn)量、效益和捕撈過(guò)度3個(gè)模型在穩(wěn)定的前提下步步深入，數(shù)學(xué)推導(dǎo)過(guò)程十分簡(jiǎn)單，卻得到了在定性關(guān)系上與實(shí)際情況完全符合的結(jié)果.如果改變對(duì)魚的自然增長(zhǎng)和人工捕撈的假設(shè)，模型及結(jié)果將隨之變化（習(xí)題1,2）.72軍備競(jìng)賽兩個(gè)國(guó)家或兩個(gè)國(guó)家集團(tuán)之間由于相互不信任和各種矛盾的存在、發(fā)展而

50、不斷增加自己的軍事力量，防御對(duì)方可能發(fā)動(dòng)的戰(zhàn)爭(zhēng).能否用一個(gè)數(shù)學(xué)模型描述這種軍備竟賽的過(guò)程，從定性和定髭的角度對(duì)競(jìng)賽的結(jié)果作出解釋或預(yù)測(cè).本節(jié)介紹LF.Richardsgl939年提出的一個(gè)模型.當(dāng)然影響軍備競(jìng)賽的因素是錯(cuò)綜復(fù)雜的，無(wú)法用數(shù)學(xué)工具給以恰當(dāng)?shù)膱A滿的描述.這個(gè)模型只不過(guò)告訴我們，一個(gè)復(fù)雜的實(shí)際過(guò)程可以被合理地簡(jiǎn)化到什么程度，得到的結(jié)果又怎樣用來(lái)解釋實(shí)際現(xiàn)象，并且如果允許軍備比賽沿著機(jī)械的本能和固定的模式發(fā)展的話.會(huì)導(dǎo)致什么樣的結(jié)局廊.模型假設(shè)與構(gòu)成為了方便起見，用軍備這個(gè)詞表示軍事力量的總和，如兵力、裝備、軍事頸算等.甲乙雙方在時(shí)刻t的軍備分別記作工3）和江。，假定它們的變化只取決于

51、下面3個(gè)因素工1,由于相互不信任及矛盾的發(fā)展，一方軍備越大，另一方軍備增加得越快.2,由于各方本身經(jīng)濟(jì)實(shí)力的限制.任一方軍備越大，對(duì)軍備增長(zhǎng)的制約作用越大.3由于相互敵視或領(lǐng)土爭(zhēng)端，每一方都存在膂增加軍備的固有潛力.進(jìn)一步假定前兩個(gè)因素的影響是線性的，第3個(gè)因素的影響是常數(shù)，那么雙上）和貝。的變化過(guò)程可用微分方程組：£（£）=-5+Gy+g（I）y（0三自=By+%表示，其中的系數(shù)均大于或等于0.是對(duì)方軍備剌激程度的度最e,©是己方經(jīng)濟(jì)實(shí)力制約程度的度量灑力是巴方軍備競(jìng)賽的固有潛力-如果我們感興趣的是軍備競(jìng)賽的結(jié)局由什么因素決定，而不關(guān)心競(jìng)賽的過(guò)程.那么只需用微分

52、方程穩(wěn)定性理論討論時(shí)間充分長(zhǎng)以后工（。，¥”）的變化趨勢(shì)，即方程（1）的平衡點(diǎn)的穩(wěn)定情況令（1）式右端等于0,容易算出平衡點(diǎn)（與，九）為fl印-H'九一4一奴方程（1）的系數(shù)矩陣為F-a*1A=>，p于是按照判斷平衡點(diǎn)穩(wěn)定性的方法計(jì)算（見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）時(shí)，平衡點(diǎn)（與，為）是穩(wěn)定的：反之.是不穩(wěn)定的.這就是說(shuō)，在（5）式的條件下，時(shí)間足夠長(zhǎng)以后雙方的軍備將分別趨向一個(gè)有限值，軍備競(jìng)賽是穩(wěn)定的、模型的定性解釋根據(jù)方程（1）和平衡點(diǎn)穩(wěn)定性的分析，

53、可以解釋幾個(gè)簡(jiǎn)單而又重要的現(xiàn)象.1.條件（5）表明，當(dāng)雙方的經(jīng)濟(jì)制約程度邛大于雙方的軍備刺激程度kl時(shí)，軍備競(jìng)賽才會(huì)趨向穩(wěn)定.反之聲U）4（f）將趨向無(wú)窮，競(jìng)賽無(wú)限地進(jìn)行下去，可能導(dǎo)致戰(zhàn)爭(zhēng).2,由（2）式，如果g=A=0,則0=0,兀=0是方程（1）的平衡點(diǎn)，并且在條件（5）下它是穩(wěn)定的，于是如果在某個(gè)時(shí)候品有=y（/0）=。，%了就永遠(yuǎn)保持為E這種情況可以解釋為雙方不存在任何敵視和爭(zhēng)端，通過(guò)裁軍可以達(dá)到持久和平.兩個(gè)友好的鄰國(guó)正是這樣.3 .如果即使由于某種原因（如裁軍協(xié)定）在某個(gè)時(shí)候雙方軍備大減，不妨設(shè)M%）=y（r0）=。，那么因?yàn)閤=h,也將使雙方重整軍備.這說(shuō)明未經(jīng)和解的裁軍（即不消

54、除敵視或領(lǐng)土爭(zhēng)端）是不會(huì)持久的.4 .如果由于某種原因（如戰(zhàn)敗或協(xié)議），在某個(gè)時(shí)候一方的軍備大減，不妨設(shè)立色）=0,那么因?yàn)槲宥?心也將使該方重整軍備.這說(shuō)明存在不信任（無(wú)*0）或固有爭(zhēng)端3*0）的單方面裁軍也不會(huì)持久.模型參數(shù)的估計(jì)為了利用（5）式判斷軍備競(jìng)賽是否會(huì)趨于穩(wěn)定，需要知道白班,AN的數(shù)直估計(jì)這些參數(shù)無(wú)疑是很困難的，下面是RichardQi提出的一種方法.1. 的估計(jì)設(shè)工（0）=0,當(dāng)二較小時(shí)，忽略g和-ox的作用,并近似地假定丁二為不變，由方程（1）得i二%220如果當(dāng)t=丁時(shí)，則由（6）式得到（7）這說(shuō)明A"是甲方軍備從0到趕上乙方軍備外所需的時(shí)間.例如，譙國(guó)從19

55、33年開始重整軍備，只用了約3年的時(shí)間就趕上了它的鄰國(guó).假設(shè)它增加軍備的固有潛力g被制約效應(yīng)ax所抵消，那么可以認(rèn)為德國(guó)的A-=3年，即上=0.3./可以類似地估計(jì)，或者合理地假定它與國(guó)家的經(jīng)濟(jì)實(shí)力成正比.這樣若乙國(guó)的經(jīng)濟(jì)實(shí)力是德國(guó)的2倍，則可以估計(jì)20.6.嘉a,尸的估計(jì)設(shè)g=0,y=0,由方程（】）可得x（t）=.（0”7以代人，算出這表示儀7是在乙方無(wú)軍備時(shí)甲方軍備減少到原來(lái)的,所需的時(shí)間.Richardson認(rèn)e為這大概是一個(gè)國(guó)家議會(huì)的任期，對(duì)于議會(huì)任期5年的國(guó)家來(lái)說(shuō)©=0.2.對(duì)模型和叁數(shù)的粗晡檢覽考察第一次世界大戰(zhàn)前夕，歐洲的兩個(gè)國(guó)家同盟法俄同盟和德奧匈同盟的軍備競(jìng)賽情況.-兩個(gè)同盟的經(jīng)濟(jì)實(shí)力大致相等,且約為德國(guó)的3倍，因?yàn)榈聡?guó)的A=0.3,所以這兩個(gè)同盟的左=/