下载掌阅APP，畅读海量书库
立即打开
种群竞争模型是一个动态的过程，种群生存期间有着出生、死亡、迁入、迁出等问题。因此，种群数量较难确定，其种群竞争的数学模型只能通过反复的修正不断地完善，从而更加接近实际。本节就种群竞争微分方程模型的求解展开讨论。
设有甲、乙两种群，当它们独自生存时，数量演变服从Logistic规律：
其中，x(t)、y(t)分别为甲、乙两种群的数量，r 1 、r 2 为它们的固有增长率，n 1 、n 2 为它们的最大容量。
当两种群在同一环境中生存时，它们之间的一种关系是为了争夺同一资源而进行竞争。考察乙消耗有限的资源对甲的增长产生的影响，可以合理地将种群甲的方程修改为：
式（2.2）中s 1 的含义：对于供养甲的资源而言，单位数量乙（相对于n 2 ）的消耗为单位数量甲（相对n 1 ）消耗的s 1 倍。类似地，如果甲的存在也影响了乙的增长，则乙的方程应改为：
式（2.2）中s 2 的含义：对于供养乙的资源而言，单位数量甲（相对于n 1 ）的消耗为单位数量乙（相对n 2 ）消耗的s 1 倍。给定种群的初始值：
并给定参数r 1 、r 2 、s 1 、s 2 、n 1 、n 2 后，由式（2.1）～式（2.4）可确定两种群数量的变化规律。
（1）设r 1 =r 2 =1，n 1 =n 2 =100，s 1 =0.5，s 2 =2，x 0 =y 0 =10，计算x(t)、y(t)，画出它们的图形及相图x(t)、y(t)，说明时间t充分大时x(t)、y(t)的变化趋势。
解：对于微分方程的求解，首先建立微分方程函数多数情况下，用数值解代替代数解进行方程的模拟，自定义种群函数程序zhongqun( )：
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1;r2=1;
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
%在此函数中，改变r 1 、r 2 、n 1 、n 2 、s 1 、s 2 的值，达到相关要求
针对题目中已知的初始条件，编写相应的MATLAB脚本文件程序，运行得图2-1、图2-2。
图2-1 x(t)、 y(t) 的变化趋势图
图2-2 解曲线
由图2-1、图2-2可知：在t=10时，x达到稳定值100，y达到稳定值0。
结论：时间t充分大时，x(t)、y(t)的值稳定在x=100,y=0。
图2-1的程序如下。
%绘制当r 1 =1，r 2 =1，n 1 =100，n 2 =100，s 1 =0.5，s 2 =2时的函数图像
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=0.5;s2=2;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-2的程序如下。
%绘制曲线向量解曲线
>>syms r1 r2 s1 s2 n1
r1=1;r2=1;s1=0.5;s2=2;n1=100;n2=100;
Xmin=0;
Xmax=140;
Ymin=0;
Ymax=100;
n=50;
% 计算切线矢量
>> [X，Y]=meshgrid(linspace(Xmin，Xmax，n)，linspace(Ymin，Ymax，n));
>> Fx=r1.*X.*(1-X./n1-s1.*Y./n2);
Fy=r2.*Y.*(1- s2.*X./n1-Y./n2);
Fx=Fx./(sqrt(Fx.^2+Fy.^2+1));
Fy=Fy./(sqrt(Fx.^2+Fy.^2+1));
%求解微分方程
>> options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
>> [T1，Y1]=ode45(@zhongqun，[0 50]，[10 10]，options);
>>% 绘制斜率场
hold on
grid on
box on
axis([Xmin，Xmax，Ymin，Ymax])
quiver(X，Y，Fx，Fy，0.5);
>>% 绘制解曲线
plot(Y1(:，1)，Y1(:，2)，'g'，'LineWidth'，2)
（2）改变r 1 、r 2 、n 1 、n 2 、x 0 、y 0 ，维持s 1 、s 2 不变，绘出r 1 =1.2，r 2 =1.1，n 1 =200，n 2 =120，x 0 =y 0 =10时的函数图像及r 1 =0.9,r 2 =1.5,n 1 =500,n 2 =800，x 0 =y 0 =10时的函数图像。
解：同第（1）问，改变初始值，程序运行后依次如图2-3、图2-4所示。图2-3为r 1 =1.2，r 2 =1 . 1,n 1 =200，n 2 =120，x 0 =y 0 =10时的函数图像；图2-4 为r 1 =0.9，r 2 =1.5，n 1 =500，n 2 =800，x 0 =y 0 =10时的函数图像。
图2-3 函数图像（1）
图2-4 函数图像（2）
图2-3的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.2;r2=1.1;
n1=200;n2=120;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
% 保持s 1 、s 2 不变，改变其他变量
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1.2;r2=1.1;s1=0.5;s2=2;n1=200;n2=120;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-4的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=0.9=1.5
n1=500;n2=800;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=0.9;r2=1.5;s1=0.5;s2=2;n1=500;n2=800;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
改变r 1 、r 2 、n 1 、n 2 、x 0 、y 0 ，维持s 1 、s 2 不变，种群甲将占优势，而种群乙变为0。
综合结论：改变r、n的初始值，甲、乙种群的最终稳定状态不会改变，都是种群甲达到环境最大承载值，而种群乙变为0。参数r、n的初始值的改变仅会影响达到稳定的速度，不会改变优势种群甲的优势地位，即最终的稳定状态情况。
若s 1 =1.5，s 2 =0.7，绘出函数图像，如图2-5所示。
图2-5 s 1 =1.5，s 2 =0 . 7
在图2-5中，当t=20时，y达到最大容量稳定值100，种群x变为零。当s 1 小而s 2 大时（s 1 与s 2 相差较大时，s 1 <1，s 2 >1），乙消耗甲的资源少，乙对甲影响小，同时甲消耗乙的资源多，甲对乙影响大，从而乙处于不利地位，甲处于有利地位，所以最后种群甲达到环境最大承载量，而种群乙则变为0。
相反，当s 1 大而s 2 小时（s 1 与s 2 相差较大时，s 1 >1，s 2 <1），乙消耗甲的资源多，乙对甲影响大，同时甲消耗乙的资源少，甲对乙影响小，乙处于有利地位，甲处于不利地位，所以最后种群乙达到环境最大承载量，而种群甲则变为0。
程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.5;r2=0.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=1.5;s2=0.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
（3）试验当s 1 =0.8，s 2 =0.7时会有什么结果；当s 1 =1.5，s 2 =1.7时又会有什么结果。你能解释这些结果吗？
解：当s 1 =0.8，s 2 =0.7时，绘制图像如图2-6所示；当s 1 =1.5，s 2 =1.7时，绘制图像，如图2-7所示。
图2-6 s 1 =0.8，s 2 =0 . 7
图2-7 s 1 =1.5，s 2 =1 . 7
由这两个图可知：
①在s 1 <1，s 2 <1，且s 1 、s 2 相近时，由于双方的抑制作用数值较小，所以任何一方都无法占据绝对优势，所以双方均无法消灭对方，只能在最大承载量以下达到稳定状态，且受到影响小的（图2-6、图2-7中乙受到的影响小）稳定值大。
②当s 1 >1，s 2 >1，且s 1 、s 2 相近时，由于s 1 、s 2 均比较大，对于种群的抑制力非常大，以至于一旦竞争处于下风，就会受到很大的抑制作用而无法生存。s 1 、s 2 相近，仍然有一方会灭绝。但这是一个较长的过程，一方无法短时间内完全抑制对方，另一方因受到抑制，故增长也会减缓。虽然r、n初值均相等，但由于s 1 <s 2 ，所以最后仍是甲受到的影响小，最终胜出。第（1）问中t=10时达到稳定，而此时t=20时才达到稳定，足足晚了10秒。
图2-6的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=0.8;r2=0.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=0.8;s2=0.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-7的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.5;r2=1.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=1.5;s2=1.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);