如何利用MATLAB建立Lotka-Volterra模型及其改进模型

前言

今天我们要说的就是利用MATLAB对常微分方程进行数值求解的实例应用-Lotka-Volterra模型及其改进模型。本文是科学计算与MATLAB语言专题六第6小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。

1 $Lotka-Volterra$模型

20世纪40年代，Lotka（1925）和Volterra（1926）奠定了种间竞争关系的理论基础，他们提出的种间竞争方程对现代生态学理论的发展有着重大影响。
可以用一阶非线性微分方程组表示：
$$\{\begin{array}{rl}\frac{dr}{dt}& =2r-\lambda rf，r(0)=r_0\\ \frac{df}{dt}& =-f+\lambda rf，f(0)=f_0\end{array}$$

其中，$t$为时间，$r（t）$为兔子数量，$f（t）$为狐狸数量，$\lambda$为一个正常数。该系统解具有周期性，其周期取决于初始条件。也就是说，对于任意数量的$r(0)$和$f(0)$，总存在一个时间$t=t_p$，使这两个种群的数量回归初始值。
①在$r_0=300、f_0=150、\lambda =0.01$的条件下，求该系统的解，结果会发现$t_p$接近5。进一步绘制$r(t)$和$f(t)$函数的曲线，以及以r和f为坐标轴画相平面图。
相平面相关概念：
设一个二阶系统可以用下列常微分方程来描述：
$$\frac{d^{2}x}{d{t}^2}+a_1(x,\frac{dx}{dt})\frac{dx}{dt}+a_0(x,\frac{dx}{dt})x=0，\ddot{x}=f(x,\dot{x})$$
令$x=x1，dx/dt=x2$
$$\frac{d{x}_2}{d{x}_1}=\frac{f(x_1,x_2)}{x_2}$$
$$\{\begin{array}{rl}\frac{d{x}_1}{dt}& =x_2\\ \frac{d{x}_2}{dt}& =f(x_1,x_2)\end{array}$$
以$x_1$为自变量，以$x_2$为因变量的一阶微分方程。二阶系统常微分方程方程的解既可用x与t的关系来表示，也可用$x_2$与$x_1$的关系来表示。实际上，看作一个质点的运动方程，则x1(t)代表质点的位置，$x_2(t)$代表质点的速度。以$x_1$为自变量，以$x_2$为因变量的一阶微分方程。二阶系统常微分方程方程的解既可用x与t的关系来表示，也可用x，与x1的关系来表示。实际上，看作一个质点的运动方程，则$x_1(t)$代表质点的位置，$x_2(t)$代表质点的速度。
用$x_1、x_2$，描述二阶系统常微分方程方程的解，也就是用质点的状态来表示该质点的运动。在物理学中，状态又称为相。把由$x_1，x_2$，所组成的平面坐标系称为相平面，系统的一个状态则对应于相平面上的一个点。
当t变化时，系统状态在相平面上移动的轨迹称为相轨迹。
而与不同初始状态对应的一簇相轨迹所组成的图叫做相平面图。

②在$r_0=15、f_0=22、\lambda =0.01$时,求解并绘制图形，发现会发现$t_p$接近6.62。
③在$r_0=102、f_0=198、\lambda =0.01$时求解并绘制图形，并确定周期$t_p$。
④点$（r_0，f_0）=（1/\lambda ，2/\lambda ）$是一个稳定平衡点。若以此为初值，那么种群数量不发生变化。若初值接近此平衡点，那么数量不会发生大的变化。试作图验证。
模型分析
$x_1(t)$:兔子在t时刻的数量
$x_2(t）$一狐狸在t时刻的数量
$r_1$一兔子独立生存时的增长率
$r_2$一狐狸独自存在时的死亡率
${\lambda }_1$，一狐狸掠取兔子的能力
${\lambda }_2$一兔子对狐狸的供养能力
这里假设$r1=2，r2=1，{\lambda }_1={\lambda }_2=0.01$，即系统模型如下：
$$\{\begin{array}{rl}\frac{d{x}_1}{dt}& =x_1(2-0.01x_2)\\ \frac{d{x}_2}{dt}& =x_2(-1+0.01x_1)\end{array}$$

第①问：兔子数量初始值$x_1（0）=300$，狐狸数量初始值$x_2(0)=150$。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl（t）','x2（t）');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on

我们可以看到，兔子和狐狸两者相互制约，当狐狸逐渐增多时，兔子逐渐减少，当狐狸增加到一定数量时，由于种群内部竞争激烈，又导致其数量减少，从而兔子天敌减少，进而导致兔子数量增加，当兔子增加到一定数量时，种间竞争加剧，从而又导致兔子数量减少，如此循环，相互制约发展。
第②问：兔子初始值$x_1（0）=15$，狐狸初始值$x_2（0）=22$。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[15,22])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl（t）','x2（t）');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on

第③问：兔子初始值$x_1（0）=102$，狐狸初始值$x_2（0）=198$。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[102,198])
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('xl（t）','x2（t）');
xlabel('时间');
ylabel('物种数量');
grid on 
subplot(1,2,2);
plot(x(:,1),x(:,2))
grid on

第④问：验证$（1/\lambda ，2/\lambda ）$是稳定平衡点。
取$\lambda$=0.01，所以稳定平衡点$（1/\lambda ,2/\lambda ）$即是$(100,200)$，以此点作为初值时，画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[100,200])
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1（t）-兔子','x2（t）-狐狸');
xlabel('时间');
ylabel('物种数量');

当将初始值变为$(98,195)$时，即向下十分接近平衡点，画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[98,195])
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1（t）-兔子','x2（t）-狐狸');
xlabel('时间');
ylabel('物种数量');

