數(shù)學(xué)建模實(shí)驗(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(t)為t時(shí)刻漁場中的魚量。r是固有增長率。N是環(huán)境容許的最大魚量。E是捕撈強(qiáng)度，即單位時(shí)間捕撈率。要求：運(yùn)行下面的m文件，并把相應(yīng)結(jié)果填空，即填入“_”。%7.1 捕魚業(yè)的持續(xù)收獲產(chǎn)量模型%文件名：p215_217.mclear; clc;%無捕撈條件下單位時(shí)間的增長量：f(x)=rx(1-x/N)%捕撈條件下單位時(shí)間的捕撈量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕撈情況下漁場魚量滿足的方程：x'(t)=F(x)%滿足F(x)=

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

3、%求F'(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(xiàn)'(x1)>0，故x0點(diǎn)穩(wěn)定，x1點(diǎn)不穩(wěn)定（根據(jù)平衡點(diǎn)穩(wěn)定性的準(zhǔn)則）；%若 E>r，則結(jié)果正好相反。%在漁場魚量穩(wěn)定在x0的前提下（E<r），求E使持續(xù)產(chǎn)量h(x0)達(dá)到最大hm。%通過分析（見教材p216圖1），只需求x0*使f(x)達(dá)到最大，且hm=f

4、(x0*)。syms r x Nfx=r*x*(1-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）時(shí)，能夠獲得最大的持續(xù)產(chǎn)量。符號簡化函數(shù)simple，變量替換函數(shù)sub的用法見提示。 給出填空后的M文件（見215217）：%7.1 捕魚業(yè)的持續(xù)收獲產(chǎn)量模型%文件名：p215_217.mclear; clc;

5、%無捕撈條件下單位時(shí)間的增長量：f(x)=rx(1-x/N)%捕撈條件下單位時(shí)間的捕撈量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕撈情況下漁場魚量滿足的方程：x'(t)=F(x)%滿足F(x)=0的點(diǎn)x為方程的平衡點(diǎn)%求方程的平衡點(diǎn)syms r x N E; %定義符號變量Fx=r*x*(1-x/N)-E*x; %創(chuàng)建符號表達(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);%符號變量x的結(jié)構(gòu)類型成為<2×

6、;1sym>%求F(x)的微分F'(x)syms x; %定義符號變量x的結(jié)構(gòu)類型為<1×1sym>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'

7、;(x0)<0，F(xiàn)'(x1)>0，故x0點(diǎn)穩(wěn)定，x1點(diǎn)不穩(wěn)定（根據(jù)平衡點(diǎn)穩(wěn)定性的準(zhǔn)則）；%若 E>r，則結(jié)果正好相反。%在漁場魚量穩(wěn)定在x0的前提下（E<r），求E使持續(xù)產(chǎn)量h(x0)達(dá)到最大hm。%通過分析（見教材p216圖1），只需求x0*使f(x)達(dá)到最大，且hm=f(x0*)。syms r x Nfx=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

8、*=N(1-E/r)，可得 E*= r/2 ，見P217(8)式%產(chǎn)量模型的結(jié)論是：%將捕撈率控制在固有增長率的一半（E=r/2）時(shí)，能夠獲得最大的持續(xù)產(chǎn)量。2.（驗(yàn)證、編程）種群的相互競爭P222228模型：其中，x1(t), x2(t)分別是甲乙兩個(gè)種群的數(shù)量。r1, r2是它們的固有增長率。N1, N2是它們的最大容量。1：單位數(shù)量乙（相對N2）消耗的供養(yǎng)甲的食物量為單位數(shù)量甲（相對N1）消耗的供養(yǎng)甲的食物量的1倍。對2可作相應(yīng)解釋。2.1（編程）穩(wěn)定性分析p224225要求：補(bǔ)充如下指出的程序段，然后運(yùn)行該m文件，對照教材上的相應(yīng)結(jié)果。%7.3 種群的相互競爭穩(wěn)定性分析%文件名：p22

9、4_225.mclear; clc;%甲乙兩個(gè)種群滿足的增長方程：% x1'(t)=f(x1,x2)=r1*x1*(1-x1/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）使用符號表達(dá)式；（2）用函數(shù)solve求解代數(shù)方程組，解放入x1, x2；（3）調(diào)整解（平衡點(diǎn)）的順序放入P中（見下面注釋所示），P的結(jié)構(gòu)類型為<4×2sym>，P的第1列對應(yīng)x1，第2列對應(yīng)x2。

10、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（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)式。syms x1 x2; %重新定義fx1=diff(f,'x1'); fx2=diff(f,'x2');gx1=diff

11、(g,'x1'); gx2=diff(g,'x2');disp(' '); A=fx1,fx2;gx1,gx2 %顯示結(jié)果p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2); %替換p=simple(p);%簡化符號表達(dá)式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp(' '); P p q %顯示結(jié)果%得到教材p225表1的前3列，經(jīng)測算可得該表的第4列，即穩(wěn)定條件。 補(bǔ)充后的程序和運(yùn)行結(jié)果（見225表1）：2341%7.3 種群的相互競爭穩(wěn)定性分

