实验08讲评参考答案稳定性模型4学时资料Word文档格式.docx
- 文档编号:7000412
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:57
- 大小:2.81MB
实验08讲评参考答案稳定性模型4学时资料Word文档格式.docx
《实验08讲评参考答案稳定性模型4学时资料Word文档格式.docx》由会员分享,可在线阅读,更多相关《实验08讲评参考答案稳定性模型4学时资料Word文档格式.docx(57页珍藏版)》请在冰点文库上搜索。
1sym>
%求F(x)的微分F'
(x)
symsx;
%定义符号变量x的结构类型为<
1×
dF=diff(Fx,'
);
%求导
dF=simple(dF)%简化符号表达式
%得F'
(x)=
%求F'
(x0)并简化
dFx0=subs(dF,x,x0);
%将x=x0代入符号表达式dF
dFx0=simple(dFx0)
(x0)=
(x1)
dFx1=subs(dF,x,x1)
(x1)=
%若E<
r,有F'
(x0)<
0,F'
(x1)>
0,故x0点稳定,x1点不稳定(根据平衡点稳定性的准则);
%若E>
r,则结果正好相反。
%在渔场鱼量稳定在x0的前提下(E<
r),求E使持续产量h(x0)达到最大hm。
%通过分析(见教材p216图1),只需求x0*使f(x)达到最大,且hm=f(x0*)。
symsrxN
fx=r*x*(1-x/N);
%fx=dF
df=diff(fx,'
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文件(见[215~217]):
%x0=-N*(-r+E)/r,x1=0,见P216(4)式
(x)=r-2*r*x/N-E
(x0)=-r+E
(x1)=r-E
%得x0*=1/2*N,见P217(6)式
%得hm=1/4*r*N,见P217(7)式
%又由x0*=N(1-E/r),可得E*=r/2,见P217(8)式
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
%甲乙两个种群满足的增长方程:
%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,'
gx2=diff(g,'
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(:
q=simple(q);
[Ppq]%显示结果
%得到教材p225表1的前3列,经测算可得该表的第4列,即稳定条件。
★补充后的程序和运行结果(见[225]表1):
2341
%编写出该程序段。
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])];
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)):
plot(x(:
1),x(:
text(0.03,0.3,'
(0.1,0.1)'
holdon;
x0=[1,2];
text(1.02,1.9,'
(1,2)'
plot([0,1],[1,0],'
:
r'
[0,1.6],[2,0],'
%画两条直线
3.(编程)种群的相互依存——稳定性分析P228~229
单位数量乙(相对N2)提供的供养甲的食物量为单位数量甲(相对N1)消耗的供养甲的食物量的σ1倍。
☆修改题2.1的程序,求模型的平衡点及稳定性。
给出程序及其运行结果(比较[229]表1,注:
只要最终结果):
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])];
A=[fx1,fx2;
gx1,gx2];
[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];
shier'
[t,x],
plot(t,x),grid,gtext('
x(t)'
),gtext('
y(t)'
),%运行中在图上标注
pause,
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
closeall;
f=@(x,b)b.*x.*(1-x);
%定义f是函数的句柄,函数b*x*(1-x)含两个变量x,b
x=ones(101,length(b));
)=0.2;
100
x(k+1,:
)=f(x(k,:
),b);
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],'
%实部为横坐标,虚部为纵坐标
x0=x1;
y0=y1;
end
title(['
图解法:
开始迭代的图形(b='
num2str(b(n))'
)'
]);
holdoff;
subplot(1,2,2);
%最后迭代的图形
xy(1:
2:
41)=x(81:
101,n)+i*x(81:
101,n);
%尽量不用循环
xy(2:
40)=x(81:
100,n)+i*x(82:
plot(xy,'
最后迭代的图形(b='
☆
(2)运行程序并给出结果(对应不同的b值)见[238]图1、2:
20倍
21倍
22倍
23倍
混沌
5.2b值下的收敛图形p238~240
下面程序是在不同b值下的收敛图形。
%b值下的收敛图形
p212tab1figure.m
%方程(6)xk+1=b*xk*(1-xk),k=0,1,2,...
b=[3.33.453.553.57];
x=0.2*ones(size(b));
%初值x0=0.2
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,...
b=2.5:
0.0001:
4;
%参数b的间隔取值
xx=[];
n=100000;
n
x=b.*x.*(1-x);
ifk>
=n-50
xx=[xx;
x];
plot(b,xx,'
.b'
'
title('
图4模型的收敛、分岔和混沌'
参数b'
x(k)(k足够大)'
☆运行以上m文件。
运行结果(比较[241]图4):
5.42n倍周期图形p240~242
在上题的程序中,修改b值,使运行后显示以下图形:
★
(1)单周期(1<
b<
3):
★
(2)2倍周期(3<
3.449):
★(3)22倍周期(3.449<
3.544):
★(4)23倍周期(3.544<
3.564):
5.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);
y
(1)=x;
%取最后2×
2^n个xk
fork=1:
2^(n+1)-1
y(k+1)=b*y(k)*(1-y(k));
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
bmax=c;
%2^n倍周期收敛的b的上限近似
编写程序,调用以上函数文件,求单倍周期、2倍周期、……、25倍周期收敛的b的上限近似值,输出要求有10位有效数字。
运行结果输出形式如下:
可使用num2str函数。
用下面的程序结构,会使运行速度加快。
functionmain()
自己编写的程序;
将上面的函数文件复制到此处。
★运行的程序及输出结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验 08 讲评 参考答案 稳定性 模型 学时 资料
![提示](https://static.bingdoc.com/images/bang_tan.gif)