1、实验08 稳定性模型（4学时）（第7章 稳定性模型）1.（验证）捕鱼业的持续收获 产量模型p215219产量模型：其中，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为
2、方程的平衡点%求方程的平衡点syms r x N E; %定义符号变量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的结构类型成为%求F(x)的微分F(x)syms x; %定义符号变量x的结构类型为dF=diff(Fx,x); %求导dF=simple(dF) %简化符号表达式%得 F(x)= %求F(x0)并简化dFx0=subs(dF,x,x0); %将x=x0代入符号表达式dFdFx0=simple(dFx0)
3、%得 F (x0)= %求F (x1)dFx1=subs(dF,x,x1)%得 F (x1)= %若 Er，有F(x0)0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若 Er，则结果正好相反。%在渔场鱼量稳定在x0的前提下（Er），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且hm=f(x0*)。syms r x Nfx=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)式%又由
4、 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)%捕捞条件下单位时间的捕捞量：h(x)=Ex%F(x)=f(x)-h(x)=rx(1-x/N)-Ex%捕捞情况下渔场鱼量满足的方程：x(t)=F(x)%满足F(x)=0的点x为方程的平衡点%求方程的平
5、衡点syms r x N E; %定义符号变量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的结构类型成为%求F(x)的微分F(x)syms x; %定义符号变量x的结构类型为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
6、=simple(dFx0) %得 F (x0)= -r+E %求F (x1)dFx1=subs(dF,x,x1)%得 F (x1)= r-E %若 Er，有F(x0)0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若 Er，则结果正好相反。%在渔场鱼量稳定在x0的前提下（Er），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p216图1），只需求x0*使f(x)达到最大，且hm=f(x0*)。syms r x Nfx=r*x*(1-x/N); %fx=dFdf=diff(fx,x);x0=solve(df,x)%得 x0*= 1/2*N ，见P217(6)式hm=subs
7、(fx,x,x0)%得 hm= 1/4*r*N ，见P217(7)式%又由 x0*=N(1-E/r)，可得 E*= r/2 ，见P217(8)式%产量模型的结论是：%将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。2.（验证、编程）种群的相互竞争P222228模型：其中，x1(t), x2(t)分别是甲乙两个种群的数量。r1, r2是它们的固有增长率。N1, N2是它们的最大容量。1：单位数量乙（相对N2）消耗的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的1倍。对2可作相应解释。2.1（编程）稳定性分析p224225要求：补充如下指出的程序段，然后运行该
8、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的结构类型为，P的第1列对应
9、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)式。syms x1 x2; %重新定义fx1=diff(f,x1); fx2=diff(f,x2);gx1=diff(g,x1); gx2=diff(g,
10、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( ); P p q %显示结果%得到教材p225表1的前3列，经测算可得该表的第4列，即稳定条件。 补充后的程序和运行结果（见225表1）：2341%7.3 种群的相互竞争稳定性分析%文件名：p224_225.mclear; clc;%甲乙两个种群满足的增长方程：% x1(t)=f(x1,x2)=
11、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%编写出该程序段。syms x1 x2 r1 r2 N1 N2 k1 k2;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
12、(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)式。syms x1 x2; %重新定义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)
13、,x1,x2,P(:,1),P(:,2); %替换p=simple(p);%简化符号表达式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);disp( ); P p q %显示结果%得到教材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
14、,x(:,1),t,x(:,2);%x(:1)是x1(t)的一组函数值，x(:,2)对应x2(t)grid on;axis(0 8 0 2);text(2.4,1.55,x1(t);text(2.2,0.25,x2(t);%函数文件名：p227.mfunction y=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
15、(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,
16、0.3,( 0.1, 0.1 );hold on;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);%画两条直线grid on;3.（编程）种群的相互依存稳定性分析P228229模型：其中，x1(t), x2(t)分别是甲乙两个种群的数量。r1, r2是它们的固有增长率。N1, N2是它们的最大容量。1：单位数量乙（相对N2）提供的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的1倍。对2可作相应解释。要求： 修改题2.1的程序，求
17、模型的平衡点及稳定性。给出程序及其运行结果（比较229表1，注：只要最终结果）：clear; clc;syms x1 x2 r1 r2 N1 N2 k1 k2;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);syms x1 x2; %重新定义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
18、(:,1),P(:,2); %替换p=simple(p);%简化符号表达式p q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);P p q %显示结果4.（验证）食饵-捕食者模型p230232模型的方程：要求：设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文件：function xdot=shier(t,x)r=1; d=0.5; a=0.1; b=0.02;xdot=(r-a*x(2
19、).*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.（验证）差分形式的阻滞增长模型p236242阻滞增长模型用微分方程描述为：也可用差分方程描述为：上式可简化为一阶非线性差分方程：考察给定b和x0值后，当k 时，xk的收敛情况（实际上取k足够大就可以了）。5.1 数值解
20、法和图解法p238240(1) 取x0=0.2，分别取b = 1.7, 2.6, 3.3, 3.45, 3.55, 3.57，对方程计算出x1 x100的值，显示x81 x100的值。观察收敛与否。（结果与教材p238239表1比较）下面是计算程序，在两处下划线的位置填入满足要求的内容。%不同b值下方程 x(k+1)=b*x(k)*(1-x(k), k=0,1,2,.的计算结果%文件名：p238table_1.mclc; clear all; format compact;b=1.7,2.6,3.3,3.45,3.55,3.57;x=zeros(100,length(b);x0=0.2; %初
21、值 ; %此处填入一条正确语句for k=1:99 ; %此处填入一条正确语句endK=(81:100); %将取81100行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; close all;f
22、=(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;for k=1:100 x(k+1,:)=f(x(k,:),b);endfor n=1:length(b) figure(n);%指定图形窗口figure n subplot(1,2,1);%开始迭代的图形 fplot(x)f(x,b(n),x,0 1 0 1);%x是自变量，画曲线y=bx(1-x)和直线y=x axis square; hold on; x0=x(1,n);
23、 y0=0; %画迭代轨迹线 for k=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) ）); hold off; subplot(1,2,2); %最后迭代的图形 fplot(x)f(x,b(n),x,0 1 0 1); axis square; hold on; xy(1:2:41)=x(81:101,n)+i*x(81:101,n);%尽量不用循环 xy(2:2:40)=x(81:
24、100,n)+i*x(82:101,n); plot(xy,r); title(图解法：最后迭代的图形（b= num2str(b(n) ）); hold off;end(2) 运行程序并给出结果（对应不同的b值）见238图1、2：20倍20倍21倍22倍23倍混沌5.2 b值下的收敛图形p238240下面程序是在不同b值下的收敛图形。%b值下的收敛图形%文件名：p212tab1figure.m%方程（6） xk+1=b*xk*(1-xk), k=0,1,2,.clear; clc; close all;b=3.3 3.45 3.55 3.57;x=0.2*ones(size(b); %初值x0
25、=0.2for k=1:100 x(k+1,:)=b.*x(k,:).*(1-x(k,:);endplot(b,x(82:101,:),.r ,markersize,5); %可改变值5，调整数据点图标的大小xlabel(b);ylabel(x (k);grid on要求： 运行以上程序。 在运行结果的图形中，从对应的b值上的“点数”判断倍周期收敛。提示：放大图形。 程序的运行结果（见238表1）：5.3 收敛、分岔和混沌p240242画出教材p241图4 模型的收敛、分岔和混沌。程序的m文件如下：%图4 模型（6）的收敛、分岔和混沌%文件名：p241fig4.m%模型：xk+1=b*xk*(
26、1-xk), k=0,1,2,.clear; clc; close all;b=2.5:0.0001:4; %参数b的间隔取值x=0.2*ones(size(b);xx=; n=100000;for k=1:n x=b.*x.*(1-x); if k=n-50 xx=xx;x; endendplot(b,xx,.b,markersize,5); %可改变值5，调整数据点图标的大小title(图4 模型的收敛、分岔和混沌);xlabel(参数 b);ylabel(x(k) （k足够大）);grid on; 运行以上m文件。运行结果（比较241图4）：5.4 2n倍周期图形p240242要求：在上
27、题的程序中，修改b值，使运行后显示以下图形：(1) 单周期（1b3）：(2) 2倍周期（3b3.449）：(3) 22倍周期（3.449b3.544）：(4) 23倍周期（3.544bbn_1，b1e-12 %使区间(c, d)足够小 b=(d+c)/2; x=0.2; %b取区间的中点 for i=1:100000 x=b*x*(1-x); end y(1)=x; %取最后22n个xk for k=1:2(n+1)-1 y(k+1)=b*y(k)*(1-y(k); end if norm(y(1:2n)-y(2n+1:2(n+1)1e-12 %范数，比较 c=b;%满足2n倍周期收敛的b给区
28、间(c, d)的下限c else d=b; %不满足2n倍周期收敛的b给区间(c, d)的上限d endendbmax=c; % 2n倍周期收敛的b的上限近似要求：编写程序，调用以上函数文件，求单倍周期、2倍周期、25倍周期收敛的b的上限近似值，输出要求有10位有效数字。运行结果输出形式如下：提示：可使用num2str函数。用下面的程序结构，会使运行速度加快。function main()自己编写的程序；将上面的函数文件复制到此处。 运行的程序及输出结果：function main()clc;n=5;%求单（赋0）、2倍（赋1）、2n倍周期收敛的b的上、下限B=ones(n+2, 1);for
29、 i=0:n B(i+2)=fun(B(i+1),i);end%单周期收敛时，B(1)bB(2)；2倍时，B(2)bbn_1，b1e-12 %使区间(c, d)足够小 b=(d+c)/2; x=0.2; %b取区间的中点 for i=1:100000 x=b*x*(1-x); end y(1)=x; %取最后22n个xk for k=1:2(n+1)-1 y(k+1)=b*y(k)*(1-y(k); end if norm(y(1:2n)-y(2n+1:2(n+1)1e-12 %范数，比较 c=b;%满足2n倍周期收敛的b给区间(c, d)的下限c else d=b; %不满足2n倍周期收敛的
30、b给区间(c, d)的上限d endendbmax=c; % 2n倍周期收敛的b的上限近似注：vpa的用法。附1：实验提示第1题提示：符号简化函数simple，变量替换函数sub符号简化函数simple的格式：simple(S)对符号表达式S尝试多种不同的算法简化，以显示S表达式的长度最短的简化形式。变量替换函数sub的格式：subs(S,OLD,NEW)将符号表达式S中的OLD变量替换为NEW变量。57附2：第7章 稳定性模型2157.1 捕鱼业的持续收获219*本节完*2227.3 种群的相互竞争225 题2.1答案227 题2.2答案2287.4 种群的相互依存229 题3答案2307.5 食饵-捕食者型231 题4程序232 题4答案2367.6 差分形式的阻滞增长模型238 题5.1、5.2答案241 题5.3答案242*本节完*
