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(
2、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= , x1= ,见P216(4)式x0=x(2);x1=x(1);%符号变量x的结构类型成为%求F(x)的微分F(x)syms x; %定义符号变量x的结构类型为dF=diff(Fx,x); %求导dF=simple(dF) %简化符号表达式%得 F(x)= %求F(x0)并简化dFx0=subs(dF,x,x0); %将x=x0代入符号表达式dFdFx
3、0=simple(dFx0) %得 F (x0)= %求F (x1)dFx1=subs(dF,x,x1)%得 F (x1)= %若 Er,有F(x0)0,故x0点稳定,x1点不稳定(根据平衡点稳定性的准则);%若 Er,则结果正好相反。%在渔场鱼量稳定在x0的前提下(Er),求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*= ,见P217(6)式hm=subs(fx,x,x0)%得 hm
4、= ,见P217(7)式%又由 x0*=N(1-E/r),可得 E*= ,见P217(8)式%产量模型的结论是:%将捕捞率控制在固有增长率的一半(E=r/2)时,能够获得最大的持续产量。符号简化函数simple,变量替换函数sub的用法见提示。 给出填空后的M文件(见215217):%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)=0的
5、点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的结构类型成为%求F(x)的微分F(x)syms x; %定义符号变量x的结构类型为dF=diff(Fx,x); %求导dF=simple(dF) %简化符号表达式%得 F(x)= r-2*r*x/N-E %求F(x0)并简化dFx0=subs(dF,x,x0); %将x=
6、x0代入符号表达式dFdFx0=simple(dFx0) %得 F (x0)= -r+E %求F (x1)dFx1=subs(dF,x,x1)%得 F (x1)= r-E %若 Er,有F(x0)0,故x0点稳定,x1点不稳定(根据平衡点稳定性的准则);%若 Er,则结果正好相反。%在渔场鱼量稳定在x0的前提下(Er),求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 ,见
7、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)时,能够获得最大的持续产量。2.(验证、编程)种群的相互竞争P222228模型:其中,x1(t), x2(t)分别是甲乙两个种群的数量。r1, r2是它们的固有增长率。N1, N2是它们的最大容量。1:单位数量乙(相对N2)消耗的供养甲的食物量为单位数量甲(相对N1)消耗的供养甲的食物量的1倍。对2可作相应解释。2.1(编程)稳定性分析p224225要求:补
8、充如下指出的程序段,然后运行该m文件,对照教材上的相应结果。%7.3 种群的相互竞争稳定性分析%文件名: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编写出该程序段。提示(1)使用符号表达式;(2)用函数solve求解代数方程组,解放入x1, x2;(3)调整解(平衡点)的顺序放入P中(见下面注释所示),
9、P的结构类型为,P的第1列对应x1,第2列对应x2。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(g,x
10、1); 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 种群的相互竞争稳定性分析%文件名:p224_225.mclear; clc;%甲乙两个种群满足的增长方程:%
11、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);P=x1(2,3,4,1),x2(2,3,4,1);x1x2=x1,x2 %显示结果disp( ); P%
12、调整位置后的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(g,x1); gx2=diff(g,x2);disp( ); A=fx1,fx2;gx1,gx2 %显示结果p=
13、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(验证、编程)计算与验证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(p22
14、7,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.mfunction 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) (
15、验证)将x1(0)=1, x2(0)=2代入(1)中的程序并运行。给出结果(比较227图2(b)):(3)(编程)在同一图形窗口内,画(1)和(2)的相轨线,相轨线是以x1(t)为横坐标,x2(t)为纵坐标所得到的一条曲线。具体要求参照下图。图中的两条“点线”直线,一条的两个端点为(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)对应
16、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,( 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可作相应解释。
17、要求: 修改题2.1的程序,求模型的平衡点及稳定性。给出程序及其运行结果(比较229表1,注:只要最终结果):clear; clc;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);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(-(f
18、x1+gx2),x1,x2,P(:,1),P(:,2); %替换p=simple(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.0
19、2;xdot=(r-a*x(2).*x(1); (-d+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的收敛情况(实际上取k足
20、够大就可以了)。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,lengt
21、h(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.mclear; c
22、lc; 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; hold
23、 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(81:101,n)+i*x(81:101,n);%尽量不用循环 x
24、y(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), k=0,1,2,.clear; clc; close all;b=3.3 3.45 3.55 3.57;x=0.2*ones
25、(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值上的“点数”判断倍周期收敛。提示:放大图形。 程序的运行结果(见238表1):5.3 收敛、分岔和混沌p240242画出教材p241图4 模型的收敛、分岔和混沌。程序的m文件如下:%图4 模型(6)的收敛、分岔和混沌%文件名:p241fig4.m
26、%模型: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,markersize,5); %可改变值5,调整数据点图标的大小title(图4 模型的收敛、分岔和混沌);xlabel(参数 b);ylabel(x(k) (k足够大));grid on; 运行以上m文件。运行结果(比较241图4):5.4 2n倍周
27、期图形p240242要求:在上题的程序中,修改b值,使运行后显示以下图形:(1) 单周期(1b3):(2) 2倍周期(3b3.449):(3) 22倍周期(3.449b3.544):(4) 23倍周期(3.544bbn_1,b1e-12 %使区间(c, d)足够小 b=(d+c)/2; x=0.2; %b取区间的中点 for i=1:100000 x=b*x*(1-x); end y(1)=x; %取最后22n个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
28、;%满足2n倍周期收敛的b给区间(c, d)的下限c else d=b; %不满足2n倍周期收敛的b给区间(c, d)的上限d endendbmax=c; % 2n倍周期收敛的b的上限近似要求:编写程序,调用以上函数文件,求单倍周期、2倍周期、25倍周期收敛的b的上限近似值,输出要求有10位有效数字。运行结果输出形式如下:提示:可使用num2str函数。用下面的程序结构,会使运行速度加快。function main()自己编写的程序;将上面的函数文件复制到此处。 运行的程序及输出结果:function main()clc;n=5;%求单(赋0)、2倍(赋1)、2n倍周期收敛的b的上、下限B=ones(n+2, 1);
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1