数学建模实验答案稳定性模型.docx
- 文档编号:11383658
- 上传时间:2023-02-28
- 格式:DOCX
- 页数:56
- 大小:2.87MB
数学建模实验答案稳定性模型.docx
《数学建模实验答案稳定性模型.docx》由会员分享,可在线阅读,更多相关《数学建模实验答案稳定性模型.docx(56页珍藏版)》请在冰豆网上搜索。
数学建模实验答案稳定性模型
实验08稳定性模型(4学时)
(第7章稳定性模型)
1.(验证)捕鱼业的持续收获——产量模型p215~219
产量模型:
其中,
x(t)为t时刻渔场中的鱼量。
r是固有增长率。
N是环境容许的最大鱼量。
E是捕捞强度,即单位时间捕捞率。
要求:
运行下面的m文件,并把相应结果填空,即填入“_________”。
%7.1捕鱼业的持续收获——产量模型
%文件名:
p215_217.m
clear;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的点x为方程的平衡点
%求方程的平衡点
symsrxNE;%定义符号变量
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)
symsx;%定义符号变量x的结构类型为<1×1sym>
dF=diff(Fx,'x');%求导
dF=simple(dF)%简化符号表达式
%得F'(x)=
%求F'(x0)并简化
dFx0=subs(dF,x,x0);%将x=x0代入符号表达式dF
dFx0=simple(dFx0)
%得F'(x0)=
%求F'(x1)
dFx1=subs(dF,x,x1)
%得F'(x1)=
%若E
%若E>r,则结果正好相反。
%在渔场鱼量稳定在x0的前提下(E %通过分析(见教材p216图1),只需求x0*使f(x)达到最大,且hm=f(x0*)。 symsrxN fx=r*x*(1-x/N);%fx=dF df=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文件(见[215~217]): %7.1捕鱼业的持续收获——产量模型 %文件名: p215_217.m clear;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的点x为方程的平衡点 %求方程的平衡点 symsrxNE;%定义符号变量 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×1sym> %求F(x)的微分F'(x) symsx;%定义符号变量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代入符号表达式dF dFx0=simple(dFx0) %得F'(x0)=-r+E %求F'(x1) dFx1=subs(dF,x,x1) %得F'(x1)=r-E %若E %若E>r,则结果正好相反。 %在渔场鱼量稳定在x0的前提下(E %通过分析(见教材p216图1),只需求x0*使f(x)达到最大,且hm=f(x0*)。 symsrxN fx=r*x*(1-x/N);%fx=dF df=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)时,能够获得最大的持续产量。 2.(验证、编程)种群的相互竞争P222~228 模型: 其中, x1(t),x2(t)分别是甲乙两个种群的数量。 r1,r2是它们的固有增长率。 N1,N2是它们的最大容量。 σ1: 单位数量乙(相对N2)消耗的供养甲的食物量为单位数量甲(相对N1)消耗的供养甲的食物量的σ1倍。 对σ2可作相应解释。 2.1(编程)稳定性分析p224~225 要求: 补充如下指出的程序段,然后运行该m文件,对照教材上的相应结果。 %7.3种群的相互竞争——稳定性分析 %文件名: p224_225.m clear;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。 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)式。 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+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('');[Ppq]%显示结果 %得到教材p225表1的前3列,经测算可得该表的第4列,即稳定条件。 ★补充后的程序和运行结果(见[225]表1): 2341 %7.3种群的相互竞争——稳定性分析 %文件名: p224_225.m clear;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 %编写出该程序段。 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(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)式。 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+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('');[Ppq]%显示结果 %得到教材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('p227',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.m functiony=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)为纵坐标所得到的一条曲线。 具体要求参照下图。 图中的两条“点线”直线,一条的两个端点为(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)'); holdon; 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');%画两条直线 gridon; 3.(编程)种群的相互依存——稳定性分析P228~229 模型: 其中, x1(t),x2(t)分别是甲乙两个种群的数量。 r1,r2是它们的固有增长率。 N1,N2是它们的最大容量。 σ1: 单位数量乙(相对N2)提供的供养甲的食物量为单位数量甲(相对N1)消耗的供养甲的食物量的σ1倍。 对σ2可作相应解释。 要求: ☆修改题2.1的程序,求模型的平衡点及稳定性。 给出程序及其运行结果(比较[229]表1,注: 只要最终结果): clear;clc; symsx1x2r1r2N1N2k1k2; 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])]; symsx1x2;%重新定义 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);%简化符号表达式p q=subs(det(A),{x1,x2},{P(: 1),P(: 2)}); q=simple(q); [Ppq]%显示结果 4.(验证)食饵-捕食者模型p230~232 模型的方程: 要求: 设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文件: 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]=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.(验证)差分形式的阻滞增长模型p236~242 阻滞增长模型用微分方程描述为: 也可用差分方程描述为: 上式可简化为一阶非线性差分方程: 考察给定b和x0值后,当k→∞时,xk的收敛情况(实际上取k足够大就可以了)。 5.1数值解法和图解法p238~240 (1)取x0=0.2,分别取b=1.7,2.6,3.3,3.45,3.55,3.57,对方程 计算出x1~x100的值,显示x81~x100的值。 观察收敛与否。 (结果与教材p238~239表1比较) 下面是计算程序,在两处下划线的位置填入满足要求的内容。 %不同b值下方程x(k+1)=b*x(k)*(1-x(k)),k=0,1,2,...的计算结果 %文件名: p238table_1.m clc;clearall;formatcompact; b=[1.7,2.6,3.3,3.45,3.55,3.57]; x=zeros(100,length(b)); x0=0.2;%初值 ;%此处填入一条正确语句 fork=1: 99 ;%此处填入一条正确语句 end K=(81: 100)’;%将取81~100行 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.m clear;clc;closeall; f=@(x,b)b.*x.*(1-x);%定义f是函数的句柄,函数b*x*(1-x)含两个变量x,b b=[1.7,2.6,3.3,3.45,3.55,3.57]; x=ones(101,length(b)); x(1,: )=0.2; fork=1: 100 x(k+1,: )=f(x(k,: ),b); end forn=1: length(b) figure(n);%指定图形窗口figuren subplot(1,2,1);%开始迭代的图形 fplot(@(x)[f(x,b(n)),x],[0101]);%x是自变量,画曲线y=bx(1-x)和直线y=x axissquare;holdon; x0=x(1,n);y0=0;%画迭代轨迹线 fork=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))')']); 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);%尽量不用循环 xy(2: 2: 40)=x(81: 100,n)+i*x(82: 101,n); plot(xy,'r'); title(['图解法: 最后迭代的图形(b='num2str(b(n))')']); holdoff; end ☆ (2)运行程序并给出结果(对应不同的b值)见[238]图1、2: 20倍 20倍 21倍 22倍 23倍 混沌 5.2b值下的收敛图形p238~240 下面程序是在不同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.2 fork=1: 100 x(k+1,: )=b.*x(k,: ).*(1-x(k,: )); end plot(b,x(82: 101,: ),'.r','markersize',5);%可改变值5,调整数据点图标的大小 xlabel('b'); ylabel('x(k)'); gridon 要求: ①运行以上程序。 ②在运行结果的图形中,从对应的b值上的“点数”判断倍周期收敛。 提示: 放大图形。 ☆程序的运行结果(见[238]表1): 5.3收敛、分岔和混沌p240~242 画出教材p241图4模型的收敛、分岔和混沌。 程序的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: n x=b.*x.*(1-x); ifk>=n-50 xx=[xx;x]; end end plot(b,xx,'.b','markersize',5);%可改变值5,调整数据点图标的大小 title('图4模型的收敛、分岔和混沌'); xlabel('参数b'); ylabel('x(k)(k足够大)'); gridon; ☆运行以上m文件。 运行结果(比较[241]图4): 5.42n倍周期图形p240~242 要求: 在上题的程序中,修改b值,使运行后显示以下图形: ★ (1)单周期(1 ★ (2)2倍周期(3 ★(3)22倍周期(3.449 ★(4)23倍周期(3.544 5.5(编程)求2n倍周期的b值区间p240~241 求2^n倍周期收敛的b的上限的程序如下: functionbmax=fun(bn_1,n) %求2^n倍周期收敛的b的上限。 %bn_1是2^(n-1)倍周期收敛的b的上限(或2^n倍周期收敛的b的下限)。 c=bn_1;d=3.57;%2^n倍周期时收敛b>bn_1,b<3.57 y=zeros(1,2*2^n);%行向量,用于存放xk最后的2×2^n个值 whiled-c>1e-12%使区间(c,d)足够小 b=(d+c)/2;x=0.2;%b取区间的中点 fori=1: 100000 x=b*x*(1-x); end y (1)=x;%取最后2×2^n个xk fork=1: 2^(n+1)-1 y(k+1)=b*y(k)*(1-y(k)); end ifnorm(y(1: 2^n)-y(2^n+1: 2^(n+1)))<1e-12%范数,比较 c=b;%满足2^n倍周期收敛的b给区间(c,d)的下限c else d=b;%不满足2^n倍周期收敛的b给区间(c,d)的上限d end end bmax=c;%2^n倍周期收敛的b的上限近似 要求: 编写程序,调用以上函数文件,求单倍周期、2倍周期、……、25倍周期收敛的b的上限近似值,输出要求有10位有效数字。 运行结果输出形式如下: 提示: 可使用num2str函数。 用下面的程序结构,会使运行速度加快。 functionmain() 自己编写的程序; 将上面的函数文件复制到此处。 ★运行的程序及输出结果: functionmain() clc; n=5;%求单(赋0)、2倍(赋1)、…、2^n倍周期收敛的b的上、下限 B=ones(n+2,1);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 建模 实验 答案 稳定性 模型