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

.PAGE.实验08稳定性模型〔4学时〔第7章稳定性模型1.〔验证捕鱼业的持续收获——产量模型p215~219产量模型：其中,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×1sym>dF=diff<Fx,'x'>;%求导dF=simple<dF>%简化符号表达式%得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>=%若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<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为方程的平衡点%求方程的平衡点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代入符号表达式dFdFx0=simple<dFx0>%得F'<x0>=-r+E%求F'<x1>dFx1=subs<dF,x,x1>%得F'<x1>=r-E%若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<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>式%又由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.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%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'>;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列,即稳定条件。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.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>：☆<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>;%简化符号表达式pq=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.mclc;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;%此处填入一条正确语句endK=<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.mclear;clc;closeall;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;fork=1:100x<k+1,:>=f<x<k,:>,b>;endforn=1:length<b>figure<n>;%指定图形窗口figurensubplot<1,2,1>;%开始迭代的图形fplot<<x>[f<x,b<n>>,x],[0101]>;%x是自变量,画曲线y=bx<1-x>和直线y=xaxissquare;holdon;x0=x<1,n>;y0=0;%画迭代轨迹线fork=2:5x1=x<k,n>;y1=x<k,n>;plot<[x0+i*y0,x0+i*y1,x1+i*y1],'r'>;%实部为横坐标,虚部为纵坐标x0=x1;y0=y1;endtitle<['图解法：开始迭代的图形〔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%方程〔6xk+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.2fork=1:100x<k+1,:>=b.*x<k,:>.*<1-x<k,:>>;endplot<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:nx=b.*x.*<1-x>;ifk>=n-50xx=[xx;x];endendplot<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<b<3：★<2>2倍周期〔3<b<3.449：★<3>22倍周期〔3.449<b<3.544：★<4>23倍周期〔3.544<b<3.564：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.57y=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:100000x=b*x*<1-x>;endy<1>=x;%取最后2×2^n个xkfork=1:2^<n+1>-1y<k+1>=b*y<k>*<1-y<k>>;endifnorm<y<1:2^n>-y<2^n+1:2^<n+1>>><1e-12%范数,比较c=b;%满足2^n倍周期收敛的b给区间<c,d>的下限celsed=b;%不满足2^n倍周期收敛的b给区间<c,d>的上限dendendbmax=c;%2^n倍周期收敛的b的上限近似要求：编写程序,调用以上函数文件,求单倍周期、2倍周期、……、25倍周期收敛的b的上限近似值,输出要求有10位有效数字。运行结果输出形式如下：提示：可使用num2str函数。用下面的程序结构,会使运行速度加快。functionmain<>自己编写的程序；将上面的函数文件复制到此处。★运行的程序及输出结果：functionm