两者图像没有交点，兔子和狐狸的数量分别以各自的平衡点进行周期性波动，但是两者的数量变化不是很大。
当将初始值变为$(70,150)$时（向下偏离平衡点比较远时），画出其图像。

可以发现兔子和狐狸的数量变化比较剧烈，但是依然成周期性变化。
当将初始值变为$(900,1600)$时（向上偏离平衡点十分远时），画出其图像。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,500],[900,1560])
plot(t,x(:,1),t,x(:,2));
legend('x1（t）-兔子','x2（t）-狐狸');
xlabel('时间');
ylabel('物种数量');

可以发现，此时兔子和狐狸的数量变化相当的剧烈，但是依然存在周期。

2 $Lotka-Volterra$改进模型

修改后的方程：
$$\{\begin{array}{rl}\frac{dr}{dt}& =2（1-\frac{r}{R})-\lambda rf，r(0)=r0\\ \frac{df}{dt}& =-f+\lambda rf，f(0)=f0\end{array}$$
其中，t为时间；r（t）为兔子数量；f（t）为狐狸数量；$\lambda$为一个正常数，R也是一个正常数。方程在$r0=300、f0=150$初始条件下的50个时间单位上求解并绘制图形。
①原模型下，狐狸数量和兔子数量的时间函数曲线
②改进模型下，狐狸数量和兔子数量的时间函数曲线。
③原模型下，狐狸数量相对兔子数量的关系曲线。
④改进模型下，狐狸数量相对兔子数量的关系曲线。
第①问：在原模型下，狐狸数量和兔子数量的时间函数曲线

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150])
plot(t,x(:,1),t,x(:,2));
legend('x1（t）-兔子','x2（t）-狐狸');
xlabel('时间');
ylabel('物种数量');
title('原模型下，狐狸和兔子数量的函数曲线');

第②问：在改进模型下，狐狸数量和兔子数量的时间函数曲线。

rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
    x(2)*(-1+0.01*x(1))];%...代表换行
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(t,x(:,1),t,x(:,2));
legend('x1（t）-兔子','x2（t）-狐狸');
xlabel('时间');
ylabel('物种数量');
title('改进模型下，狐狸和兔子数量的函数曲线');

第③问：在原模型下，狐狸数量相对兔子数量的关系曲线。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150])
plot(x(:,1),x(:,2));
xlabel('兔子数量');
ylabel('狐狸数量');
title('原模型下，狐狸数量相对于兔子数量的关系曲线');

第④问：在改进模型下，狐狸数量相对于兔子数量的关系曲线。

rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
    x(2)*(-1+0.01*x(1))];%...代表换行
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(x(:,1),x(:,2));
xlabel('时间');
ylabel('物种数量');
title('改进模型下，狐狸数量相对于兔子数量的关系曲线');

模型的进一步改进：
考虑狐狸的环境容纳量
考虑到兔子其他天敌和狐狸天敌的影响
考虑自然灾害、人为捕捉等因素

结语

欢迎大家点赞👍，收藏⭐，转发🚀，
如有问题、建议，请您在评论区留言💬哦。