两种群间的相互竞争
摘要
本文针对两种群间的竞争问题作了详细的论述，主体分为两部分，第一部分主要通过理论分析的方法来阐述模型，第二部分主要利用MATLAB通过数值分析的方法从另一个角度来阐述模型，两个部分相辅相成，从不同的角度对同一个模型进行分析，并在最后得到一致的结果。另外本文在第一部分主要以理论的方式对模型进行数学上的描述，在第二部分主要以生物间的角度对模型进行描述，与此同时对第一部分作一个总结。
关键词： 稳定性 平面动力系统 增广相空间 轨线
一、问题提出
两种群竞争模型很好的描述了种群间的各种关系，而如果从发展的眼光来看待问题，我们不禁对两种群在未来很长一段时间内的状态产生兴趣，换句话说，我们要研究的是在无穷远的将来，两个种群的数量变化关系，这对我们进一步研究生物学的各种问题是有意义的。
二、基本假设
假设1： 有甲乙两个种群，它们独自生存时的数量变化服从Logistic规律。 假设2： 两种群一起生存时，乙种群对甲种群增长的阻滞作用与乙种群的数
量成正比，甲种群对乙种群增长的阻滞作用与甲种群的数量也成正比。
三、问题分析
根据“假设1”，我们容易得到方程组如下
x⎧dx(t)=rx(1-)1⎪dtn1⎪
（1） ⎨
⎪dy(t)=ry(1-y)
2
⎪n2⎩dt
其中x(t)，y(t)分别为甲乙两种群随时间变化的数量；r1，r2为它们的固有增长率；n1和n2为环境允许条件下，甲乙两种群的最大数量。
再由“假设2”，对方程组（1）变形，我们得到方程组如下
xy⎧dx(t)
=rx(1--s)11⎪dtnn⎪12
（2） ⎨
dy(t)xy⎪=r2y(1-s2-)⎪dtn1n2⎩
其中s1的含义是，对于供养甲种群的资源而言，单位数量乙（相对于n2）的消耗为单位数量甲（相对于n1）消耗的s1倍；s2的含义是，对于供养乙种群的资源而言，单位数量甲（相对于n1）的消耗为单位数量乙（相对于n2）消耗的s2倍。
我们所要研究的问题是当t→+∞时，x(t)与y(t)的极限状态，即稳定性，这是本文
的第一步，即通过理论方法予以研究。由于方程组（2）不存在解析解，因此我们只能求出它的数值解，并利用MATLAB中先画出增广相空间中的积分曲线，然后画出相空间中的轨线（即积分曲线沿t轴在相空间中的投影），进一步支持上述的理论研究，这是本文的第二步。
四、模型的建立与求解
第一部分：理论分析
现在我们回过头来看方程组（2）
xy⎧dx(t)
=rx(1--s)11⎪dtnn⎪12
⎨
dy(t)xy⎪=r2y(1-s2-)⎪n1n2⎩dt
可以看到，常微分方程组不显含自变量，因此方程组是自治的，现在我们令
xy⎧f(x,y)=rx(1--s)11⎪nn⎪12
（3） ⎨
xy⎪g(x,y)=ry(1-s-)22
⎪n1n2⎩
易知，f(x,y)和g(x,y)在xoy平面上连续，并且
r1s1x⎧∂f
=-⎪∂yn2⎪
⎨
sx∂g2y⎪=r(1-2-)2
⎪∂yn1n2⎩
也在xoy平面内连续，因此f(x,y)与g(x,y)满足Lipschitz条件，这就保证了柯西问题解的存在唯一性。因此易见，方程组（2）描绘的是一个平面上的动力系统。 方程组（2）作为一个平面上的动力系统，它不具有解析解，造成这个现象的主要原因是方程组（2）的右端非常系数且非线性。为了使分析得以进行下去，我们有必要对方程组（2）的右端做一些改变，其中一个想法就是将非线性的函数近似线性化，因为平面线性动力系统的理论已经比较完善，并且较易于判断稳定性。接下来我们将这一过程用数学语言描述出来。
首先，对方程组（2）的右端在奇点(x0,y0)处进行二元函数的泰勒展开（取一阶），将
右端近似线性化，如下
∂f⎧f(x,y)=f(x,y)+00⎪∂x⎪⎨
⎪g(x,y)=g(x,y)+∂g
00
⎪∂x⎩
∂f
∂y∂g
(x-x)+(x0,y0)0
∂y
(x0,y0)(x-x0)+
(x0,y0)
(y-y0)+o(ρ2)
(x0,y0)
(y-y0)+o(ρ2)
∂f⎧
f(x,y)≈⎪∂x⎪⇒⎨
⎪g(x,y)≈∂g⎪∂x⎩
其中
∂f
∂y∂g⋅x+(x0,y0)
∂y
(x0,y0)⋅x+
(x0,y0)
⋅y+C
（4）
(x0,y0)
⋅y+C'
ρ=∂f
∂x∂g
C'=g(x0,y0)-
∂xC=f(x0,y0)-
(x0,y0)⋅x0-
∂f
(x,y)⋅y0 ∂y00∂g⋅x-(x0,y0)0(x,y)⋅y0
∂y00
方程组（4）是去掉高阶项后近似线性化得到的。现在可以将方程组（2）写成如下向量矩阵的形式
dv
=Av+c （5）
dt
T
其中v=(x,y)
c=(C,C')T
∂f∂y∂g∂y⎫(x0,y0)⎪
⎪ ⎪(x0,y0)⎪
⎭
⎛∂f ∂xA=
∂g ∂x⎝
(x0,y0)
(x0,y0)
经过简单的计算，我们将（5）写成如下形式
2x0s1y0r1s1x0⎛⎫r(1--)- 1⎪xn1n2n2⎛⎫⎛C⎫d⎛x⎫ ⎪ ⎪= ⎪+ ⎪ （6）
r2s2y0s2x02y0⎪⎝y⎭⎝C'⎭dt⎝y⎭
-r2(1--)⎪ nnn⎝112⎭
现在令
p=-tr(A)=(2r1+s2r2)q=det(A)
（7） x0y
+(s1r1+2r2)0-(r1+r2)n1n2
=r1r2(1-
2x0s1y0sx2yrrssxy （8）-)(1-20-0)-121200n1n2n1n2n1n2
根据常微分方程的理论知，当p>0且q>0时，奇点(x0,y0)是稳定的，当p
点(x0,y0)是不稳定的。现在我们来求出方程组（2）的所有奇点，具体如下
xy⎧rx(1--s)=01⎪1
n1n2⎪
令 ⎨
⎪ry(1-sx-y)=022⎪n1n2⎩
得到四个奇点以并利用（7）和（8）算得相应的p值和q值如下
11
.奇点(x0,y0)=(0,0)
p1=-(r1+r2)
q1=r1r2
22.奇点(x0,y0)=(n1,0)
p2=r1-r2(1-s2)
q2=-r1r2(1-s2)
33.奇点(x0,y0)=(0,n2)
p3=-r1(1-s1)+r2
q3=-r1r2(1-s1)
44. 奇点(x0,y0)=(
n1(1-s1)n2(1-s2)
,)
1-s1s21-s1s2
p4=
r1(1-s1)+r2(1-s2)
1-s1s2
q4=
r1r2(1-s1)(1-s2)
1-s1s2
首先由模型可知，参数r1、r2、s1和s2都是大于0的。
现在我们开始来逐个分析四个奇点，先来讨论.根据常微分方程稳定性理论，
11
p1=-(r1+r2)
时灭绝的。
其次来讨论.
2⎧⎪p>0令 ⎨2
⎪⎩q>0
r1⎧
s>1-⎪2
r2⇒s2>1 ⇔⎨
⎪s>1⎩2
22
这说明在s2>1的条件下，奇点(x0,y0)=(n1,0)是稳定的。换句话说，在s2>1的条件下，
甲种群达到环境所允许的最大数量，乙种群最终会灭亡。 接下来讨论.
3
⎧⎪p>0令 ⎨3
⎪⎩q>0
r2⎧s>1-⎪1
r1⇒s1>1 ⇔⎨
⎪s>1⎩1
33
这说明在s1>1的条件下，奇点(x0,y0)=(0,n2)是稳定的。换句话说，在s1>1的条件下，
乙种群达到环境所允许的最大数量，甲种群最终会灭亡 最后来讨论.
4
⎧⎪p>0令 ⎨4
⎪⎩q>0
⎧s,s>1(不满足q>0,因此舍去)⇒⎨12⇒0
44
这说明在0
n1(1-s1)n2(1-s2)
,)是稳定的。换句话
1-s1s21-s1s2
说，在0
第二部分：数值解法
根据第一部分的理论分析，我们在这里根据第一部分的结果，分别给出不同s1、s2条
件下的数值解。在分情况之前，需对模型中的其他参数给定一个初值，这样更有利于我们接下来的分析，现作如下规定：r1=r2=1，n1=n2=1000，初值x(0)=y(0)=50，接下来我们开始分情况讨论。 情况（一）：s1=0.5,s2=2
利用MATLAB，编写M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=0.5;s2=2;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
然后在命令窗口输入如下命令（这里用5级4阶Runge-Kutta-Fehberg公式，并描绘从第0年到第50年的变化趋势） >> ts=0:0.1:50; >> x0=[50,50];
>> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1200
1000
800
600
400
200
-200
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)') 得到相空间中的轨线如下
200
150
100
50
y
0-500
200400
600x
[1**********]
可以看到，s1=0.5,s2=2对应于理论分析中的. 这很好的验证了之前的理论分析，
并且可以看出，在第10年左右的时候甲种群趋于环境所允许的最大数量，而乙种群趋于灭亡。
情况（二）：s1=1.9,s2=0.4
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.9;s2=0.4;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1200
1000
800
600
400
200
-200
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
1200
1000
800
600
y
400
200
0-50
050
100x
150200250
可以看到，s1=1.9,s2=0.4对应于理论分析中的. 这很好的验证了之前的理论分析，并且可以看出，在第10年左右的时候乙种群趋于环境所允许的最大数量，而甲种群趋于灭亡。
情况（三）：s1=0.8,s2=0.7
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=0.8;s2=0.7;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
700
600
500
400
300
200
100
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
700
600
500
400
y
[1**********]0
[1**********]0
300x
[**************]
s1=0.8,s2=0.7对应于理论分析中的. 这很好的验证了之前的理论分析，可以看到，
并且可以看出，两个种群一直保持相互依存的关系，由于s1>s2，所以乙种群的数量一直大于甲种群的数量。
情况（四）：s1=1.28,s2=1.3
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.28;s2=1.3;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1000900
[***********]20010000
5
10
15
20
25
30
35
40
45
50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
450400
350300250
y
[1**********]0
[1**********]00
500x
[**************]0
可以看到，在s1=1.28,s2=1.3条件下，甲乙两种群先相互竞争，但是由于s1
因此甲的竞争力大于乙，从图中可以看到，在第45年左右时，甲种群最终会获胜并趋于环境所允许的最大容量，而乙种群会趋于灭亡。同时，另一方面，这种结果是在理论分析中无法得到的，这也是我们进行数值分析方法的另外一个原因。
情况（五）：s1=1.33,s2=1.25
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.33;s2=1.25;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1000900
[***********]20010000
5
10
15
20
25
30
35
40
45
50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
1000900
[***********]20010000
50
100
150
200x
250
300
350
400
y
可以看到，在s1=1.33,s2=1.25的条件下，甲乙两种群先相互竞争，但是由于s1>s2，
因此乙的竞争力大于甲，从图中可以看到，在第35年左右时，乙种群最终会获胜并趋于环境所允许的最大容量，而甲种群会趋于灭亡。同时，另一方面，这种结果是在理论分析中无法得到的，这也是我们进行数值分析方法的另外一个原因。
五、参考文献
[1]姜启源等 《大学数学实验》（第2版） 清华大学出版社 2010年12月 [2]丁同仁等 《常微分方程教程》（第2版）高等教育出版社 2004年07月
两种群间的相互竞争
摘要
本文针对两种群间的竞争问题作了详细的论述，主体分为两部分，第一部分主要通过理论分析的方法来阐述模型，第二部分主要利用MATLAB通过数值分析的方法从另一个角度来阐述模型，两个部分相辅相成，从不同的角度对同一个模型进行分析，并在最后得到一致的结果。另外本文在第一部分主要以理论的方式对模型进行数学上的描述，在第二部分主要以生物间的角度对模型进行描述，与此同时对第一部分作一个总结。
关键词： 稳定性 平面动力系统 增广相空间 轨线
一、问题提出
两种群竞争模型很好的描述了种群间的各种关系，而如果从发展的眼光来看待问题，我们不禁对两种群在未来很长一段时间内的状态产生兴趣，换句话说，我们要研究的是在无穷远的将来，两个种群的数量变化关系，这对我们进一步研究生物学的各种问题是有意义的。
二、基本假设
假设1： 有甲乙两个种群，它们独自生存时的数量变化服从Logistic规律。 假设2： 两种群一起生存时，乙种群对甲种群增长的阻滞作用与乙种群的数
量成正比，甲种群对乙种群增长的阻滞作用与甲种群的数量也成正比。
三、问题分析
根据“假设1”，我们容易得到方程组如下
x⎧dx(t)=rx(1-)1⎪dtn1⎪
（1） ⎨
⎪dy(t)=ry(1-y)
2
⎪n2⎩dt
其中x(t)，y(t)分别为甲乙两种群随时间变化的数量；r1，r2为它们的固有增长率；n1和n2为环境允许条件下，甲乙两种群的最大数量。
再由“假设2”，对方程组（1）变形，我们得到方程组如下
xy⎧dx(t)
=rx(1--s)11⎪dtnn⎪12
（2） ⎨
dy(t)xy⎪=r2y(1-s2-)⎪dtn1n2⎩
其中s1的含义是，对于供养甲种群的资源而言，单位数量乙（相对于n2）的消耗为单位数量甲（相对于n1）消耗的s1倍；s2的含义是，对于供养乙种群的资源而言，单位数量甲（相对于n1）的消耗为单位数量乙（相对于n2）消耗的s2倍。
我们所要研究的问题是当t→+∞时，x(t)与y(t)的极限状态，即稳定性，这是本文
的第一步，即通过理论方法予以研究。由于方程组（2）不存在解析解，因此我们只能求出它的数值解，并利用MATLAB中先画出增广相空间中的积分曲线，然后画出相空间中的轨线（即积分曲线沿t轴在相空间中的投影），进一步支持上述的理论研究，这是本文的第二步。
四、模型的建立与求解
第一部分：理论分析
现在我们回过头来看方程组（2）
xy⎧dx(t)
=rx(1--s)11⎪dtnn⎪12
⎨
dy(t)xy⎪=r2y(1-s2-)⎪n1n2⎩dt
可以看到，常微分方程组不显含自变量，因此方程组是自治的，现在我们令
xy⎧f(x,y)=rx(1--s)11⎪nn⎪12
（3） ⎨
xy⎪g(x,y)=ry(1-s-)22
⎪n1n2⎩
易知，f(x,y)和g(x,y)在xoy平面上连续，并且
r1s1x⎧∂f
=-⎪∂yn2⎪
⎨
sx∂g2y⎪=r(1-2-)2
⎪∂yn1n2⎩
也在xoy平面内连续，因此f(x,y)与g(x,y)满足Lipschitz条件，这就保证了柯西问题解的存在唯一性。因此易见，方程组（2）描绘的是一个平面上的动力系统。 方程组（2）作为一个平面上的动力系统，它不具有解析解，造成这个现象的主要原因是方程组（2）的右端非常系数且非线性。为了使分析得以进行下去，我们有必要对方程组（2）的右端做一些改变，其中一个想法就是将非线性的函数近似线性化，因为平面线性动力系统的理论已经比较完善，并且较易于判断稳定性。接下来我们将这一过程用数学语言描述出来。
首先，对方程组（2）的右端在奇点(x0,y0)处进行二元函数的泰勒展开（取一阶），将
右端近似线性化，如下
∂f⎧f(x,y)=f(x,y)+00⎪∂x⎪⎨
⎪g(x,y)=g(x,y)+∂g
00
⎪∂x⎩
∂f
∂y∂g
(x-x)+(x0,y0)0
∂y
(x0,y0)(x-x0)+
(x0,y0)
(y-y0)+o(ρ2)
(x0,y0)
(y-y0)+o(ρ2)
∂f⎧
f(x,y)≈⎪∂x⎪⇒⎨
⎪g(x,y)≈∂g⎪∂x⎩
其中
∂f
∂y∂g⋅x+(x0,y0)
∂y
(x0,y0)⋅x+
(x0,y0)
⋅y+C
（4）
(x0,y0)
⋅y+C'
ρ=∂f
∂x∂g
C'=g(x0,y0)-
∂xC=f(x0,y0)-
(x0,y0)⋅x0-
∂f
(x,y)⋅y0 ∂y00∂g⋅x-(x0,y0)0(x,y)⋅y0
∂y00
方程组（4）是去掉高阶项后近似线性化得到的。现在可以将方程组（2）写成如下向量矩阵的形式
dv
=Av+c （5）
dt
T
其中v=(x,y)
c=(C,C')T
∂f∂y∂g∂y⎫(x0,y0)⎪
⎪ ⎪(x0,y0)⎪
⎭
⎛∂f ∂xA=
∂g ∂x⎝
(x0,y0)
(x0,y0)
经过简单的计算，我们将（5）写成如下形式
2x0s1y0r1s1x0⎛⎫r(1--)- 1⎪xn1n2n2⎛⎫⎛C⎫d⎛x⎫ ⎪ ⎪= ⎪+ ⎪ （6）
r2s2y0s2x02y0⎪⎝y⎭⎝C'⎭dt⎝y⎭
-r2(1--)⎪ nnn⎝112⎭
现在令
p=-tr(A)=(2r1+s2r2)q=det(A)
（7） x0y
+(s1r1+2r2)0-(r1+r2)n1n2
=r1r2(1-
2x0s1y0sx2yrrssxy （8）-)(1-20-0)-121200n1n2n1n2n1n2
根据常微分方程的理论知，当p>0且q>0时，奇点(x0,y0)是稳定的，当p
点(x0,y0)是不稳定的。现在我们来求出方程组（2）的所有奇点，具体如下
xy⎧rx(1--s)=01⎪1
n1n2⎪
令 ⎨
⎪ry(1-sx-y)=022⎪n1n2⎩
得到四个奇点以并利用（7）和（8）算得相应的p值和q值如下
11
.奇点(x0,y0)=(0,0)
p1=-(r1+r2)
q1=r1r2
22.奇点(x0,y0)=(n1,0)
p2=r1-r2(1-s2)
q2=-r1r2(1-s2)
33.奇点(x0,y0)=(0,n2)
p3=-r1(1-s1)+r2
q3=-r1r2(1-s1)
44. 奇点(x0,y0)=(
n1(1-s1)n2(1-s2)
,)
1-s1s21-s1s2
p4=
r1(1-s1)+r2(1-s2)
1-s1s2
q4=
r1r2(1-s1)(1-s2)
1-s1s2
首先由模型可知，参数r1、r2、s1和s2都是大于0的。
现在我们开始来逐个分析四个奇点，先来讨论.根据常微分方程稳定性理论，
11
p1=-(r1+r2)
时灭绝的。
其次来讨论.
2⎧⎪p>0令 ⎨2
⎪⎩q>0
r1⎧
s>1-⎪2
r2⇒s2>1 ⇔⎨
⎪s>1⎩2
22
这说明在s2>1的条件下，奇点(x0,y0)=(n1,0)是稳定的。换句话说，在s2>1的条件下，
甲种群达到环境所允许的最大数量，乙种群最终会灭亡。 接下来讨论.
3
⎧⎪p>0令 ⎨3
⎪⎩q>0
r2⎧s>1-⎪1
r1⇒s1>1 ⇔⎨
⎪s>1⎩1
33
这说明在s1>1的条件下，奇点(x0,y0)=(0,n2)是稳定的。换句话说，在s1>1的条件下，
乙种群达到环境所允许的最大数量，甲种群最终会灭亡 最后来讨论.
4
⎧⎪p>0令 ⎨4
⎪⎩q>0
⎧s,s>1(不满足q>0,因此舍去)⇒⎨12⇒0
44
这说明在0
n1(1-s1)n2(1-s2)
,)是稳定的。换句话
1-s1s21-s1s2
说，在0
第二部分：数值解法
根据第一部分的理论分析，我们在这里根据第一部分的结果，分别给出不同s1、s2条
件下的数值解。在分情况之前，需对模型中的其他参数给定一个初值，这样更有利于我们接下来的分析，现作如下规定：r1=r2=1，n1=n2=1000，初值x(0)=y(0)=50，接下来我们开始分情况讨论。 情况（一）：s1=0.5,s2=2
利用MATLAB，编写M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=0.5;s2=2;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
然后在命令窗口输入如下命令（这里用5级4阶Runge-Kutta-Fehberg公式，并描绘从第0年到第50年的变化趋势） >> ts=0:0.1:50; >> x0=[50,50];
>> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1200
1000
800
600
400
200
-200
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)') 得到相空间中的轨线如下
200
150
100
50
y
0-500
200400
600x
[1**********]
可以看到，s1=0.5,s2=2对应于理论分析中的. 这很好的验证了之前的理论分析，
并且可以看出，在第10年左右的时候甲种群趋于环境所允许的最大数量，而乙种群趋于灭亡。
情况（二）：s1=1.9,s2=0.4
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.9;s2=0.4;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1200
1000
800
600
400
200
-200
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
1200
1000
800
600
y
400
200
0-50
050
100x
150200250
可以看到，s1=1.9,s2=0.4对应于理论分析中的. 这很好的验证了之前的理论分析，并且可以看出，在第10年左右的时候乙种群趋于环境所允许的最大数量，而甲种群趋于灭亡。
情况（三）：s1=0.8,s2=0.7
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=0.8;s2=0.7;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
700
600
500
400
300
200
100
[***********]50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
700
600
500
400
y
[1**********]0
[1**********]0
300x
[**************]
s1=0.8,s2=0.7对应于理论分析中的. 这很好的验证了之前的理论分析，可以看到，
并且可以看出，两个种群一直保持相互依存的关系，由于s1>s2，所以乙种群的数量一直大于甲种群的数量。
情况（四）：s1=1.28,s2=1.3
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.28;s2=1.3;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1000900
[***********]20010000
5
10
15
20
25
30
35
40
45
50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
450400
350300250
y
[1**********]0
[1**********]00
500x
[**************]0
可以看到，在s1=1.28,s2=1.3条件下，甲乙两种群先相互竞争，但是由于s1
因此甲的竞争力大于乙，从图中可以看到，在第45年左右时，甲种群最终会获胜并趋于环境所允许的最大容量，而乙种群会趋于灭亡。同时，另一方面，这种结果是在理论分析中无法得到的，这也是我们进行数值分析方法的另外一个原因。
情况（五）：s1=1.33,s2=1.25
修改M文件如下
function aaa=bbb(t,x)
r1=1;r2=1;n1=1000;n2=1000;s1=1.33;s2=1.25;
aaa=diag([r1*(1-x(1)/n1-s1*x(2)/n2),r2*(1-x(2)/n2-s2*x(1)/n1)])*x;
保存后直接执行如下命令 >> [t,x]=ode45(@bbb,ts,x0); >> plot(t,x),grid,
>> gtext('\fontsize{20}x(t)'),gtext('\fontsize{20}y(t)')
得到增广相空间中的积分曲线如下（蓝色是x(t)，绿色是y(t)）
1000900
[***********]20010000
5
10
15
20
25
30
35
40
45
50
继续在命令窗口中输入如下命令 >> pause,plot(x(:,1),x(:,2)),grid, >> xlabel('x'),ylabel('y')
>> gtext('\fontsize{20}y=y(x)')
得到相空间中的轨线如下
1000900
[***********]20010000
50
100
150
200x
250
300
350
400
y
可以看到，在s1=1.33,s2=1.25的条件下，甲乙两种群先相互竞争，但是由于s1>s2，
因此乙的竞争力大于甲，从图中可以看到，在第35年左右时，乙种群最终会获胜并趋于环境所允许的最大容量，而甲种群会趋于灭亡。同时，另一方面，这种结果是在理论分析中无法得到的，这也是我们进行数值分析方法的另外一个原因。
五、参考文献
[1]姜启源等 《大学数学实验》（第2版） 清华大学出版社 2010年12月 [2]丁同仁等 《常微分方程教程》（第2版）高等教育出版社 2004年07月