河北大学数学模型实验 实验报告一、实验目的 1.学会编写程序段。2.能根据m文件的结果进行分析。3.根据图像进行比较和分析。二、实验要求 8-1捕鱼业的持续收获运行下面的m文件，并把相应结果填空，即填入“_”。clear;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为方程的平衡点 %求方程的平衡点 syms r x N E; %定义符号变量 Fx=r*x*(1-x/N)-E*x; %创建符号表达式 x=solve(Fx,x) %求解F(x)=0（求根） %得到两个平衡点，记为： % x0=_ , x1=_ 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代入符号表达式dF dFx0=simple(dFx0) %得 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 N fx=r*x*(1-x/N); df=diff(fx,x); x0=solve(df,x) %得 x0*=_ hm=subs(fx,x,x0) %得 hm=_ %又由 x0*=N(1-E/r)，可得 E*=_ %产量模型的结论是： %将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。 8-2种群的相互竞争（1）补充如下指出的程序段，然后运行该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) %求方程的平衡点，即解代数方程组 % f(x1,x2)=0 g(x1,x2)=08-3种群的相互竞争（2）补充如下指出的程序段，然后运行该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) %求方程的平衡点，即解代数方程组 % f(x1,x2)=0 g(x1,x2)=0 编写出该程序段。三、实验内容8-1捕鱼业的持续收获产量模型%文件名：p178.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为方程的平衡点%求方程的平衡点syms r x N E; %定义符号变量Fx=r*x*(1-x/N)-E*x; %创建符号表达式x=solve(Fx,x) %求解F(x)=0（求根） %得到两个平衡点，记为：% x0=N(1-x/N) , x1=0x0=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)%得 F(x0)=E-r%求F(x1)=r-EdFx1=subs(dF,x,x1)%得 F(x1)=r-E%若 Er，有F(x0)0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；%若 Er，则结果正好相反。%在渔场鱼量稳定在x0的前提下（Er），求E使持续产量h(x0)达到最大hm。%通过分析（见教材p178图1），只需求x0*使f(x)达到最大，且hm=f(x0*)。syms r x Nfx=r*x*(1-x/N);a) df=diff(fx,x);x0=solve(df,x)%得 x0*=N/2hm=subs(fx,x,x0)%得 hm=(r*N)/4%又由 x0*=N(1-E/r)，可得 E*=r/2%产量模型的结论是：%将捕捞率控制在固有增长率的一半（E=r/2）时，能够获得最大的持续产量。提示符号简化函数simple的格式：simple(S)对符号表达式S尝试多种不同的算法简化，以显示S表达式的长度最短的简化形式。变量替换函数sub的格式：Subs(S,OLD,NEW)将符号表达式S中的OLD变量替换为NEW变量。8-2种群的相互竞争%文件名：p186.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)%求方程的平衡点，即解代数方程组% f(x1,x2)=0% g(x1,x2)=0%得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。%判断平衡点的稳定性（参考教材p200）fx1=diff(f,x1);fx2=diff(f,x2);gx1=diff(g,x1);gx2=diff(g,x2);A=fx1,fx2;gx1,gx2syms x1 x2;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);P p q%得到教材p186表2的前3列，经测算可得该表的第4列，即稳定条件。8-3种群的相互竞争（2）求微分方程组的数值解，分别画出教材p189中的图5-1，图5-2，图5-3。有关数据参见教材p188中“计算与验证”。 提示（1）求微分方程组的数值解可参考教材p139的程序。（2）在figure(1)中画图5-1，在figure(2)中画图5-2，在figure(3)中画图5-3。在程序中，figure(图形编号)用于定位对应图形。（3）使用text(x,y,标识文本)，坐标点(x,y)在“标识文本”的左边，调整(x,y)值，使“标识文本”放在图中的适当位置。（4）用axis(xmin xmax ymin ymax)控制坐标的刻度范围。（5）用grid on打开网格，grid off关闭网格。（6）用hold on把要画的图形保持在之前在同一figure上所画的图形中（同一坐标系）。（7）图5-3中的两“点线”直线，一条的两个端点为(0,1)和(1,0)，另一条的两个端点为(0,2)和(1.6,0)。四、实验结果及其分析8-1捕鱼业的持续收获 产量模型8-2 种群的相互竞争（1）1.实验代码：syms r1 r2 N1 N2 k1 k2 x1 x2;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,x2;fx1=diff(f,x1);fx2=diff(f,x2);gx1=diff(g,x1);gx2=diff(g,x2);A=fx1,fx2;gx1,gx2syms x1 x2;p=subs(-(fx1+gx2),x1,x2,P(:,1),P(:,2);p=simple(p);q=subs(det(A),x1,x2,P(:,1),P(:,2);q=simple(q);P p q2.实验结果8-3 种群的相互竞争（2）1.实验代码%脚本文件function f=fun(t,x) r1=2.5;r2=1.8; N1=1.6;N2=1; k1=0.5;k2=1.6; f=r1*x(1).*(1-x(1)./N1-k1.*(x(2)./N2);r2*x(2).*(1-k2*x(1)./N1-x(2)./N2); %t,x=ode45(fun,0 8,0.1;0.1); plot(t,x(:,1),t,x(:,2); axis(0 8 0 2); grid on; t,x=ode45(fun,0 8,1;2); plot(t,x(:,1),t,x(:,2); axis(0 8 0 2); grid on; t,x=ode45(fun,0 8,0.1;0.1); plot(x(:,1),x(:,2); hold on; t,x=ode45(fun,0 8,1;2); plot(x(:,1),x(:,2); x1=0:0.1:1; y1=-x1+1; x2=0:0.1:2; y2=-1.25*x2+2;2.实验结果