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

1实验08稳定性模型（4学时）（第7章稳定性模型）1.（验证）捕鱼业的持续收获——产量模型p215~219产量模型：()()1xxtFxrxExN其中，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)=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×1symdF=diff(Fx,'x');%求导dF=simple(dF)%简化符号表达式2%得F'(x)=%求F'(x0)并简化dFx0=subs(dF,x,x0);%将x=x0代入符号表达式dFdFx0=simple(dFx0)%得F'(x0)=%求F'(x1)dFx1=subs(dF,x,x1)%得F'(x1)=%若Er，有F'(x0)0，F'(x1)0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若Er，则结果正好相反。%在渔场鱼量稳定在x0的前提下（Er），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且hm=f(x0*)。symsrxNfx=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文件（见[215~217]）：%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的点x为方程的平衡点%求方程的平衡点3symsrxNE;%定义符号变量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×1symdF=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%若Er，有F'(x0)0，F'(x1)0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若Er，则结果正好相反。%在渔场鱼量稳定在x0的前提下（Er），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且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)式4%又由x0*=N(1-E/r)，可得E*=r/2，见P217(8)式%产量模型的结论是：%将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。2.（验证、编程）种群的相互竞争P222~228模型：)1()()1()(2211222222111111NxNxxrtxNxNxxrtx其中，x1(t),x2(t)分别是甲乙两个种群的数量。r1,r2是它们的固有增长率。N1,N2是它们的最大容量。σ1：单位数量乙（相对N2）消耗的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的σ1倍。对σ2可作相应解释。2.1（编程）稳定性分析p224~225要求：补充如下指出的程序段，然后运行该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中（见下面注释所示），P的结构类型为4×2sym，P的第1列对应x1，第2列对应x2。x1x2=[x1,x2]%显示结果disp('');P%调整位置后的4个平衡点：%P(1,:)=P1（N1，0）%P(2,:)=P2（0，N2）5%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);%简化符号表达式pq=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.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%编写出该程序段。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');6gx1=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);%简化符号表达式pq=subs(det(A),{x1,x2},{P(:,1),P(:,2)});q=simple(q);disp('');[Ppq]%显示结果%得到教材p225表1的前3列，经测算可得该表的第4列，即稳定条件。2.2（验证、编程）计算与验证p227微分方程组1211111212222212121212()(1)()(1)0.5,1.6,2.5,1.8,1.6,1xxxtrxNNxxxtrxNNrrNN(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]);7text(2.4,1.55,'x1(t)');text(2.2,0.25,'x2(t)');%函数文件名：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)/N2),r2*x(2)*(1-k2*x(1)/N1-x(2)/N2)]';☆(1)运行程序并给出结果（比较[227]图2(a)）：0246800.511.52x1(t)x2(t)☆(2)（验证）将x1(0)=1,x2(0)=2代入(1)中的程序并运行。给出结果（比较[227]图2(b)）：0246800.511.52x1(t)x2(t)(3)（编程）在同一图形窗口内，画(1)和(2)的相轨线，相轨线是以x1(t)为横坐标，x2(t)为纵坐标所得到的一条曲线。具体要求参照下图。800.511.5200.511.52(0.1,0.1)(1,2)图中的两条“点线”直线，一条的两个端点为(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;900.511.5200.511.52(0.1,0.