%通过分析（见教材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.5445.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（
升级会员 