文章目录
前言
今天我们要说的就是利用MATLAB对常微分方程进行数值求解的实例应用-Lotka-Volterra模型及其改进模型。本文是科学计算与MATLAB语言专题六第6小节的学习笔记,如果大家有时间的话,可以去听听课,没有的话,可以看看我的笔记,还是很不错的。
1 L o t k a − V o l t e r r a Lotka-Volterra Lotka−Volterra模型
20世纪40年代,Lotka(1925)和Volterra(1926)奠定了种间竞争关系的理论基础,他们提出的种间竞争方程对现代生态学理论的发展有着重大影响。
可以用一阶非线性微分方程组表示:
{
d
r
d
t
=
2
r
−
λ
r
f
,
r
(
0
)
=
r
0
d
f
d
t
=
−
f
+
λ
r
f
,
f
(
0
)
=
f
0
\left\{ \begin{aligned} \frac{dr}{dt}&=2r-\lambda rf,r(0)=r_0\\ \frac{df}{dt}&=-f+\lambda rf,f(0)=f_0 \end{aligned} \right.
⎩⎪⎪⎨⎪⎪⎧dtdrdtdf=2r−λrf,r(0)=r0=−f+λrf,f(0)=f0
其中,
t
t
t为时间,
r
(
t
)
r(t)
r(t)为兔子数量,
f
(
t
)
f(t)
f(t)为狐狸数量,
λ
\lambda
λ为一个正常数。该系统解具有周期性,其周期取决于初始条件。也就是说,对于任意数量的
r
(
0
)
r(0)
r(0)和
f
(
0
)
f(0)
f(0),总存在一个时间
t
=
t
p
t=t_p
t=tp,使这两个种群的数量回归初始值。
①在
r
0
=
300
、
f
0
=
150
、
λ
=
0.01
r_0=300、f_0=150、\lambda =0.01
r0=300、f0=150、λ=0.01的条件下,求该系统的解,结果会发现
t
p
t_p
tp接近5。进一步绘制
r
(
t
)
r(t)
r(t)和
f
(
t
)
f(t)
f(t)函数的曲线,以及以r和f为坐标轴画相平面图。
相平面相关概念:
设一个二阶系统可以用下列常微分方程来描述:
d
2
x
d
t
2
+
a
1
(
x
,
d
x
d
t
)
d
x
d
t
+
a
0
(
x
,
d
x
d
t
)
x
=
0
,
x
¨
=
f
(
x
,
x
˙
)
\frac{d^2x} {dt^2}+a_1(x,\frac{dx}{dt})\frac{dx}{dt}+a_0(x,\frac{dx}{dt})x=0 ,\ddot{x}=f(x,\dot{x})
dt2d2x+a1(x,dtdx)dtdx+a0(x,dtdx)x=0,x¨=f(x,x˙)
令
x
=
x
1
,
d
x
/
d
t
=
x
2
x=x1,dx/dt=x2
x=x1,dx/dt=x2
d
x
2
d
x
1
=
f
(
x
1
,
x
2
)
x
2
\frac{dx_2}{dx_1}=\frac{f(x_1,x_2)}{x_2}
dx1dx2=x2f(x1,x2)
{
d
x
1
d
t
=
x
2
d
x
2
d
t
=
f
(
x
1
,
x
2
)
\left\{ \begin{aligned} \frac{dx_1}{dt}&=x_2\\ \frac{dx_2}{dt}&=f(x_1,x_2)\\ \end{aligned} \right.
⎩⎪⎪⎨⎪⎪⎧dtdx1dtdx2=x2=f(x1,x2)
以
x
1
x_1
x1为自变量,以
x
2
x_2
x2为因变量的一阶微分方程。二阶系统常微分方程方程的解既可用x与t的关系来表示,也可用
x
2
x_2
x2与
x
1
x_1
x1的关系来表示。实际上,看作一个质点的运动方程,则x1(t)代表质点的位置,
x
2
(
t
)
x_2(t)
x2(t)代表质点的速度。以
x
1
x_1
x1为自变量,以
x
2
x_2
x2为因变量的一阶微分方程。二阶系统常微分方程方程的解既可用x与t的关系来表示,也可用x,与x1的关系来表示。实际上,看作一个质点的运动方程,则
x
1
(
t
)
x_1(t)
x1(t)代表质点的位置,
x
2
(
t
)
x_2(t)
x2(t)代表质点的速度。
用
x
1
、
x
2
x_1、x_2
x1、x2,描述二阶系统常微分方程方程的解,也就是用质点的状态来表示该质点的运动。在物理学中,状态又称为相。把由
x
1
,
x
2
x_1,x_2
x1,x2,所组成的平面坐标系称为相平面,系统的一个状态则对应于相平面上的一个点。
当t变化时,系统状态在相平面上移动的轨迹称为相轨迹。
而与不同初始状态对应的一簇相轨迹所组成的图叫做相平面图。
②在
r
0
=
15
、
f
0
=
22
、
λ
=
0.01
r_0=15、f_0=22、\lambda=0.01
r0=15、f0=22、λ=0.01时,求解并绘制图形,发现会发现
t
p
t_p
tp接近6.62。
③在
r
0
=
102
、
f
0
=
198
、
λ
=
0.01
r_0=102、f_0=198、\lambda=0.01
r0=102、f0=198、λ=0.01时求解并绘制图形,并确定周期
t
p
t_p
tp。
④点
(
r
0
,
f
0
)
=
(
1
/
λ
,
2
/
λ
)
(r_0,f_0)=(1/\lambda,2/\lambda)
(r0,f0)=(1/λ,2/λ)是一个稳定平衡点。若以此为初值,那么种群数量不发生变化。若初值接近此平衡点,那么数量不会发生大的变化。试作图验证。
模型分析
x
1
(
t
)
x_1(t)
x1(t):兔子在t时刻的数量
x
2
(
t
)
x_2(t)
x2(t)一狐狸在t时刻的数量
r
1
r_1
r1一兔子独立生存时的增长率
r
2
r_2
r2一狐狸独自存在时的死亡率
λ
1
\lambda_1
λ1,一狐狸掠取兔子的能力
λ
2
\lambda_2
λ2一兔子对狐狸的供养能力
这里假设
r
1
=
2
,
r
2
=
1
,
λ
1
=
λ
2
=
0.01
r1=2,r2=1,\lambda_1=\lambda_2=0.01
r1=2,r2=1,λ1=λ2=0.01,即系统模型如下:
{
d
x
1
d
t
=
x
1
(
2
−
0.01
x
2
)
d
x
2
d
t
=
x
2
(
−
1
+
0.01
x
1
)
\left\{ \begin{aligned} \frac{dx_1}{dt}&=x_1(2-0.01x_2)\\ \frac{dx_2}{dt}&=x_2(-1+0.01x_1) \end{aligned} \right.
⎩⎪⎪⎨⎪⎪⎧dtdx1dtdx2=x1(2−0.01x2)=x2(−1+0.01x1)
第①问:兔子数量初始值 x 1 ( 0 ) = 300 x_1(0)=300 x1(0)=300,狐狸数量初始值 x 2 ( 0 ) = 150 x_2(0)=150 x2(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_1(0)=15
x1(0)=15,狐狸初始值
x
2
(
0
)
=
22
x_2(0)=22
x2(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_1(0)=102 x1(0)=102,狐狸初始值 x 2 ( 0 ) = 198 x_2(0)=198 x2(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
/
λ
,
2
/
λ
)
(1/\lambda,2/\lambda)
(1/λ,2/λ)是稳定平衡点。
取
λ
\lambda
λ=0.01,所以稳定平衡点
(
1
/
λ
,
2
/
λ
)
(1/\lambda,2/\lambda)
(1/λ,2/λ)即是
(
100
,
200
)
(100,200)
(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
)
(98,195)
(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
)
(70,150)
(70,150)时(向下偏离平衡点比较远时),画出其图像。
可以发现兔子和狐狸的数量变化比较剧烈,但是依然成周期性变化。
当将初始值变为
(
900
,
1600
)
(900,1600)
(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 L o t k a − V o l t e r r a Lotka-Volterra Lotka−Volterra改进模型
修改后的方程:
{
d
r
d
t
=
2
(
1
−
r
R
)
−
λ
r
f
,
r
(
0
)
=
r
0
d
f
d
t
=
−
f
+
λ
r
f
,
f
(
0
)
=
f
0
\left\{ \begin{aligned} \frac{dr}{dt}&= 2(1-\frac{r}{R})-\lambda rf,r(0)=r0 \\ \frac{df}{dt}&= -f+\lambda rf,f(0)=f0 \end{aligned} \right.
⎩⎪⎪⎨⎪⎪⎧dtdrdtdf=2(1−Rr)−λrf,r(0)=r0=−f+λrf,f(0)=f0
其中,t为时间;r(t)为兔子数量;f(t)为狐狸数量;
λ
\lambda
λ为一个正常数,R也是一个正常数。方程在
r
0
=
300
、
f
0
=
150
r0=300、f0=150
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('改进模型下,狐狸数量相对于兔子数量的关系曲线');
模型的进一步改进:
考虑狐狸的环境容纳量
考虑到兔子其他天敌和狐狸天敌的影响
考虑自然灾害、人为捕捉等因素
结语
欢迎大家点赞👍,收藏⭐,转发🚀,
如有问题、建议,请您在评论区留言💬哦。