数学建模实验答案_稳定性模型

1、实验08 稳定性模型（4学时）（第7章 稳定性模型）1.（验证）捕鱼业的持续收获 产量模型p215219产量模型：其中，x(t)为t时刻渔场中的鱼量。r是固有增长率。N是环境容许的最大鱼量。E是捕捞强度，即单位时间捕捞率。要求：运行下面的m文件，并把相应结果填空，即填入“_”。%7.1 捕鱼业的持续收获产量模型%文件名：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'(t)=F(x)%满足F(x)=

2、0的点x为方程的平衡点%求方程的平衡点syms r x N E; %定义符号变量Fx=r*x*(1-x/N)-E*x; %创建符号表达式x=solve(Fx,x) %求解F(x)=0（求根）%得到两个平衡点，记为：% x0= , x1= ，见P216(4)式x0=x(2);x1=x(1);%符号变量x的结构类型成为<2×1sym>%求F(x)的微分F'(x)syms x; %定义符号变量x的结构类型为<1×1sym>dF=diff(Fx,'x'); %求导dF=simple(dF) %简化符号表达式%得 F'(x)=

3、%求F'(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点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若 E>r，则结果正好相反。%在渔场鱼量稳定在x0的前提下（E<r），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且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)式%产量模型的结论是：%将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。符号简化函数simple，变量替换函数sub的用法见提示。 给出填空后的M文件（见215217）：%7.1 捕鱼业的持续收获产量模型%文件名：p215_217.mclear; clc;

5、%无捕捞条件下单位时间的增长量：f(x)=rx(1-x/N)%捕捞条件下单位时间的捕捞量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕捞情况下渔场鱼量满足的方程：x'(t)=F(x)%满足F(x)=0的点x为方程的平衡点%求方程的平衡点syms r x N E; %定义符号变量Fx=r*x*(1-x/N)-E*x; %创建符号表达式x=solve(Fx,x) %求解F(x)=0（求根）%得到两个平衡点，记为：% x0= -N*(-r+E)/r , x1= 0 ，见P216(4)式x0=x(2);x1=x(1);%符号变量x的结构类型成为<2×

6、;1sym>%求F(x)的微分F'(x)syms x; %定义符号变量x的结构类型为<1×1sym>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'

7、;(x0)<0，F'(x1)>0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若 E>r，则结果正好相反。%在渔场鱼量稳定在x0的前提下（E<r），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且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)式%产量模型的结论是：%将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。2.（验证、编程）种群的相互竞争P222228模型：其中，x1(t), x2(t)分别是甲乙两个种群的数量。r1, r2是它们的固有增长率。N1, N2是它们的最大容量。1：单位数量乙（相对N2）消耗的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的1倍。对2可作相应解释。2.1（编程）稳定性分析p224225要求：补充如下指出的程序段，然后运行该m文件，对照教材上的相应结果。%7.3 种群的相互竞争稳定性分析%文件名：p22

9、4_225.mclear; clc;%甲乙两个种群满足的增长方程：% 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)%求方程的平衡点，即解代数方程组（见P224的(5)式）% f(x1,x2)=0% g(x1,x2)=0编写出该程序段。提示（1）使用符号表达式；（2）用函数solve求解代数方程组，解放入x1, x2；（3）调整解（平衡点）的顺序放入P中（见下面注释所示），P的结构类型为<4×2sym>，P的第1列对应x1，第2列对应x2。

10、x1x2=x1,x2 %显示结果disp(' '); P%调整位置后的4个平衡点：% 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）%平衡点位于第一象限才有意义，故要求P3：k1, k2同小于1，或同大于1。%判断平衡点的稳定性参考教材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 %显示结果p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2); %替换p=simple(p);%简化符号表达式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp(' '); P p q %显示结果%得到教材p225表1的前3列，经测算可得该表的第4列，即稳定条件。 补充后的程序和运行结果（见225表1）：2341%7.3 种群的相互竞争稳定性分

12、析%文件名：p224_225.mclear; clc;%甲乙两个种群满足的增长方程：% 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)%求方程的平衡点，即解代数方程组（见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 %显示结果disp(' '); P%调整位置后的4个平衡点：% 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）%平衡点位于第一象限才有意义，故要求P3：k1, k2同小于1，或同大于1。%判断平衡点的稳定性参考教材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 %显示结果p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2); %替换p=simple(p);%简化符号表达式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp(' '); P p q %显示结果%得到教材p225表1的前3列，经测算可得该表的第4列，即稳定条件。2.2（验证、编程）计算与验

