实验08 稳定性模型（4学时）（第7章 稳定性模型）1.（验证）捕鱼业的持续收获 ——产量模型p215~219产量模型：()()1x x t F x rx Ex N ⎛⎫==-- ⎪⎝⎭其中，x(t)为t 时刻渔场中的鱼量。
r 是固有增长率。
N 是环境容许的最大鱼量。
E 是捕捞强度，即单位时间捕捞率。
要求：运行下面的m 文件，并把相应结果填空，即填入“_________”。
符号简化函数simple，变量替换函数sub的用法见提示。
★给出填空后的M文件（见[215~217]）：2.（验证、编程）种群的相互竞争P222~228 模型：⎪⎪⎩⎪⎪⎨⎧--=--=••)1()()1()(2211222222111111N x N x x r t x N x N x x r t x σσ 其中，x 1(t ), x 2(t )分别是甲乙两个种群的数量。
r 1, r 2是它们的固有增长率。
N 1, N 2是它们的最大容量。
σ1：单位数量乙（相对N 2）消耗的供养甲的食物量为单位数量甲（相对N 1）消耗的供养甲的食物量的σ1倍。
对σ2可作相应解释。
2.1（编程）稳定性分析p224~225要求：补充如下指出的程序段，然后运行该m 文件，对照教材上的相应结果。
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(' '); [P p q] %显示结果%得到教材p225表1的前3列，经测算可得该表的第4列，即稳定条件。
2.2（验证、编程）计算与验证p227微分方程组1211111212222212121212()(1)()(1)0.5, 1.6, 2.5, 1.8, 1.6, 1x xx t r xN Nx xx t r xN Nr r N Nσσσσ••⎧=--⎪⎪⎨⎪=--⎪⎩======(1)（验证）当x1(0)=x2(0)=0.1时，求微分方程的数值解，将解的数值分别画出x1(t)和x2(t)的曲线，它们同在一个图形窗口中。
程序：☆(2) （验证）将x1(0)=1, x2(0)=2代入(1)中的程序并运行。
给出(3)（编程）在同一图形窗口内，画(1)和(2)的相轨线，相轨线是以x1(t)为横x2(t)为纵坐标所得到的一条曲线。
具体要求参照下图。
坐标，3.（编程）种群的相互依存——稳定性分析P228~229模型：1211111212222212()(1)()(1)x x x t r x N N x x x t r x N N σσ••⎧=-+⎪⎪⎨⎪=-+-⎪⎩其中，x1(t), x2(t)分别是甲乙两个种群的数量。
r1, r2是它们的固有增长率。
N1, N2是它们的最大容量。
σ1：单位数量乙（相对N2）提供的供养甲的食物量为单位数量甲（相对N1）消耗的供养甲的食物量的σ1倍。
对σ2可作相应解释。
要求：☆修改题2.1的程序，求模型的平衡点及稳定性。
给出程序及其运行结果（比较[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(:,1),P(:,2)}); %替换p=simple(p);%简化符号表达式pq=subs(det(A),{x1,x2},{P(:,1),P(:,2)});q=simple(q);[P p q] %显示结果4.（验证）食饵-捕食者模型p230~232模型的方程：00()(0)()(0)x t rx axy x x y t dy bxy y y =-=⎧⎨=-+=⎩要求：设r =1，d =0.5，a =0.1，b =0.02，x 0=25，y 0=2。
输入p231的程序并运行，结果与教材p232的图1和图2比较。
☆ 给出2个M 文件（见[231]）和程序运行后输出的图形（比较[232]图1、2）：函数M 文件：命令M 文件：x (t ), y (t )图形：相轨线y (x )图形：5.（验证）差分形式的阻滞增长模型p236~242阻滞增长模型用微分方程描述为：)1()(Nx rx t x -=•也可用差分方程描述为：1(1),0,1,2,kk k k y y y r y k N+-=-=上式可简化为一阶非线性差分方程：1(1),0,1,2,k k k x b x x k +=-=考察给定b 和x 0值后，当k → ∞时，x k 的收敛情况（实际上取k 足够大就可以了）。
5.1 数值解法和图解法p238~240(1) 取x 0=0.2，分别取b = 1.7, 2.6, 3.3, 3.45, 3.55, 3.57，对方程1(1),0,1,2,k k k x b x x k +=-=计算出x 1 ~ x 100的值，显示x 81 ~ x 100的值。
观察收敛与否。
（结果与教材p238~239表1比较）下面是计算程序，在两处下划线的位置填入满足要求的内容。
★(1) 对程序正确填空，然后运行。
填入的正确语句和程序的运行结果（对应不同的b 值）见[238]表1：(2) 运行以下程序，观察显示的图形，与题(1)的数据对照，注意收敛的倍周期。
%图解法，图1和图2%文件名：p238fig_1_2.mclear; clc; close all;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;for k=1:100x(k+1,:)=f(x(k,:),b);endfor n=1:length(b)figure(n);%指定图形窗口figure n5.2 b值下的收敛图形p238~240 下面程序是在不同b值下的收敛图形。
plot(b,x(82:101,:),'.r' ,'markersize',5); %可改变值5，调整数据点图标的大小xlabel('b');ylabel('x (k)');grid on要求：①运行以上程序。
②在运行结果的图形中，从对应的b值上的“点数”判断倍周期收敛。
提示：放大图形。
程序的运行结果（见[238]表5.3 收敛、分岔和混沌p240~242画出教材p241图4 模型的收敛、分岔和混沌。
程序的m文件如下：☆运行以上m文件。
运行结果（比较[241]图4）：5.4 2n倍周期图形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的上限的程序如下：function bmax=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个值while d-c>1e-12 %使区间(c, d)足够小b=(d+c)/2; x=0.2; %b取区间的中点for i=1:100000x=b*x*(1-x);endy(1)=x; %取最后2×2^n个xkfor k=1:2^(n+1)-1y(k+1)=b*y(k)*(1-y(k));endif norm(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)的上限d endendbmax=c; % 2^n倍周期收敛的b的上限近似要求：编写程序，调用以上函数文件，求单倍周期、2倍周期、……、25倍周期收敛的b的上限近似值，输出要求有10位有效数字。
运行结果输出形式如下：提示：可使用num2str函数。
用下面的程序结构，会使运行速度加快。
function main()自己编写的程序；将上面的函数文件复制到此处。
function main()clc;注：vpa的用法。
附1：实验提示第1题提示：符号简化函数simple，变量替换函数sub符号简化函数simple的格式：simple(S)对符号表达式S尝试多种不同的算法简化，以显示S表达式的长度最短的简化形式。
变量替换函数sub的格式：subs(S,OLD,NEW)将符号表达式S中的OLD变量替换为NEW变量。
附2：第7章稳定性模型[215]7.1 捕鱼业的持续收获[219]****本节完****[222]7.3 种群的相互竞争[225] 题2.1答案[227] 题2.2答案[228]7.4 种群的相互依存[229] 题3答案[230]7.5 食饵-捕食者型。