数学建模试验答案稳定性模型

1、实验08稳定性模型(4学时)(第7章稳定性模型)1.(验证)捕鱼业的持续收获一一产量模型p215219产量模型：X二F(x)=rx1一年-Ex其中，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为方程的平衡点%求方程的平衡点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的结构类型成为<2Msym>%求F(x)的微分F'(x)symsx;%定义符号变量x的结构类型为<1Msym>dF=diff(Fx,'x');%求导dF=simple(dF)%简化符号表达式%得F'(x)=%求F'(x0)并简化dFx

3、0=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(x0M到最大hm%通过分析(见教材p216图1),只需求x0*使f(x)达到最大，且hm=f(x0*)。symsrxNfx=r*x*(1-x/N);%fx

4、=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;%无捕捞条件下单位时间的增长量：f(x)=rx(1-x/N)%捕捞条件下单位时间的捕捞

5、量：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的结构类型成为<2Msym>%求F(x)的微分F'(x)symsx;%定义符号变量x的结构类型为<1Msym>

6、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的前提

7、下(E<r),求E使持续产量h(x0M到最大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.(验证、编程)种群的相互竞争P22222

8、8模型：xi(t)=riXi(1。1),NiN2XiX2、X2(t)-r2X2(i-2)NiN2其中，Xl(t),X2(t)分别是甲乙两个种群的数量。ri,r2是它们的固有增长率。Ni,N2是它们的最大容量。E：单位数量乙(相对N2)消耗的供养甲的食物量为单位数量甲(相对N1)消耗的供养甲的食物量的E倍。对d2可作相应解释。2.i(编程)稳定性分析p224225要求：补充如下指出的程序段，然后运行该m文件，对照教材上的相应结果%7.3种群的相互竞争一一稳定性分析%文件名：p224_225.mclear;clc;%甲乙两个种群满足的增长方程：%x1,(t)=f(x1,x2)=r1*x1*(1-x

9、1/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的结构类型为<4X2sym>,P的第1列对应x1,第2列对应x2。x1x2=x1,x2%显示结果disp('');P%调整位置后的4个平衡点：%P(1,:)=P1(N1,0)%P(2,:)=P2(0,N2)

10、%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(-(fx

11、1+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,x

12、2)=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(N

13、1*(-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,

14、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微分方程组xi(t)=r1x1(1-X2NiN2X2(t)=许(1-二X1X22N1N201二0.5,二1.6,1=2.5,R=1.8,M=1.6,N2=1(1)(验证)当x1(0)=x2(0)=0.1时，求微分方程的数值解，将解的数值分别画出x1(t)和x2(t)的曲线，它们同

15、在一个图形窗口中。程序：%命令文件ts=0:0.2:8;x0=0.1,0.1;t,x=ode45cp227',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

16、)/N2),r2*x(2)*(1-k2*x(1)/N1-x(2)/N2)'(1)运行程序并给出结果(比较227图2(a):(2)(验证)将xi(0)=1,X2(0)=2代入(1)中的程序并运行。给出结果(比较227图2(b):(3)(编程)在同一图形窗口内，画(1)和(2)的相轨线，相轨线是以xi(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=od

17、e45cp227',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=ode45cp227',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.(编程)种群的相互依存一一稳定性分析P228229模型：xi(t)=1为(1

18、-7+。1告)jNiN2iX1X2X2(t)=2X2(7)JNiN2其中，X1（t）,X2（t）分别是甲乙两个种群的数量。ri,r2是它们的固有增长率。Ni,N2是它们的最大容量。01：单位数量乙（相对N2）提供的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的E倍。对d2可作相应解释。要求：修改题2.1的程序，求模型的平衡点及稳定性。给出程序及其运行结果（比较2291表1,注：只要最终结果）：clear;clc;symsx1x2r1r2N1N2k1k2;f=r1*x1*(1-x1/N1+k1*x2/N2);g=r2*x2*(-1+k2*x1/N1-x2/N2);x1,x2=sol

19、ve(f,g);P=x1(2,4,1,3),x2(2,4,1,3);symsx1x2;%Uffcr义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%显示结果j匚4|bbj1ndif!PX|匕I.fi

20、ai/吨"i#duh*?iJid"rIfvll?十HL屯r)irS-IdTrN11心叱-(kL-Dj旧产2-7r4rl腔-ME*k>r：)/fJi+UE1J*W+r拄皿-L)>/(k+k?-L)10,0,r!-ri*-rl»r2)Q-IE,kl'rl-r：-rL-rl>fT2r(kl-1)1fit»|4.(验证)食饵-捕食者模型p230232模型的方程：X(t)=rx-axyx(0)=%y(t)=-dybxyy(0)=Vo要求：设r=1,d=0.5,a=0.1,b=0.02,xo=25,yo=2。输入p231的程序并运行，结果与

21、教材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=ode45fshier',ts,x0);t,x,plot(t,x),grid,gtext('x(t)'),gtext('y(t)'),%运行中在图上标注pause,plot(x(:,1),x(:,2),grid,x(t

22、),y(t)图形：相轨线y(x)图形:5.(验证)差分形式的阻滞增长模型p236242阻滞增长模型用微分方程描述为：xx(t)=rx(1-)N也可用差分方程描述为：yk1-yk=ryk(17),k=0,1,2,111N上式可简化为一阶非线性差分方程：Xk.1=bXk(1-Xk),k=0,1,2JH考察给定b和x0值后，当k一训寸,xk的收敛情况(实际上取k足够大就可以了)。5.1数值解法和图解法p238240(1)取xo=0.2,分别取b=1.7,2.6,3.3,3.45,3.55,3.57对方程Xki=bXk(1-Xk),k=0,1,2,H|计算出X1X100的值，显示X81X100的值。观

23、察收敛与否。(结果与教材p238239表1比较)下面是计算程序，在两处下划线的位置填入满足要求的内容。%不同b值下方程x(k+1)=b*x(k)*(1-x(k),k=0,1,2,.的计算结果%文件名：p238table_1.mclc;clearall;formatcompactb=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)'%将取81100行disp(num2str(NaN,b;K,x(K,:),4);%取4位有效数字，

24、NaN表示不是数值(1)对程序正确填空，然后运行。填入的正确语句和程序的运行结果(对应不同的b值)见238表1:x(1,:)=b*x0*(1-x0);x(k+1,:)=b.*x(k,:).*(1-x(k,:);CouandVin.dovE3叵区1.91.72.63.33.4E3.553.57310.4113D.6L5447940,4475-0.MS0.475gS20.U18CU0.82360.8530.88740.B903出0.411M口54Qt4794U.43Z&.35二U.MSBK40.411K0.8236D.R4SS0,E12T0.nQ6950.4118D.1540.4T网Ol4

25、4750.G40St).5死鸵Cl41ISD.61540.82360,8530.83170.羽，S370.11130,47940.43270.37030.转5口S30.4113U6LE40.82360.82730.E27B39.41IBD.6154Ci.47940.44打CJ.506U.5DF9WH.4118OiLE40.8236凡函mQ.S3740.B9Q2910.41IS61540474Cl4327(J.犯4HCLM二口常0.41130161sd0,825o.三处g0.31270.亚930.41IEDl6L540.47940l447j0.5M50.5607%0.41130,S23CG.阴M

26、0.8S1T。冏柒第0.41ISCk61B40,47940.43270.37030.3TBB960,*11180,61540.8236以总归0.8278.,84四0.41ISD,6154fi.4740,4471r,EC60.479?第0.4113口540.82560.S&30.距X0.89L的0.4113-G1540.47M0l432?0.354B0.34661000.1113(X61B4Q,鸵蒸0.812T0.3055(2)运行以下程序，观察显示的图形，与题(1)的数据对照，注意收敛的倍周期。%图解法，图1和图2%文件名：p238fig_1_2.mclear;clc;closeall

27、;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;%画迭代轨迹线for

28、k=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)

29、;plot(xy,'r');title(图解法：最后迭代的图形(b='num2str(b(n)')');holdoff;end(2)运行程序并给出结果(对应不同的b值)见238图1、2:图解法：最后迭代的图形(b=1.7)10.80.60.40.200.5120倍0图解法：开始迭代的图形（b=2.6）10.80.60.40.20图解法：最后迭代的图形（b=2.6）10.80.60.40.2000.5100.5120倍00.5121倍图解法：最后迭代的图形（b=3.3）10.80.60.40.20图解法：开始迭代的图形（b=3.45）图解法：最后迭代的图形

30、（b=3.45）图解法：开始迭代的图形（b=3.55）23倍图解法：最后迭代的图形(b=3.55)1.0.80.60.40.2000.510.80.60.40.200.5图解法：开始迭代的图形（b=3.57）1图解法：最后迭代的图形（b=3.57）10.80.60.40.2000.51混沌5.2b值下的收敛图形p238240下面程序是在不同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);%初

31、值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收敛、分岔和混沌p240242画出教材p241图4模型的收敛、分岔和混沌程序的m文件如下:%图4模型(6)的收敛、分岔和混

32、沌%文件名：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文件。

33、运行结果(比较241图4):fjehebiL'r|x1=1.XdaiVaw.ZxixhtIT0:IvD-aiJirtcrp工ndnr世1窜?白N、工"？*公口口国I国5.42n倍周期图形p240242要求：在上题的程序中，修改b值，使运行后显示以下图形:(1)单周期(1<b<3):>FLgUZEt匚巴摩EileEditVie*lasertloolsIesJrtottindo*HelpLfdA玲K1忑为立1X二总0"圉日模型的收被、分岔和混他U>；111r；：!:；,.,1',_：-：_r：；：:，Q.D-11Jitlllll；!II!

34、I1IBIIIIIIB|Je|IIIIIIJI(2)2倍周期(3<b<3.449):/FLgUZEt口叵仅.91.1e*Insert工。工与Rudrt<xptindo*Help，k鼠一2立130因口模型的收敛“分岔和您随3.053.13.153.23.253.33.353.43.4536参抽b576E6了0.0&Doo(1<卷嗖)等(3)22倍周期(3.449<b<3.544):>FlgUEE1fZT|XEileIditInsertIoo>ls取drt叮tindwrHelp日白U)珞-：A-QF®配3ns"0圉日模型的

35、收敛、分岔和混沌(4)23倍周期(3.544<b<3.564):F£tuie1L叵|区EditVie*Insert工口o1工上Rudrta;ptindo*HelpLRJU，玲.鼠-T盟总S图b模型的收散、分茗和混沌5.5（编程）求2n倍周期的b值区间p240241求2An倍周期U敛的b的上限的程序如下：functionbmax=fun(bn_1,n)%求2入门倍周期收敛的b的上限。%bn_1是2A(n-1)倍周期收敛的b的上限(或2=倍周期收敛的b的下限)c=bn_1;d=3.57;%2倍周期时收敛b>bn_1,b<3.57y=zeros(1,2*2An);%

36、行向量，用于存放xk最后的2X2n个值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八门个xkfork=1:2A(n+1)-1y(k+1)=b*y(k)*(1-y(k);endifnorm(y(1:2An)-y(2An+1:2A(n+1)<1e-12%范数，比较c=b;%满足2An倍周期收敛的b给区间(c,d)的下限celsed=b;%不满足2八门倍周期收敛的b给区间(c,d)的上限dendendbmax=c;%2八门倍周期收敛的b的上限近

37、似要求:编写程序，调用以上函数文件，求单倍周期、2倍周期、25倍周期收敛的b的上限近似值，输出要求有10位有效数字。运行结果输出形式如下:提示：可使用num2str函数。用下面的程序结构，会使运行速度加快。functionmain()自己编写的程序；将上面的函数文件复制到此处。运行的程序及输出结果：functionmain()clc;n=5;%求单(赋0)、2倍(赋1)、2倍周期收敛的b的上、下限B=ones(n+2,1);fori=0:nB(i+2)=fun(B(i+1),i);end%单周期收敛时，B(1)<b<B(2);2倍时，B(2)<b<B(3);disp(n

38、um2str(-1:n)',B,10);%10位有效数字functionbmax=fun(bn_1,n)%求2入门倍周期收敛的b的上限。%bn_1是2A(n-1)倍周期收敛的b的上限(或2=倍周期收敛的b的下限)c=bn_1;d=3.57;%2倍周期时收敛b>bn_1,b<3.57y=zeros(1,2*2An);%行向量，用于存放xk最后的2X2n个值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八门个xkfork=1:2A

39、(n+1)-1y(k+1)=b*y(k)*(1-y(k);endifnorm(y(1:2An)-y(2An+1:2A(n+1)<1e-12%范数，比较c=b;%满足2An倍周期收敛的b给区间(c,d)的下限celsed=b;%不满足2八门倍周期收敛的b给区间(c,d)的上限dendendbmax=c;%2八门倍周期收敛的b的上限近似conaiid一.rz-iFn'irx£d七口也。耳3-1102.99976909713.44935596723.54402147133.56435787943.56873999353.569679465A»|注：vpa的用法附1:

40、实验提示第1题提示：符号简化函数simple,变量替换函数sub符号简化函数simple的格式：simple(S)对符号表达式S尝试多种不同的算法简化，以显示S表达式的长度最短的简化形式变量替换函数sub的格式:subs(S,OLD,NEW)将符号表达式S中的OLD变量替换为NEW变量附2:第7章稳定性模型2157.1捕鱼业的持续收获7第/邕稳定性模型虽然动态过程的变化规律一般要用微分方程建立的动态模型来描述，但是对于某些实际问题，建模的主要目的并不是要寻求动态过程每个瞬时的性态，而是研究某种意义下稳定状态的特征，特别是当时间充分长以后动态过程的变化趋势.比如，在什么条件下描述过程的变量会越来

41、越接近某些确定的数值，在什么情况下又会越来越远离这些数值而导致过程不稳定.为了分析这种稳定与不稳定的规律，常常不需要求解微分方程（并且我们将会看到，即使对于不太复杂的方程，解析解也不是总能得到的），而可以利用微分方程稳定性理论，直接研究平衡状态的稳定性就行了.在下面建模过程中，大多数情况下我们将直接引用微分方程稳定性理论的结果.不熟悉这方面内容的读者可以参阅17节.差分方程的稳定性与微分方程有很多相似之处.也有明显的差别，工5和工7节将介绍这方面的内容.7.1捕鱼业的持续收获可持续发展是一项基本国策，对于像渔业、林业这样的再生资源，一定要注意适度开发，不能为了一时的高产去“竭泽而渔”，应该在持

42、续稳产的前提下追求产量或效益的最优化.考察一个渔场，其中的鱼量在天然环境下按一定规律增及，如果捕捞量恰好等于增长最，那么渔场鱼量将保持不变，这个捕捞量就可以持续下去本节要建立在捕播情况下渔场鱼鼠遵从的方程，分析鱼址稳定的条件，并且在稳定的前提下讨论如何控制捕捞使持续产量或经济效益达到最大，最后研究所谓捕捞过度的问题产模型记时刻t渔场中鱼量为X。，关于的自然增长和人工捕捞作如下假设：L在无捕捞条件下的增氏服从1。机比规律（见5.6节的阻滞增长模型），即*(0=f(x)=rx|-N/(1)是固有增长率，川是环境容许的最大鱼最，用人工）表示单位时间的增长量.2.单位时间的捕捞量（即产量）与渔场鱼量工

43、成正比，比例常数E表示单位时间捕捞率,又称为捕捞强度，可以用比如捕鱼网眼的大小或出海渔船数量来控制其大小.于是单位时间的捕捞量为(2)h（x）-Ex根据以上假设并记尸（#）=/（）-A（x）得到捕措情况下渔场鱼量满足的方程x（t）=F（富）=1_木）-Ex我们并不需要解方程（3）以得到工（。的动态变化过程.同希望知道渔场的稳定鱼量和保持稳定的条件，即时间上足够长以后渔场鱼最工。）的趋向，并由此确定最大持续产量.为此可以直接求方程（3）的平衡点并分析其稳定性.令Ex-0得到两个平衡点2=N(1-g,J?!=0不难算出F'(不)=E-r.k(盯)-r-E所以若£<r（5）有

44、尸（如）<0,尸（盯）>0,故通点稳定百点不稳定（判断平衡点稳定性的准则见7.7节）彳若E>r,则结果正好相反.y-f(x)图】最大持续产拼的图解法xf=M2xDNE是捕捞率，r是最大的增长率，上述分析表明，只要捕捞适度（£<1）.就可使渔场鱼量稳定在质，从而获得持续产最M不）二取°而当捕捞过度时（E>力，渔场鱼量将趋向盯=0,当然谈不上获得持续产量了.进一步讨论渔场鱼量稳定在与的前提下，如何控制捕捞强度E使持续产量最大的问瓯用图解法可以非常简单地得到216；结果.根据（1）,（2）式作抛物线耳）和直线=乂*）=&,如图1.注意到y=/

45、«）在原点的切线为？二%所以在条件（5）下,二株必与有交点P,P的横坐标就是稳定平衡点小.根据假设2,尸点的纵坐标为为稳定条件下单位时间的持续产量.由图1立刻知道，当y=&与在抛物线顶点P*相交时可获得最大的持续产量,此时的稳定平衡点为(6)且单位时间的最大持续产量为而由（4）式不难算出保持渔场鱼量稳定在工；的捕捞率为E=2综上所述，产量模型的结论是将捕捞率控制在固有增长率r的一半，更简单一些，可以说使渔场鱼量保持在最大鱼量N的一半时，能够获得最大的持续产量.效益模型从经济角度看不应追求产量最大，而应考虑效益最佳.如果经济效益用从捕捞所得的收入中扣除开支后的利润来衡量，并且简

46、单地假设：鱼的销售单价为常数P,单位捕捞率（如每条出海渔船）的费用为常数c,那么单位时间的收入T和支出S分别为T-ph(x)=pEx,S=cE(9)单位时间的利润为RT-S-pEx-cE在稳定条件*=褊下，以（4）代入（10）式，得(10)/?(£)=T(E)-S(£)=pNE(T-(11)用微分法容易求出使利润K（E）达到最大的捕捞强度为(12)将4代入（4）式，可得最大利润下的渔场稳定鱼量5及单位时间的持续产量K为(13)(14)217将（12）14）式与产量模型中的（6）（8）式相比较可以看出,在最大效益原则下捕捞率和持续产量均有所减少，而渣场应保持的稳定鱼量有所增加

47、，并且减少或增加的部分随着捕捞成本e的增长而变大，循着销售价格p的增长而变小.请读者分析这些结果是符合实际情况的.捕措过度上面的效益模型是以计划捕捞（或称封闭式捕捞）为基础的，即渔场由单独的经营者有计划地捕捞，可以追求最大利润.如果渔场向众多盲目的经营者开放，比如在公海上无规则地捕捞,那么即使只有微薄的利润，经营者也会蟒拥而去，这种情况称为盲目捕捞（或开放式捕捞）.这种捕捞方式将导致捕捞过度，下面讨论这个模型.（11）式给出了利润与捕捞强度的关系令R（E）=0的解为跖，可得(15)当E风时，利润靛（E）0,盲目的经营者们会加大捕措强度；若EaE-利润R（E）<0*他们当然要减小强度,所以

48、乙是盲目捕捞下的临界强度M也可由图解法确定.在图2中以£为横坐标，按（11）式画出T（E）和S（E）,它们交点的横坐标即为名（图2中的号1或后加）,由（15）式或图2容易知道，酌存在的必要条件（即4>。）是图2盲目捕捞强度的图解法即售价大于（相对于总量而言）成本.并且由（15）式可知，成本越低，售价越高，则瓦越大,将（15）代人（4）式，得到盲目捕捞下的渔场稳定鱼量为（17）%完全由成本-价格比决定，随着价格的上升和成本的下降%将迅速减少，出现捕捞过度.比校32）和（15）式可知=2号,即盲目捕捞强度比最大效益下捕捞噩度大一倍.从（15）式和图2还可以得到，当竟<p<

49、;2时，（昂<）/</，如图2中218凰，称经济学捕捞过度；当P>2右时，曷>£，如图2中”，称生态学捕捞过度.:*本节完*评注为了研究渔业的产量、效益及捕捞过度问题，首先在对鱼的自然增长和捕捞情况的合理假设F,建立渔场鱼量的基本方程（3）,并利用平衡点稳定性分析确定了保持渔场鱼鼠楫定的条件.产量、效益和捕捞过度3个模型在稳定的前提下步步深入，数学推导过程十分简单，却得到了在定性关系上与实际情况完全符合的结果.如果改变对鱼的自然增长和人工捕捞的假设，模型及结果将随之变化（习题1,2）.72军备竞赛两个国家或两个国家集团之间由于相互不信任和各种矛盾的存在、发展而

50、不断增加自己的军事力量，防御对方可能发动的战争.能否用一个数学模型描述这种军备竟赛的过程，从定性和定髭的角度对竞赛的结果作出解释或预测.本节介绍LF.Richardsgl939年提出的一个模型.当然影响军备竞赛的因素是错综复杂的，无法用数学工具给以恰当的圆满的描述.这个模型只不过告诉我们，一个复杂的实际过程可以被合理地简化到什么程度，得到的结果又怎样用来解释实际现象，并且如果允许军备比赛沿着机械的本能和固定的模式发展的话.会导致什么样的结局廊.模型假设与构成为了方便起见，用军备这个词表示军事力量的总和，如兵力、装备、军事颈算等.甲乙双方在时刻t的军备分别记作工3）和江。，假定它们的变化只取决于

51、下面3个因素工1,由于相互不信任及矛盾的发展，一方军备越大，另一方军备增加得越快.2,由于各方本身经济实力的限制.任一方军备越大，对军备增长的制约作用越大.3由于相互敌视或领土争端，每一方都存在膂增加军备的固有潜力.进一步假定前两个因素的影响是线性的，第3个因素的影响是常数，那么双上）和贝。的变化过程可用微分方程组：£（£）=-5+Gy+g（I）y（0三自=By+%表示，其中的系数均大于或等于0.是对方军备剌激程度的度最e,©是己方经济实力制约程度的度量洒力是巴方军备竞赛的固有潜力-如果我们感兴趣的是军备竞赛的结局由什么因素决定，而不关心竞赛的过程.那么只需用微分

52、方程稳定性理论讨论时间充分长以后工（。，¥”）的变化趋势，即方程（1）的平衡点的稳定情况令（1）式右端等于0,容易算出平衡点（与，九）为fl印-H'九一4一奴方程（1）的系数矩阵为F-a*1A=>，p于是按照判断平衡点稳定性的方法计算（见7.7节（9）-（13）式）p=-（an+a12）=a>0（3）qdetA=ap-kl（4）由稳定性准则（见7.7节（15）式八当印奴（5）时，平衡点（与，为）是稳定的：反之.是不稳定的.这就是说，在（5）式的条件下，时间足够长以后双方的军备将分别趋向一个有限值，军备竞赛是稳定的、模型的定性解释根据方程（1）和平衡点稳定性的分析，

53、可以解释几个简单而又重要的现象.1.条件（5）表明，当双方的经济制约程度邛大于双方的军备刺激程度kl时，军备竞赛才会趋向稳定.反之声U）4（f）将趋向无穷，竞赛无限地进行下去，可能导致战争.2,由（2）式，如果g=A=0,则0=0,兀=0是方程（1）的平衡点，并且在条件（5）下它是稳定的，于是如果在某个时候品有=y（/0）=。，%了就永远保持为E这种情况可以解释为双方不存在任何敌视和争端，通过裁军可以达到持久和平.两个友好的邻国正是这样.3 .如果即使由于某种原因（如裁军协定）在某个时候双方军备大减，不妨设M%）=y（r0）=。，那么因为x=h,也将使双方重整军备.这说明未经和解的裁军（即不消

54、除敌视或领土争端）是不会持久的.4 .如果由于某种原因（如战败或协议），在某个时候一方的军备大减，不妨设立色）=0,那么因为五二行+心也将使该方重整军备.这说明存在不信任（无*0）或固有争端3*0）的单方面裁军也不会持久.模型参数的估计为了利用（5）式判断军备竞赛是否会趋于稳定，需要知道白班,AN的数直估计这些参数无疑是很困难的，下面是RichardQi提出的一种方法.1. 的估计设工（0）=0,当二较小时，忽略g和-ox的作用,并近似地假定丁二为不变，由方程（1）得i二%220如果当t=丁时，则由（6）式得到（7）这说明A"是甲方军备从0到赶上乙方军备外所需的时间.例如，谯国从19

55、33年开始重整军备，只用了约3年的时间就赶上了它的邻国.假设它增加军备的固有潜力g被制约效应ax所抵消，那么可以认为德国的A-=3年，即上=0.3./可以类似地估计，或者合理地假定它与国家的经济实力成正比.这样若乙国的经济实力是德国的2倍，则可以估计20.6.嘉a,尸的估计设g=0,y=0,由方程（】）可得x（t）=.（0”7以代人，算出这表示仪7是在乙方无军备时甲方军备减少到原来的,所需的时间.Richardson认e为这大概是一个国家议会的任期，对于议会任期5年的国家来说©=0.2.对模型和叁数的粗晡检览考察第一次世界大战前夕，欧洲的两个国家同盟法俄同盟和德奥匈同盟的军备竞赛情况.-两个同盟的经济实力大致相等,且约为德国的3倍，因为德国的A=0.3,所以这两个同盟的左=/