实验8-1 捕鱼业的持续收获 ——产量模型
一、实验目的
熟悉MATLAB编程，学会稳定性模型的分析，分析实际问题。
二、实验要求
填空后的M文件。
三、实验内容
6.1 捕鱼业的持续收获——产量模型
%文件名：p178.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=N(1-x/N) , x1=0
x0=x(2);
x1=x(1);%符号变量x的结构类型成为<2×1sym>
%求F(x)的微分F'(x)
syms x; %定义符号变量x的结构类型为<1×1sym>
dF=diff(Fx,'x');
dF=simple(dF) %简化符号表达式
%得 F'(x)=___________
%求F'(x0)并简化
dFx0=subs(dF,x,x0); %将x=x0代入符号表达式dF
dFx0=simple(dFx0)
%得 F’(x0)=E-r
%求F’(x1)=r-E
dFx1=subs(dF,x,x1)
%得 F’(x1)=r-E
%若 E<r，有F'(x0)<0，F'(x1)>0，故x0点稳定，x1点不稳定（根据平衡点稳定性的准则）；
%若 E>r，则结果正好相反。
%在渔场鱼量稳定在x0的前提下（E<r），求E使持续产量h(x0)达到最大hm。
%通过分析（见教材p178图1），只需求x0*使f(x)达到最大，且hm=f(x0*)。
syms r x N
fx=r*x*(1-x/N);
a) df=diff(fx,'x');
x0=solve(df,x)
%得 x0*=N/2
hm=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 种群的相互竞争（1）
一、 实验目的
熟悉MATLAB编程，学会稳定性模型的分析，分析实际问题。
二、 实验要求
1．补充的程序段。
2．运行结果。
三、 实验内容
%6.3 种群的相互竞争
%文件名：p186.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
%得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,gx2]
syms 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列，即稳定条件。
四、 实验结果与分析
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,gx2]
syms 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 q]
2.实验结果
实验8-3 种群的相互竞争（2）
一、 实验目的
熟悉MATLAB编程，学会稳定性模型的分析，分析实际问题。
二、 实验要求
1．补充的程序段。
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)。
四、 实验结果与分析
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.实验结果
第二篇：数学模型试验报告册1
课程设计名称： 设计一：MATLAB软件入门 指导教师：
课程设计时数： 课程设计设备：安装了Matlab、C++软件的计算机 课程设计日期： 实验地点：
课程设计目的：
1. 熟悉MATLAB软件的用户环境；
2. 了解MATLAB软件的一般目的命令；
3. 掌握MATLAB数组操作与运算函数；
4. 掌握MATLAB软件的基本绘图命令；
4. 掌握MATLAB语言的几种循环、条件和开关选择结构。
课程设计准备：
1. 在开始本实验之前，请回顾相关内容；
2. 需要一台准备安装Windows XP Professional操作系统和装有数学软件的计算机。
课程设计内容及要求
要求：设计过程必须包括问题的简要叙述、问题分析、实验程序及注释、实验数据及结果分析和实验结论几个主要部分。
1. 采用向量构造符得到向量[1,4,7, 程序： ,31]。
x=1:3:31
2. 随机产生一向量x，求向量x的最大值。
程序：
x=rand(10,1) %产生一个10维向量
a=max(x)
3. 利用列向量(1,2,3,,6)T建立一个范德蒙矩阵A，并利用位于矩阵A的奇数行偶数列的元素建立一个新的矩阵B，须保持这些元素的相对位置不变。
程序：
T=[1 2 3 4 5 6]';
Vandermonde=zeros(length(T),length(T));
for i=1:1:length(T);
Vandermonde(i,:)=(T).^(i-1);
end
Vandermonde
Vandermonde(1:2:6,2:2:6)
4. 按水平和竖直方向分别合并下述两个矩阵：
?100?
A???110??234?
??,B???567??
?001????8910??
程序：
a=[1 0 0
b=[2 3 4
c=[a;b]
d=[a ,b]
结果：
c =
1 0 0
1 1 0
0 0 1
2 3 4
5 6 7
8 9 10
d =
1 0 0 2 3 4 1 1 0 5 6 7 0 0 1 8 9 10
5. 当n?100时，求y??n1
?12i?1的值。
i
程序：
n=100;
s=0;
for i=1:n
end
s 1 1 0 0 0 1]; 5 6 7 8 9 10]; s=s+1/(2*i-1);
结果：
s =
3.2843
6. 一个三位整数各位数字的立方和等于该数本身则称该数为水仙花数。输出全部水仙花数。
程序：
x=[];
for i=100:999
n1=fix(i/100);
n2=fix((i-n1*100)/10);%取出十位数
end
x n3=i-n1*100-n2*10;%取出个位数 if i==n1^3+n2^3+n3^3 x=[x i]; end
结果：
x =
153 370 371 407
7. 求[1000,2000]之间第一个被17整除的整数。
程序：
for i=1000:2000
if mod(i,17)==0
i
break;
end
end
结果：
i =
1003
8. 用MATLAB绘制两条曲线，x?[0,2?]，以?为步长，一条是正弦曲线，一条是余弦曲线，线宽10
为6个象素，正弦曲线为绿色，余弦曲线为红色，线型分别为实线和虚线，并给所绘的两条曲线增添图例，分别为“正弦曲线”和“余弦曲线”。
程序：
x=0:pi/10:2*pi;
plot(x,sin(x),'-g ')
hold on
plot(x,cos(x),':r','linewidth',6)
legend('正弦函数','余弦函数')
1
0.8
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-101234567
9. 在同一坐标内，分别用不同线型和颜色绘制曲线y1?0.2e?0.5xcos(4?x)和y2?2e?0.5xcos(?x),x?[0,2?]，并标记两曲线交叉点。
程序：
x=linspace(0,2*pi,1000);
for i=1:1:1000
y1(i)=0.2*cos(4*pi*x(i)).*(exp(-0.5*x(i))); y2(i)=2*cos(pi*x(i)).*(exp(-0.5*x(i)));
end
plot(x,y1,'r');
hold on
plot(x,y2,'g');
for i=1:1:1000
delt=abs(y1(i)-y2(i));
if delt<=1e-2
hold on
plot(x(i),y1(i),'k*');
end
end
2
1.5
1
0.5
-0.5
-1-1.
01234567
10. 分别在同一窗口的不同子图绘制y1?sint,y2?cost,y3?sin3t,y4?sint在区间[0,2?]上的图像。
程序：t=0:pi/100:2*pi;
subplot(2,2,1);
plot(t,sin(t))
subplot(2,2,2);
plot(t,cos(t))
subplot(2,2,3);
plot(t,sin(3*t))
subplot(2,2,4);
plot(t,abs(sin(t)))
课程设计总结：