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

实验08稳定性模型(4学时)

(第7章稳定性模型)

1.(验证)捕鱼业的持续收获一产量模型p215~219

产量模型：

x(t)=F(x)=rx1区-Ex,IN)

其中，

x(t)为t时刻渔场中的鱼量。

r是固有增长率。

N是环境容许的最大鱼量。

E是捕捞强度，即单位时间捕捞率。

要求：

运行下面的m文件，并把相应结果填空，即填入

%7.1捕鱼业的持续收获一■产量模型--------------------------------

%文件名：p215_217.m

clear;cic;

%无捕捞条件下单位时间的增长量：f(x)=rx(l-x/N)

%捕捞条件下单位时间的捕捞量：h(x)=Ex

%F(x)=f(x)-h(x)=rx(l-x/N)-Ex

%捕捞情况下渔场鱼量满足的方程：x'(t)=F(x)

%满足F(x)=0的点x为方程的平衡点

%求方程的平衡点

symsrxNE;%定义符号变量

Fx=r*x*(l-x/N)-E*x;%创建符号表达式

x=solve(Fx,x)%求解F(x)=O(求根)

%得到两个平衡点，记为:

%xO=,xl=_,见直16(4)式

x0=x(2);

米1=(1);%符号变量x的结构类型成为V2xlsym>

%求¥值)的微分F，(x)

symsx;%定义符号变量x的结构类型为<1*1§眄>

dF=diff(Fx,'x');%求导

dF=simple(dF)%简化符号表达式

%得F'(x)=

%求F(xO)并简化

dFxO=subs(dFaxO);%将x=xO代入符号表达式dF

dFxO=simple(dFxO)

%得F'(xO)=

%求F'(xl)

dFxl=subs(dF,x,xl)

%得F'(xl)=

%若£<1•，有F'(x0)<0,F'(xl)>0,故xO点稳定，xl点不稳定(根据平衡点稳

定性的准则)；

%若E>r,则结果正好相反。

%在渔场鱼量稳定在xO的前提下(E<r)，求E使持续产量h(xO)达到最大hm。

%通过分析(见教材p216图1),只需求xO*使f(x)达到最大，且hm=f(xO*)0

symsrxN

fx=r*x*(l-x/N)；%fx=dF