15、证p227微分方程组(1)（验证）当x1(0)=x2(0)=0.1时，求微分方程的数值解，将解的数值分别画出x1(t)和x2(t)的曲线，它们同在一个图形窗口中。程序：%命令文件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)的一组函数值，x(:,2)对应x2(t)grid on;axis(0 8 0 2);text(2.4,1.55,'x1(t)');text(2.2,0.25,'x2(t)');%函数文件名：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) 运行程序并给出结果（比较227图2(a)）：(2) （验证）将x1(0)=1, x2(0)=2代入(1)中的程序并运行。给出结果（比较227图2(b)）：(3)（编程）在同一图形窗口内，画(1)和(2)的相轨线，相轨线是以x1(t)为横坐标，x2(t)为纵坐标所得到的一条曲线。具体要求参照下图。图中的两条“点线”直线，一条的两个端点

17、为(0,1)和(1,0)，另一条的两个端点为(0, 2)和(1.6, 0)。(3) 给出程序及其运行结果（比较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)的一组函数值，x(:,2)对应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.（编程）种群的相互依存稳定性分析P228229模型：其中，x1(t), x2(t)分别是甲乙两个种群的数量。r1, r2是它们的固有增长率。N1, N2是它们的最大容量。1：单位数量乙（相对N2）提供的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的1倍。对2可作相应解释。要求： 修改题2.1的程序，求模型的平衡点及稳定性。给出程序及其运行结果（比较229表1，注：只要最终结果）：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);%简化符号表达式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);P p q %显示结果4.（验证）食饵-捕食者模型p230232模型的方程：要求：设r =1，d =0.5，a=0.1，b=0.02，x0=25，y0=2。输入p231的程序并运行，结果与教材p232的图1和图2比较。 给出2个M文件（见231）和程序运行后输出的图形（比较232图1、2）：函数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)'), %运行中在图上标注pause, plot(x(:,1),x(:,2), grid, x(t), y(t)图形：相轨线y(x)图形：5.（验证）差分形式的阻滞增长模型p236242阻滞增长模型用微分方程描述为：也可用差分方程描述为：上式可简化为一阶非线性差分方程：考察给定b和x0值后，当k 时，xk的收敛情况（实际

22、上取k足够大就可以了）。5.1 数值解法和图解法p238240(1) 取x0=0.2，分别取b = 1.7, 2.6, 3.3, 3.45, 3.55, 3.57，对方程计算出x1 x100的值，显示x81 x100的值。观察收敛与否。（结果与教材p238239表1比较）下面是计算程序，在两处下划线的位置填入满足要求的内容。%不同b值下方程 x(k+1)=b*x(k)*(1-x(k), k=0,1,2,.的计算结果%文件名：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位有效数字，NaN表示不是数值(1) 对程序正确填空，然后运行。填入的正确语句和程序的运行结果（对应不同的b值）见238表1：x(1,:)=b*x0*(1-x0);x(k+1,:)=b.*x(k,:).*(1-x(k,:);(2) 运行以下程序，观察显示的图形，与题(1)的数据对照，注意收敛的倍周期。%图解法，图1和图2%文件名：p238fig_1_2.mclea

24、r; clc; close all;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;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');%实部为横坐标，虚部为纵坐标 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);%尽量不用循环 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) 运行程序并给出结果（对应不同的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，调整数据点图标的大小xlabel('b');ylabel('x (k)');grid on要求： 运行以上程序。 在运行结果的图形中，从对应的b值上的“点数”判断倍周期收敛。提示：放

28、大图形。 程序的运行结果（见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; %参数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，调整数据点图标的大小title('图4 模型的收敛、分岔和混沌');xlabel('参数 b');ylabel('x(k) （k足够大）');grid on; 运行以上m文件。运行结果（比较241图4）：5.4 2n倍周期图形p240242要求：在上题的程序中，修改b值，使运行后显示以下图形：(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值区间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倍周期时收敛b>bn_1，b<3.57y=zeros(1,2*2n);%行向量，用于存放xk最后的2×2n个值while d-c>1e-12 %使区间(c, d)足够小 b=(d+c)/2; x=0.2; %b取区间的中点 for i=1:100000 x=b*x*(1-x);

31、end y(1)=x; %取最后2×2n个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 %范数，比较 c=b;%满足2n倍周期收敛的b给区间(c, d)的下限c else d=b; %不满足2n倍周期收敛的b给区间(c, d)的上限d endendbmax=c; % 2n倍周期收敛的b的上限近似要求：编写程序，调用以上函数文件，求单倍周期、2倍周期、25倍周期收敛的b的上限近似值，输出要求有10位有效数字。运行结果输出形式如下：提示：可使用num2str函数。用下面的程序结构，会使运行速度加