种群相互竞争的Matlab程序

两种群相互竞争模型如下：
1112
2221(1)(1)dx x y r x s dt n n dy y x r y s dt n n =--=-- 其中x(t),y(t)分别为甲乙两种群的数量，1r ，2r 为它们的固有增长率，1n ，2
n 为它们的最大容量。

1s 的含义是，对于供养甲的资源来说，单位数量的乙（相对2n ）的消耗为单位数量甲（相对1n ）消耗的1s 倍，对2s 可以作相应解释。

经过计算，该模型无解析解，故用数值方法研究，为此提出以下问题：
（1） 设r1=r2=1,n1=n2=100,s1=0.5,s2=2,初值x0=y0=10,计算x(t),y(t),画出
它们的图形及图（x,y ），说明时间t 充分大了以后x(t),y(t)的变化趋
势。

（2） 改变r1,r2,n1,n2,x0,y0,但s1,s2不变（或保持s1<1,s2>1），计算并分
析所得结果，若s1=1.5(>1),s2=0.7(<1),再分析结果。

由此可以得到
什么结论，请作出解释。

（3） 试验当s1=0.8,s2=0.7时会有什么结果，当s1=1.5,s2=1.7时，又会有
什么结果。

模型求解：
程序如下：
fun.m:
function dx=fun(t,x,r1,r2,n1,n2,s1,s2)
dx=[r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)];
p3.m:
h=0.1;%所取时间点间隔
ts=[0:h:30];%时间区间
x0=[10,10];%初始条件
opt=odeset('reltol',1e-6,'abstol',1e-9);%相对误差1e-6，绝对误差1e-9
[t,x]=ode45(@fun,ts,x0,opt,1,1,100,100,0.5,2);%使用5级4阶龙格—库塔公式计算%后面的参数传给fun,分别是r1,r2,n1,n2,s1,s2
[t,x]%输出t,x(t),y(t)
plot(t,x,'.-'),grid%输出x1(t), x2(t)的图形
gtext('x1(t)'),gtext(' x2(t)'),pause
plot(x(:,1),x(:,2),'.-'),grid,%作相轨线
gtext('x1'),gtext('x2');
运行结果[t,x]为：
ans =
0 10.0000 10.0000
0.1000 10.8805 10.7120
0.2000 11.8235 11.4454
0.3000 12.8309 12.1962
0.4000 13.9044 12.9595
0.5000 15.0453 13.7295
……
29.4000 100.0000 0.0000
29.5000 100.0000 0.0000
29.6000 100.0000 0.0000
29.7000 100.0000 0.0000
29.8000 100.0000 0.0000
29.9000 100.0000 0.0000
30.0000 100.0000 0.0000
最后数值稳定在x=100,y=0上，即物种甲达到最大值，物种乙灭绝。

改变参数进一步讨论：
下面在保持s1,s2不变的基础上，分别改变r1,r2；n1,n2；x0,y0观察变化趋
势：
(1)改变r1,r2:
先不同的是变化速度减缓了，这是由于自然增长率r1,r2变小的缘故（相当于变化率减小）。

(2)改变n1,n2:
一度达到90以上，但最终仍然灭绝。

物种容量的改变并不能影响最终谁会灭绝。

下面的情况证明了这一点：
(3)改变x0,y0:
x0=10,y0=100:
量最大且乙物种灭绝的结果。

下面再改变s1,s2观察变化趋势：
(1)s1>1,s2<1
限。

如果这时改变r1,r2,n1,n2,x0,y0这些参数，变化趋势和上面列举的相同（甲乙相反），这从方程的对称性上可以求证。

现在得出结论，由s1,s2的物理意义，当某个s1或者s2大于1时（另一个小于1），它将严重消耗其作用的物种的生存资源，最终的结果是致使此物种灭绝。

(2)s1<1,s2<1
s1=0.8,s2=0.7
综上所述，s1,s2小于1时消耗生存资源的严重程度较轻，所以甲乙物种可以共存，但两者都达不到最大值；当其中之一大于1时，对应作用的物种就会由于生存资源的过度消耗而灭绝；当s1,s2都大于1时，两物种竞争激烈，最后s1,s2中更大者对应作用的物种灭绝。

所谓物尽天择，自然资源是有限的，需要更少资源就能生存的物种在竞争中将占有优势。

还有一种情况，当s1,s2都大于1但相等时，由于方程的对称性，甲乙两物种都能生存下来，但都不能达到最大值。