df=diff(fx,'x);

xO=soive(df,x)

%得x0*=,见P217(6)式

hm=subs(fx,x,xO)

%得11111=，见P217(7)式

%又由xO*=N(l-E/r),可得E*=,见P217(8A

%产量模型的结论是：

%将捕捞率控制在固有增长率的一半(E=r/2)时，能够获得最大的持续产量。

符号简化函数simple,变量替换函数sub的用法见提示。

★给出填空后的M文件(见[215〜217]):

%7.1捕鱼业的持续收获一一产量模型

%文件名:p215_217.m

clear;clc;

%无捕捞条件下单位时间的增长量：f(x)=rx(l-x/N)

%捕捞条件下单位时间的捕捞量：h(x)=Ex

%F(x)=f(x)-h(x)=rx(l-x/N)-Ex

%捕捞情况下渔场鱼量满足的方程：x'(t)=F(x)

%满足F(x)=O的点x为方程的平衡点

%求方程的平衡点

symsrxNE;%定义符号变量

Fx=r*x*(l-x/N)-E*x;%创建符号表达式

x=solve(Fx,x)%求解F(x)=O(求根)

%得到两个平衡点，记为：

%xO=-N*(-r+E)/r,xl=O_,见物16(4)式

x0=x(2);

*1=(1);%符号变量x的结构类型成为V2*16加〉

%求¥6)的微分F'(x)

symsx;%定义符号变量x的结构类型为＜1*1§内＞

dF=diff(Fx,'x');%求导

dF=simple(dF)%简化符号表达式

%得F'(x)=r-2*r*x/N-E

%求F，(xO)并简化

dFxO=subs(dF,x,xO);%将x=xO代入符号表达式dF

dFxO=simple(dFxO)

%得F'(xO)=

%求F'(xl)

dFxl=subs(dF,x,xl)

%得F'(xl)=r-E

%若£<「，有F'(x0)<0,F'(xl)>0,故xO点稳定，xl点不稳定(根据平衡点稳

定性的准则);

%若E>i•，则结果正好相反。

%在渔场鱼量稳定在xO的前提下(E<r)，求E使持续产量h(xO)达到最大hm。

%通过分析(见教材p216图1),只需求xO*使f(x)达到最大，且hm=f(xO*)。

symsrxN

fx=r*x*(l-x/N);%fx=dF

df=diff(fX,'x');

xO=solve(df,x)

%得x0*=1/2*N,见P217(6)S

hm=subs(fxx,xO)

%得hm=l/4*r*N,见P217(7)式

%又由xO*=N(l-E/r),可得E*=r/2,见P217(8)S

%产量模型的结论是：

%将捕捞率控制在固有增长率的一半(E=r/2)时，能够获得最大的持续产量。

2.(验证、编程)种群的相互竞争P222-228

模型:

22

rx(1--u-cA_)11NiN

nxxrx(

12

其中，

X{t},X2(f)分别是甲乙两个种群的数量。

12

rik2是它们的固有增长率。

Ni,N2是它们的最大容量。

O，:单位数量乙(相对N2)消耗的供养甲的食物量为单位数量甲(相对N1)消耗的

供养甲的食物量的G7倍。对。2可作相应解释。

2.1(编程)稳定性分析p224~225

要求：

补充如下指出的程序段，然后运行该m文件，对照教材上的相应结果。

%7.3种群的相互竞争一一稳定性分析

%文件名：p224_225.ni

clear;clc;

%甲乙两个种群满足的增长方程：

%xr(t)=f(xl,x2)=rl*xl*(l-xl/Nl-kl*x2/N2)

%x2'(t)=g(xl,x2)=r2*x2*(l-k2*xl/Nl-x2/N2)

%求方程的平衡点，即解代数方程组(见P224的(5)式)

%f(xl,x2)=0

%g(xl,x2)=0

编写出该程序段。

[提示]

(1)使用符号表达式；

(2)用函数solve求解代数方程组，解放入[xl,x2];

(3)调整解(平衡点)的顺序放入P中(见下面注释所示)，P的结构类型为

<4x2sym>,P的第1列对应xl,第2列对应x2。

xlx2=[xl,x2]%显示结果

disp('力P

%调整位置后的4个平衡点：

%P(1,:)=P1(N1,0)

%P(2,:)=P2(0,N2)

%P(3,:)=P3(Nl*(-l+kl)/(-l+k2*kl),N2*(-l+k2)/(-l+k2*kl))

%P(4,:)=P4(0,0)

%平衡点位于第一象限才有意义，故要求P3:kl,k2同小于1,或同大于1。

%判断平衡点的稳定性参考教材p245的(18),(19)式。

symsxlx2;%重新定义

fXl=diff(f,'xT);fx2=diff(f,'x2');

gxl=diff(g,'xl');gx2=diff(g,'x2');

disp('');A=[fXlf2;gxl,gx2]%显示结果

p=subs(-(fXl+gx2),{xl,x2},{P(:,l),P(:,2)});%替换

p=simple(p);%简化符号表达式p

q=subs(det(A),{xl,x2},{P(:,l),P(:,2)});

q=simple(q);

disp(");[Ppq]%显示结果

%得到教材p225表1的前3列，经测算可得该表的第4歹U，即稳定条件。

★补充后的程序和运行结果(见[225]表1)：2341

%7.3种群的相互竞争一一稳定性分析

%文件名：p224_225.m

clear;clc;

%甲乙两个种群满足的增长方程：

%xl'(t)=f(xl,x2)=rl*xl*(l-xl/Nl-kl*x2/N2)

%x2'(t)=g(xl,x2)=r2*x2*(1-k2*x1/Nl-x2/N2)

%求方程的平衡点，即解代数方程组(见P224的(5)式)

%f(xl,x2)=0

%g(xl,x2)=0

%编写出该程序段。

symsxlx2rlr2N1N2klk2;

f=rl*xl*(l-xl/Nl-kl*x2/N2);

g=r2*x2*(l-k2*xl/Nl-x2/N2);

[xl,x2]=solve(f,g,xl,x2);

P=[Xl([234,l]),x2([2,3,4,l])l；

xlx2=[xl,x2J%显示结果

disp('');P

%调整位置后的4个平衡点：

%P(1,:)=P1(Nl,0)

%P(2,:)=P2(0,N2)

%P(3,:)=P3(Nl*(-l+kl)/(-l+k2*kl),N2*(-l+k2)/(-l+k2*kl))

%P(4,:)=P4(0,0)

%平衡点位于第一象限才有意义，故要求P3：kl,k2同小于1,或同大于L

%判断平衡点的稳定性参考教材p245的(18),(19)式。

symsxlx2;%重新定义

fxl=diff(f,'xr);fx2=diff(f,'x2');

gxl=diff(g,'xr);gx2=diff(g,'x2');

disp(A=[fxl,fx2;gxl,gx2]%显示结果

p=subs(-(fXl+gx2),{xl,x2},{P(:,l),P(:,2)});%替换p=simple(p);%简化符号表达式p

q=subs(det(A),{xl,x2},{P(:,l),P(:,2)});q=simple(q);

disp(");[Ppqj%显示结果

%得到教材p225表1的前3列，经测算可得该表的第4歹U，即稳定条件。

1C<iB*«>dVladM

£il«Jdit1t[de，”firrknU«lj

xlx2=

[0.0]

[NL.0]

[0,N2J

I(XL*(kl-1))/O:l*h2-1),(«2*(k2-D)/(kl*fc2-1)1

P;

(ra.o]

[o»N2]

I<KL*(kl-D)/Ql*k2-1),Ltt2»(k2-l»/(kl«fc2-L)J

I0,0]

(-rl*(xl/N]4(kl»«2)/X2-D-(rl»xD/Nl.-(ki»rl«Tl)/M2]

[-(k2*r2»x2)/Hl.-r2»(ic2JTl2*(k2»xl)/Nl-1)-Cr2»i2)/K2】

3ns=

[Nl,Url-r2*k2*r2„rl»r2*(k2-1)3

I。N2,r2-rt+kl*rl,rl*r2*(kl-1)J

[(XI*(kl-]))/ai*k2-D,(H24(k2-l))/(klMc2-I),"(rl4-x2-kl*rl-k2*r2)/(kl*k2-1),-Crl*r2*(kl-D*0c2-1»/Ckl*k2-1)J

[0,a-rl-T2.rl«r2]

A»|

2.2(验证、编程)计算与验证P227

微分方程组

工串)二丁「HV-。iN)

<12

XX

x(f)=rx(1—0--)

2221NN

Op0.5,O2-1.6,n=2.5,n-1.8,Ni=1.6,M=/

⑴(验证)当XI(0)R2(0)=0.1时，求微分方程的数值解，将解的数值分别画出xi(f)和

X2(。的曲线，它们同在一个图形窗口中。

程序:

%命令文件

ts=0:0.2:8;

xO=[O.l,O.l];

[t,x]=ode45('p227',ts,x0);

plot(t,x(:,l),t,x(:,2));%x(:l)是xl(t)的一组函数值,x(:,2)对应x2(t)

gridon;

axis([0802]);

text(2.4,l-55,'xl⑴')；

text(2.2,0.25；x2(t)');

%函数文件名：p227.m

functiony=p227(t,x)

kl=0.5;k2=1.6;rl=2.5;r2=1.8;Nl=1.6;N2=l;

y=[rl*x(l)*(Lx(l)/Nl-kl*x(2)/N2),r2*x(2)*(l-k2*x(l)/NLx(2)/N2)]';

☆⑴运行程序并给出结果(比较［227］图2(a))：

☆(2)(验证)将xi(O)=l,X2(0)=2代入(1)中的程序并运行。给出

结果(比较［227］图2(b))：

(3)(编程)在同一图形窗口内，画(1)和(2)的相轨线，相轨线是以xiQ)为横

坐标，为纵坐标所得到的一条曲线。具体要求弁照下图。

图中的两条“点线”直线，一条的两个端点为(0,1)和(1,0),另一条的两个端点

为(0,2)和(1.6,0)»

★⑶给出程序及其运行结果(比较［227］图2(c))：

%命令文件

ts=0:0,2:8;

xO=[O.l,O.l];

[t,x]=ode45(*p2271,ts,xO);

plot(x(:,l),x(:,2));%x(:l)是xl(t)的一组函数值，x(:,2)对应x2(t)

text(0.03,0.3；(0.1,0.1)');

holdon;

x0=[l,2];

[t,x]=ode45(,p227\ts,x0);

plot(x(:,l),x(：＞2));

text(1.02,1.9,(1,2),);

plot([0,l],[l,0],':r',[0,1.6],[2,0],':r');%商两条直线

gridon;

3.（编程）种群的相互依存一稳定性分析P228-229

模型:

XX、

x(t)-rx(1——+o——)

1—NiN

12

•XX

x(t)=rx(-1+0—T-)

222NN

12

其中,

X1（E）,X2（。分别是甲乙两个种群的数量。

门k2是它们的固有增长率。

Ni,N2是它们的最大容量。

0]：单位数量乙(相对N2)提供的供养甲的食物量为单位数量甲(相对N1)消耗的供

养甲的食物量的0,倍。对02可作相应解释。

要求：

☆修改题2.1的程序，求模型的平衡点及稳定性。给出程

序及其运行结果(比较[229]表1,注：只要最终结果):

clear;clc;

symsxlx2rlr2N1N2klk2;

f=rl*xl*(l-xl/Nl+kl*x2/N2);

g=r2*x2*(-l+k2*xl/Nl-x2/N2);

[xl,x2]=solve(f,g);

P=[xl([2,4,l,3]),x2([2,4,l,3])]；

symsxlx2;%重新定义

fxl=diff(f,'xr);fx2=diff(f,'x2');

gxl=diff(g,'xr);gx2=diff(g,'X2');

A=[fxl,fx2;gxl,gx2];

p=subs(-(fXl+gx2),{xl,x2},{P(:,l),P(:,2)});%替换

p=simple(p);%简化符号表达式p

q=subs(det(A),{xl,X2},{P(:,l),P(:,2)));

q=simple(q);

[Ppq]%显示结果

[HvE。hHFLiDKLrdnrHiL!?ans=

[■-4-r2-->■,T-'-It--LJ]

[■i--la.I--])■■1>,-HI--1-■!]4-k2*r2)/(k,-LJ.(H*r2*■-1),■?•LJ.■-LJ]

..：；.-it

(.::r.'rl-E-n.:v:■:•H]

4.（验证）食饵-捕食者模型p230~232

模型的方程:

x(t)-rx-axyx(0)二x

<o

y(0)二

y(0-—dy+bxy

要求：

设r=l,d=0.5,a=0.1,6=0.02,x0=25,y>=2。输入p231的程序并运行，结

果与教材p232的图1和图2比较。

☆给出2个M文件（见［231］）和程序运行后输出的图形

（比较［232］图1、2）:

函数M文件：

functionxdot=shier(t,x)r=l;d=0.5;a=0.1;b=0.02;

xdot=[(r-a*x(2)).*x(l);(-d+b*x(l)).*x(2)];

命令M文件：

ts=O: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(:,l),x(:,2)),grid,

（。图形:

相轨线j(r)图形:

5.(验证)差分形式的阻滞增长模型p236~242

阻滞增长模型用微分方程描述为：

也可用差分方程描述为：

y-即，*-

上式可简化为一阶非线性差分方程:

x=bx(1-x),k-0,1,2,

考察给定8和、值后，当A-8时；X&的收敛情况(实际上取上足够大就可以了)。

5.1数值解法和图解法p238~240

⑴取xo=O.2,分别取b=1.7,2.6,3.3,3.45,3.55,3.57,对方程

x-bx(l-x),k-0,1,2,

计算出的值，显示小的值。观察收敛与否。(结果与教材p238~239表

1比较)

下面是计算程序，在两处下划线的位置填入满足要求的内容。

%不同b值下方程x(k+l)=b*x(k)*(l-x(k)),k=0,1,2,…的计算结果

%文件名：p238table_l.m

clc;clearall;formatcompact;

b=[1.7,26,33,345,355,3-57];

x=zeros(100,length(b));

x0=0.2;%初值

；%此处填入一条正确语句

fork=l:99

;%此处填入一条正确语句

end

K=(81:100)';%将取81~100行

disp(num2str([NaN,b;K,x(K,:)],4));%取4位有效数字，NaN表示不是数值

★⑴对程序正确填空，然后运行。

填入的正确语句和程序的运行结果(对应不同的b值)见[238]表1:

x(l,:)=b*x0*(l-x0);

x(k+l,:)=b.*x(k,:).*(l-x(k,:));

!Co*AandVIadov口叵区)

FslfiSditRovktopYiedavHelp

NaN1.72.63.33.453.553.57

810.41180.61540.47940.44750.5060.4754

820.41180.615g0.82360.8530.88740.8903

830.41180.61540.47940.43260.35480.3486

840.41180.61540.82360.84680.21270.8106

850.41180.61540.47940.0.54050.548

860.41180.61540.82560.8530.881T0.8843

870.41180.61540.47940.43270.37030.3654

880.41180.61540.82360.84690.82780.8278

890.41180.6L590.47940.44740.5060.5089

900.41180.61540.82360.8530.83740.B922

910.41180.61540.47940.432?0.35480.3433

920.41180.61540.82360.24690.812T0.8049

930.41180.61540.47940.M740.54050.5607

940.41180.61540.82360.8530.88LT0.8793

950.41180.6L540.47940.93270.37030.3788

960.41180・6凤0.82360.3469S82780.84

970.41120.61540.47940.44740.5060.479?

980.41180.61540.82360.8530.88740.891

990.41180.6L540.47940.43270.35480.3466

1000.41180.61540.82360.84690.812T0.8085

A»I

(2)运行以下程序，观察显示的图形，与题(1)的数据对照，注意收敛的倍周

期。

%图解法，图1和图2

%文件名：p238fig_l_2.m

clear;clc;closeall;

f=@(x,b)b.*x.*(l-x);%定义f是函数的句柄，函数产瓦*(1-又)含两个变量x,b

b=[1.7,2,6,3.3,345,3-55,3-57J;

x=ones(101,length(b));

x(l,:)=0.2;

fork=l:100

x(k+l,:)=f(x(k,:),b);

end

forn=l:length(b)

立部「©3);%指定图形窗口figuren

subplot(l,2,l)；%开始迭代的图形

fpk>t(@(x)[f(x,b(n)),x],[O10l]);%x是自变量,画曲线y=bx(l-xA直线y=xaxis

square;holdon;

xO=x(l,n);y0=0;%商迭代轨迹线

fork=2:5

xl=x(k,n);yl=x(k,n);

plot([xO+i*yO,xO+i*yl,xl+i*yl]；「')；％实部为横坐标，虚部为纵坐标

xO=xl;yO=yl;

end

titled'图解法:开始迭代的图形(b='num2str(b(n)),)'D;

holdoff;

subplot(l,2,2);%最后迭代的图形

fplot(@(x)[f(x,b(n)),x],[0101]);

axissquare;holdon;

xy(l: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');

。加(「图解法：最后迭代的图形(b='num2str(b(n))')']);

holdoff;

end

☆(2)运行程序并给出结果(对应不同的b值)见［238］图1、2：

图解法：最后迭代的图形（b=1.7）

1

0.8

0.6

0.4

0.2

0

00.51

2。倍

图解法：开始迭代的图形(b=3.55)

图解法：开始迭代的图形(b=3.57)

5.2b值下的收敛图形p238~240

下面程序是在不同b值下的收敛图形。

%1)值下的收敛图形

%文件名：p212tablfigure.m

%方程(6)xk+l=b*xk*(l-xk),k=0,l,2,..

clear;clc;closeall;

b=[3.33.453.553.57];

x=0.2*ones(size(b));%初值x0=0.2

fork=l:100

x(k+l,:)=b,*x(k,:).*(l-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+l=b*xk*(l・xk),k=0,l,2„..

clear;clc;closeall;

b=2.5:0.0001:4;%参数b的间隔取值

x=0.2*ones(size(b));

xx=[];n=100000;

fork=l:n

x=b.*x.*(l-x);

ifk>=n-50xx=[xx;x];

endend

plot(b,xx「.b','markersize',5);%可改变值5,调整数据点图标的大小

titleC图4模型的收敛、分岔和混沌')；

xlabel(‘参数b*);

ylabel(*x(k)(k足够大)*)；gridon;

☆运行以上m文件。运行结果(比较［241］图

JFicuro1匕[Q|快|

Inx«rtTaalvopVmdtv¥・lp

-k♦「、*⑨,&/•3口司・ED

图4瘠型的收敛、分缜和泥沌

Q9

Q0

0

OQ.7

R6

必OO.5

室

A

》

音O.

Q.m

O.

2

1

Q

*I

5.42n倍周期图形p240~242

要求：

在上题的程序中，修改b值，使运行后显示以下图形:

★(1)单周期(14<3)：

•>FigureI

|FileEditVie»Insert£oolsDesktopWindowHelp

□诗⑤要□国I■国

图8模型的收敛、分岔和混沌

07

06

K(

绘

e叱

茗

x

★⑵2倍周期(3<b<3.449)：

SB®

£il«£ditVitwloolsQ«sktop£indov{{tip

口「H2k♦、，C黝S/•@□国■0

图8模型的收敛、分岔和混沌

0.9

o6

.0a8

o75

(G7

长

卷5

叫o

*

).e6

w

xo3.5

o

5

0,

45

o

3053.13.153.23.253.33.353.43.453.5

参数b

★(3)22倍周期(3.449<b<3.544)：

Figure131130

FiliEditViewInsertToolEDesktopWindowHelp

旨a久踮⑤□国i■国

图8模型的收敛、分岔和混沌

09

0E

(

K7

徐00.

或

*

)6

xg0

5

3443.463483.53.523.643.563.58

参数b

★(4)23倍周期(3.544<b<3.564)：

FigureI|Z州0X|

FiliEditViewInsertTOOlEDesktopWindowHelp

bT洪kN踮⑥*£▼二［口口出国

R

m7

o.

B。

嗖

二6

to

5

图8模型的收敛、分岔和混沌

1

09

08

5.5(编程)求2n倍周期的b值区间P240-241

求2人口倍周期收敛的b的上限的程序如下：

functionbmax=fUn(bn_l,n)

%求2八n倍周期收敛的b的上限。

3.543.553.56

3.55

%bn_l是2A(口-1)倍周期收敛的b的上限(或2An倍周期收敛的b的下限)。

c=bn_l;d=3.57;%2An倍周期时收敛b>bn_l,b<3.57

y=zeros(l,2*2An);%行向量，用于存放xk最后的2x2An个值

whiled-c>le-12%使区间(c,d)足够小

b=(d+c)/2;x=0.2;%b取区间的中点

fori=l:100000

x=b*x*(l-x);

end

y(l)=x;%取最后2*2An个xk

fork=l:2A(n+l)-l

y(k+l)=b*y(k)*(l-y(k));

end

ifnorm(y(l:2An)-y(2An+l：2A(n+l)))<le-12%范数，比较

c=b;%满足2人口倍周期收敛的b给区间(c,d)的下限c

else

d=b;%不满足2An倍周期收敛的b给区间(c,d)的上限dend

end

bmax=c;%2人n倍周期收敛的b的上限近似

要求：

编写程序，调用以上函数文件，求单倍周期、2倍周期、......、25倍周期收

敛的b的上限近似值，输出要求有10位有效数字。运行结果输出形式如下：

提示：

可使用num2str函数。

用下面的程序结构，会使运行速度加快。

functionmain()

自己编写的程序：

将上面的函数文件复制到此处。

★运行的程序及输出结果：

functionmain()clc;

n=5;%求单(赋0)、2倍(赋1)...................2人口倍周期收敛的b的上、卜限

B=ones(n+2,1);fori=0:n

B(i+2)=fUn(B(i+l),i)；

end%单周期收敛时,B(l)<b<B(2)；2倍时,B(2)<b<B(3);...

disp(num2str([(-l:n),,BJ,10));%10位有效数字

functionbmax=fUn(bn_l,n)

%求2人n倍周期收敛的b的上限。

%bn_l是2人(口-1)倍周期收敛的b的上限(或2八n倍周期收敛的b的下限)。

c=bn_l;d=3.57;%2八11倍周期时收敛b>bn_l,b<3.57

y=zeros(l,2*2An);%行向量，用于存放xk最后的2*2An个值

whiled-c>le-12%使区间(c,4)足够小

b=(d+c)/2;x=0.2;%b取区间的中点

fori=l:100000

x=b*x*(l-x);end

y(l)=x;%取最后2*2An个xk

fork=l:2A(n+l)-l

y(k+l)=b*y(k)*(l-y(k));end

ifnorm(y(l:2An)-y(2An+l:2A(n+l)))<le-12%范数，比较

c=b;%满足2An倍周期收敛的b给区间(c,d)的下限c

else

d=b;%不满足2An倍周期收敛的b给区间(c,d)的上限d

endend

bmax=c;%2An倍周期收敛的b的上限近似

JCowandR］叵］区］

EditDebug»勺

-11

02.999769097

13.449355967

23.544021471

33.564357879

43.568739993

53.569679465

AA>|

注：vpa的用法。

附1:实验提示第1题

提示：符号简化函数5而03,变量替换函数SUb

符号简化函数simple的格式：

simple(S)

对符号表达式S尝试多种不同的算法简化以显示S表达式的长度最短的简

化形式。

变量替换函数sub的格式：

subs(S,OLD,NEW)

将符号表达式S中的OLD变量替换为NEW变量。

附2：第7章稔定性模型

[215]7.I捕鱼业的持续收获

第7章

稳定性模型

虽然动态过程的变化规律一般要用微分方程建立的动态模型来描述，但是对于某些实际问

题，建模的主要目的并不是要寻求动态过程每个瞬时的性态，而是研究某种意义下稳定状态的特

征，特别是当时间充分长以后动态过程的变化趋势.比如，在什么条件下描述过程的变量会越来

越接近某些确定的数值，在什么情况下又会越来越远离这些数值而导致过程不稳定，为了分析这

种稳定与不稳定的规律，常常不需要求解微分方程（并且我们将会看到，即使对于不太复杂的

方程，解析解也不是总能得到的八而可以利用微分方程稳定性理论，直接研究平衡状态的稳定性

就行了.

在下面建模过程中，大多数情况下我们将直接引用微分方程稳定性理论的结果.不熟悉这

方面内容的读者可以参阅7门节.

差分方程的稳定性与微分方程有很多相似之处.也有明显的差别，7.6和7.7节将介绍这

方面的内容.

7.1捕鱼业的持续收获

可持续发展是一项基本国策，对于像渔业、林业这样的再生资源，一定要注意适度开发，

不能为了一时的高产去“竭泽而渔”，应该在持续稳产的前提下追求产量或效益的最优化.

考察一个渔场，其中的鱼量在天然环境下按一定规律增长，如果捕捞量恰好等于增长量，

那么渔场鱼量将保持不变，这个捕捞量就可以持续下去•本节要建立在捕捞情况下油场鱼量遵从

的方程，分析鱼量稳定的条件，并且在检定的前提下讨论如何控制捕捞使持续产量或经济效益达

到最大•最后研究所谓捕捞过度的问就门叫

产■模型记时刻士渔场中鱼量为4£）.关于的自然增长和人工捕捞作如下假设：

L在无捕捞条件下工的增长服从logistic规律（地5.6节的阻滞增长模型），即

=/（4）=rx

r是固有增长率小是环境容许的最大鱼量，用人式）表示单位时间的增长量.

久单位时间的捕捞量（即产量）与渔场鱼量工（E）成正比，比例常数£表示单位时间捕捞

率，又称为捕捞强度，可以用比如捕鱼网眼的大小或出海渔船数量来控制其大小.于是单位时间的

捕捞量为

h（x）-Ex

根据以上假设并记

=f（x）-h（x）

得到捕捞情况下渔场鱼最满足的方程

金（力=尸的=rx

我们并不需要解方程（3）以得到工《）的动态变化过程，只希望知道渔场的稳定鱼量和保持稳定

的条件，即时间士足够长以后渔场鱼量武士）的趋向，并由此确定最大持续产量，为此可以直接求

方程（3）的平衡点并分析其稳定性.

得到两个平衡点

与=1-，物二0

不难算出所以若

有尸Yw）《。，尸，（IHJ）>0,故0点稳定，跖点不稳定（判断平衡点稳定性的准则见7.7

节）；若E>r,则结果正好相反.

E是捕捞率，「是最大的增长率，上述

分析表明，只要捕捞适度（总</）,就可使渔，及巾，

场鱼量稳定在入，从而获得持续产展/_

从与）二而当捕捞过度时（£>「）・渔场”二二二2才

鱼量将趋向如二。，当然谈不上获得持续产"一一方才丁

进一步讨论渔场鱼量稳定在外的前提下，如何控制捕捞强度E使持续产量最大的问题.用图

解法可以非常简单地得到­?

*T

图】最大持续产量的图解法

结果.

根据(1),(2)式作抛物线,=/(/)和直线7=乂#)=氐，如图L注意到y=人工)在原点的切线为

y=所以在条件(5)下，=&必与7=汽动有交点尸，P的横坐标就是稳定平衡点

根据假设2,尸点的纵坐标力为稳定条件下单位时间的持续产量.由图1立刻知道,当y=&

与在抛物线顶点P■相交时可获得最大的持续产量“此时的稳定平衡点为

且单位时间的最大持续产量为

院二不⑺

而由(4)式不难算出保持渔场鱼量稳定在工；的捕捞率为

=y⑻

综上所述，产量模型的结论是将捕捞率控制在固有增长率r