本节书摘来自异步社区出版社《MATLAB智能算法超级学习手册》一书中的第2章，第2.1节，作者：MATLAB技术联盟 , 高飞 , 许玢更多章节内容可以访问云栖社区“异步社区”公众号查看。
2.1 种群竞争微分方程模型
MATLAB智能算法超级学习手册
种群竞争模型是一个动态的过程，种群生存期间有着出生、死亡、迁入、迁出等问题。因此，种群数量较难确定，其种群竞争的数学模型只能通过反复的修正不断地完善，从而更加接近实际。本节就种群竞争微分方程模型的求解展开讨论。
并给定参数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);
%在此函数中，改变r1、r2、n1、n2、s1、s2的值，达到相关要求
针对题目中已知的初始条件，编写相应的MATLAB脚本文件程序，运行得图2-1、图2-2。
由图2-1、图2-2可知：在t=10时，x达到稳定值100，y达到稳定值0。
结论：时间t充分大时，x(t)、y(t)的值稳定在x=100,y=0。
图2-1的程序如下。
%绘制当r1=1，r2=1，n1=100，n2=100，s1=0.5，s2=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的程序如下。
%自定义种群函数
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);
运行的脚本文件程序如下。
% 保持s1、s2不变，改变其他变量
>> 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中，当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所示。
由这两个图可知：
① 在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
图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);