12、析%文件名：p224_225.mclear; clc;%甲乙兩個(gè)種群滿足的增長方程：% x1'(t)=f(x1,x2)=r1*x1*(1-x1/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%編寫出該程序段。syms x1 x2 r1 r2 N1 N2 k1 k2;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)

13、;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（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)式。syms x1 x2; %重新定義fx1=diff(f,'x1'); fx2=di

14、ff(f,'x2');gx1=diff(g,'x1'); gx2=diff(g,'x2');disp(' '); A=fx1,fx2;gx1,gx2 %顯示結(jié)果p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2); %替換p=simple(p);%簡化符號表達(dá)式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp(' '); P p q %顯示結(jié)果%得到教材p225表1的前3列，經(jīng)測算可得該表的第4列，即穩(wěn)定條件。2.2（驗(yàn)證、編程）計(jì)算與驗(yàn)

15、證p227微分方程組(1)（驗(yàn)證）當(dāng)x1(0)=x2(0)=0.1時(shí)，求微分方程的數(shù)值解，將解的數(shù)值分別畫出x1(t)和x2(t)的曲線，它們同在一個(gè)圖形窗口中。程序：%命令文件ts=0:0.2:8;x0=0.1,0.1;t,x=ode45('p227',ts,x0);plot(t,x(:,1),t,x(:,2);%x(:1)是x1(t)的一組函數(shù)值，x(:,2)對應(yīng)x2(t)grid on;axis(0 8 0 2);text(2.4,1.55,'x1(t)');text(2.2,0.25,'x2(t)');%函數(shù)文件名：p227.mfunct

16、ion y=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)/N2), r2*x(2)*(1-k2*x(1)/N1-x(2)/N2)'(1) 運(yùn)行程序并給出結(jié)果（比較227圖2(a)）：(2) （驗(yàn)證）將x1(0)=1, x2(0)=2代入(1)中的程序并運(yùn)行。給出結(jié)果（比較227圖2(b)）：(3)（編程）在同一圖形窗口內(nèi)，畫(1)和(2)的相軌線，相軌線是以x1(t)為橫坐標(biāo)，x2(t)為縱坐標(biāo)所得到的一條曲線。具體要求參照下圖。圖中的兩條“點(diǎn)線”直線，一條的兩個(gè)端點(diǎn)

17、為(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=ode45('p227',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 )');hold on;x0=1,2;t,x=ode45('p227',ts,x0);plot(x(:,1),x(:,2);text(1.02,1.9,'(

18、1, 2 )');plot(0,1,1,0,':r',0,1.6,2,0,':r');%畫兩條直線grid on;3.（編程）種群的相互依存穩(wěn)定性分析P228229模型：其中，x1(t), x2(t)分別是甲乙兩個(gè)種群的數(shù)量。r1, r2是它們的固有增長率。N1, N2是它們的最大容量。1：單位數(shù)量乙（相對N2）提供的供養(yǎng)甲的食物量為單位數(shù)量甲（相對N1）消耗的供養(yǎng)甲的食物量的1倍。對2可作相應(yīng)解釋。要求： 修改題2.1的程序，求模型的平衡點(diǎn)及穩(wěn)定性。給出程序及其運(yùn)行結(jié)果（比較229表1，注：只要最終結(jié)果）：clear; clc;syms x1 x2 r

19、1 r2 N1 N2 k1 k2;f=r1*x1*(1-x1/N1+k1*x2/N2);g=r2*x2*(-1+k2*x1/N1-x2/N2);x1,x2=solve(f,g);P=x1(2,4,1,3),x2(2,4,1,3);syms x1 x2; %重新定義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=simp

20、le(p);%簡化符號表達(dá)式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);P p q %顯示結(jié)果4.（驗(yàn)證）食餌-捕食者模型p230232模型的方程：要求：設(shè)r =1，d =0.5，a=0.1，b=0.02，x0=25，y0=2。輸入p231的程序并運(yùn)行，結(jié)果與教材p232的圖1和圖2比較。 給出2個(gè)M文件（見231）和程序運(yùn)行后輸出的圖形（比較232圖1、2）：函數(shù)M文件：function xdot=shier(t,x)r=1; d=0.5; a=0.1  b=0.02 xdot=(r-a*x(2).*x(1); (-d+

21、b*x(1).*x(2);命令M文件：ts=0 :0.1 :15;x0=25, 2;t,x=ode45('shier',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), y(t)圖形：相軌線y(x)圖形：5.（驗(yàn)證）差分形式的阻滯增長模型p236242阻滯增長模型用微分方程描述為：也可用差分方程描述為：上式可簡化為一階非線性差分方程：考察給定b和x0值后，當(dāng)k 時(shí)，xk的收斂情況（實(shí)際

22、上取k足夠大就可以了）。5.1 數(shù)值解法和圖解法p238240(1) 取x0=0.2，分別取b = 1.7, 2.6, 3.3, 3.45, 3.55, 3.57，對方程計(jì)算出x1 x100的值，顯示x81 x100的值。觀察收斂與否。（結(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; clear all; format compact;b=1.7,2.6,3.3,3.45,3.55,3.57;x=zeros(100,l

23、ength(b);x0=0.2; %初值 ; %此處填入一條正確語句for k=1:99 ; %此處填入一條正確語句endK=(81:100); %將取81100行disp(num2str(NaN,b; K,x(K,:),4);%取4位有效數(shù)字，NaN表示不是數(shù)值(1) 對程序正確填空，然后運(yùn)行。填入的正確語句和程序的運(yùn)行結(jié)果（對應(yīng)不同的b值）見238表1：x(1,:)=b*x0*(1-x0);x(k+1,:)=b.*x(k,:).*(1-x(k,:);(2) 運(yùn)行以下程序，觀察顯示的圖形，與題(1)的數(shù)據(jù)對照，注意收斂的倍周期。%圖解法，圖1和圖2%文件名：p238fig_1_2.mclea

24、r; clc; close all;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;for k=1:100 x(k+1,:)=f(x(k,:),b);endfor n=1:length(b) figure(n);%指定圖形窗口figure n subplot(1,2,1);%開始迭代的圖形 fplot(x)f(x,b(n),x,0 1 0 1);%x是自變量，畫曲線y=bx(1-x)和直線y=x axis square;

25、hold on; x0=x(1,n); y0=0; %畫迭代軌跡線 for k=2:5 x1=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; end title('圖解法：開始迭代的圖形（b=' num2str(b(n) '）'); hold off; subplot(1,2,2); %最后迭代的圖形 fplot(x)f(x,b(n),x,0 1 0 1); axis square; hold on; xy(1:2:41)=x

26、(81:101,n)+i*x(81:101,n);%盡量不用循環(huán) xy(2:2:40)=x(81:100,n)+i*x(82:101,n); plot(xy,'r'); title('圖解法：最后迭代的圖形（b=' num2str(b(n) '）'); hold off;end(2) 運(yùn)行程序并給出結(jié)果（對應(yīng)不同的b值）見238圖1、2：20倍20倍21倍22倍23倍混沌5.2 b值下的收斂圖形p238240下面程序是在不同b值下的收斂圖形。%b值下的收斂圖形%文件名：p212tab1figure.m%方程（6） xk+1=b*xk*(1-xk)

27、, k=0,1,2,.clear; clc; close all;b=3.3 3.45 3.55 3.57;x=0.2*ones(size(b); %初值x0=0.2for k=1:100 x(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)');grid on要求： 運(yùn)行以上程序。 在運(yùn)行結(jié)果的圖形中，從對應(yīng)的b值上的“點(diǎn)數(shù)”判斷倍周期收斂。提示：放

28、大圖形。 程序的運(yùn)行結(jié)果（見238表1）：5.3 收斂、分岔和混沌p240242畫出教材p241圖4 模型的收斂、分岔和混沌。程序的m文件如下：%圖4 模型（6）的收斂、分岔和混沌%文件名：p241fig4.m%模型：xk+1=b*xk*(1-xk), k=0,1,2,.clear; clc; close all;b=2.5:0.0001:4; %參數(shù)b的間隔取值x=0.2*ones(size(b);xx=; n=100000;for k=1:n x=b.*x.*(1-x); if k>=n-50 xx=xx;x; endendplot(b,xx,'.b','ma

29、rkersize',5); %可改變值5，調(diào)整數(shù)據(jù)點(diǎn)圖標(biāo)的大小title('圖4 模型的收斂、分岔和混沌');xlabel('參數(shù) b');ylabel('x(k) （k足夠大）');grid on; 運(yùn)行以上m文件。運(yùn)行結(jié)果（比較241圖4）：5.4 2n倍周期圖形p240242要求：在上題的程序中，修改b值，使運(yùn)行后顯示以下圖形：(1) 單周期（1<b<3）：(2) 2倍周期（3<b<3.449）：(3) 22倍周期（3.449<b<3.544）：(4) 23倍周期（3.544<b<3.

30、564）：5.5（編程）求2n倍周期的b值區(qū)間p240241求2n倍周期收斂的b的上限的程序如下：function bmax=fun(bn_1,n)%求2n倍周期收斂的b的上限。%bn_1是2(n-1)倍周期收斂的b的上限（或2n倍周期收斂的b的下限）。c=bn_1; d=3.57;% 2n倍周期時(shí)收斂b>bn_1，b<3.57y=zeros(1,2*2n);%行向量，用于存放xk最后的2×2n個(gè)值while d-c>1e-12 %使區(qū)間(c, d)足夠小 b=(d+c)/2; x=0.2; %b取區(qū)間的中點(diǎn) for i=1:100000 x=b*x*(1-x);

31、end y(1)=x; %取最后2×2n個(gè)xk for k=1:2(n+1)-1 y(k+1)=b*y(k)*(1-y(k); end if norm(y(1:2n)-y(2n+1:2(n+1)<1e-12 %范數(shù)，比較 c=b;%滿足2n倍周期收斂的b給區(qū)間(c, d)的下限c else d=b; %不滿足2n倍周期收斂的b給區(qū)間(c, d)的上限d endendbmax=c; % 2n倍周期收斂的b的上限近似要求：編寫程序，調(diào)用以上函數(shù)文件，求單倍周期、2倍周期、25倍周期收斂的b的上限近似值，輸出要求有10位有效數(shù)字。運(yùn)行結(jié)果輸出形式如下：提示：可使用num2str函數(shù)。用下面的程序結(jié)構(gòu)，會使運(yùn)行速度加