双时滞四维捕食网络的分析【基于matlab的动力学模型学习笔记_5】
作者:互联网
/*仅当作学习笔记,若有纰漏欢迎友好交流指正,此外若能提供一点帮助将会十分荣幸*/
本系列谈论过单时滞,但还没提及过双时滞,本文将着重介绍一种双时滞系统并对其进行简单处理分析。
摘 要:本文针对一个捕食网络,根据其特性建立起一个双时滞四维系统。首先分析其稳定性,找到其正平衡点,然后通过调试时滞量τ来模拟预测网络模型的发展前景。并根据所求得临界滞量τ0,带入实值仿真出其Hopf分岔图形。
0引言
自然界中,广泛存在着种群之间的捕食与被捕食关系,他们一起构成了生物链,为流传的“大鱼吃小鱼,小鱼吃虾米”就通俗易懂的道出了生物链之间的关系。而为了保护生态链的完整,甚至保证整个生态系统的有序发展,对其中捕食系统的研究就存在相当的意义。
在捕食系统中,最经典的就是Lotka-Volterra捕食模型(简称L-V模型),其模型方程为:
-V模型中,x1(t)、x2(t)分别表示再t时刻食饵和捕食者的密度,而x1、x2则分别表示食饵和捕食者在t时刻的种群密度变化量。a10代表食饵种群内自然增长率、a12代表捕食者对食饵捕食的所造成的影响(或者说威胁)、a20代表捕食者种群内部因素带来的数量变化,比如生老病死亦或同族之间的相互争斗带来的影响、a21代表食饵对捕食者所提供的积极影响,即提供食物来源[2]。
显而易见,L-V模型相当的粗略,考虑的因素也很少。因此,在此基础上本文通过增加模型的维数,来加长食物链的长度;通过加入更多系数来模拟捕食者捕食时间、捕食者成长期等因素的影响。试图通过考虑到更多的因素来提高模型的真实性,然后可以人为的调整控制因素来使生物链达到一个动态平衡,即达到一个良性。
1.建立系统模型
自然界中,广泛存在着种群之间的捕食与被捕食关系,他们一起构成了生物链,为流传的“大鱼吃小鱼,小鱼吃虾米”就通俗易懂的道出了生物链之间的关系。而为了保护生态链的完整,甚至保证整个生态系统的有序发展,对其中捕食系统的研究就存在相当的意义。
根据实际情况,并借鉴文献[1]我们这里设置的捕食模型应该具备如下前提条件:
该双时滞的四维捕食模型为:
在该模型中,x1(t)、x2(t)、x3(t)、x4(t)分别代表一二三级食饵,和食物链顶端的捕食者。其中,一级食饵为食物链最低端,会被二级食饵捕食,二三级食饵以此类推(即生物链的概念,但捕食者只捕食比自己低一级的食饵,这里不考虑跨级捕食)。
而dx1/dt、dx2/dt、dx3/dt、dx4/dt分别代表一二三级食饵,和食物链顶端的捕食者在时刻的数量变化情况。
1.1参数设置
①对一级食饵,即x1(t)而言:
②对二级食饵,即x2(t)而言:
③对三级食饵,即x3(t)而言
④对最高级捕食者,即x4(t)而言
⑤而对于时间常数τ
1.2约束条件
2.正平衡点及稳定性分析
求正平衡点,即各个种群之间的发展达到一个动态平衡,也就是各个种群的数量达到一个稳态值,变化量为0,即
2.1初始条件下的正平衡点
我们设平衡点为E*由于各个种群的变化率为0,由上文中的模型可得:
求解该方程组,我们可以用到solve函数来求解,其对应程序如下:
%计算正平衡点
syms x1 x2 x3 x4%定义自变量
syms a10 a11 a12 a20 a22 a21 a23 a30 a33 a32 a34 a40 a43%定义系数
[solx1,solx2,solx3,solx4]=solve(x1*(a10-a11*x1-a12*x2)==0,x2*(-a20-a22*x2+a21*x1-a23*x3)==0,x3*(-a30-a33*x3+a32*x2-a34*x4)==0,x4*(-a40+a43*x3)==0,x1,x2,x3,x4)
solutions=[solx1,solx2,solx3,solx4]%方程组
m=double(solx1)%x1*的值
y=double(solx2)%x2*的值
z=double(solx3)%x3*的值
h=double(solx4)%x4*的值
由于x1*、x2*、x3*、x4*均为非负,因此为了表达的简洁,这里只把符合条件的解列出来:
可以看出,符合条件的解就只有一组,因此这也佐证了该系统只有唯一的一个正平衡点E*
2.2特征矩阵求解
利用solve函数可以求得模型(1)的雅可比矩阵,其程序为:
%计算正平衡点
syms x1 x2 x3 x4%定义自变量
x=[x1,x2,x3,x4]
syms a10 a11 a12 a20 a22 a21 a23 a30 a33 a32 a34 a40 a43%定义系数
f=[a10*x1-a11*x1^2-a12*x1*x2,-a20*x2-a22*x2^2+a21*x1*x2-a23*x2*x3,-a30*x3-a33*x3^2+a32*x2*x3-a34*x3*x4,-a40*x4+a43*x3*x4]%定义函数,以矩阵形式展示
j=jacobian(f,x)%计算雅可比矩阵
但上述公式并未考虑到时滞因素τ,加上时滞因素后,求得模型(1)的完整雅可比矩阵为:
由雅可比矩阵可以推出其对应特征矩阵X为:
通过特征矩阵X我们又可以推出其特征方程(其中为了在MATLAB里表示方便,τ用y代替),这里可以用到det和collect函数,其具体程序为:
%求特征方程
syms x1 x2 x3 x4%定义自变量
syms a10 a11 a12 a20 a22 a21 a23 a30 a33 a32 a34 a40 a43 y t1 t2 e%定义系数
f=[-a11*x1,-a12*x1*e^(-y*t1),0,0;a21*x2*e^(-y*t2),-a22*x2,-a23*x2*e^(-y*t1),0;0,a32*x3*e^(-y*t2),-a33*x3,-a34*x3*e^(-y*t1);0,0,a43*x4*e^(-y*t2),0];
m=det(f)%计算矩阵
k=collect(m,y)%得特征方程
根据初始条件τ=τ1+τ2,将计算结果整合可得出其特征方程为:
其中:
2.3稳定性分析
2.3.1 τ1=τ2=τ=0时
特征方程(2)可以改写为:
根据文献[3]对正平衡点的初值情况分析思路,我们可得以劳斯-赫尔维茨判据为依据,对于四阶方程来说,为了保证特征方程M的系统稳定,即本文的系统(1)。其需要满足假设H1:
2.3.2存在一个临界τ0,使得τ2=τ0-τ1
即假设存在一个临界滞量τ0,使得当τ在τ0附近取值时(但τ<τ0),系统(1)仍稳定。这里就近似的把系统(1)转化为一个单时滞系统。
τ>0,假设λ=iw(w>0)作为特征方程(2)的一个纯虚根带入(2)中,分离出实部和虚部,得:
根据上式,整理出sinwτ和coswτ的表达式:
其中:
2.4 总结
3.仿真
3.1参数设置
3.2系数求解
在正平衡点E*处,将参数带入模型(1)中得:
求解上式方程组,对应程序为:
syms S L A M
[solS,solL,solA,solM]=solve(S*(1.4-0.2*S-0.4*L)==0,L*(-0.2-0.2*L+0.3*S-0.4*A)==0,A*(-0.2-0.2*A+0.3*L-0.4*M)==0,M*(-0.2+0.3*A)==0,S,L,A,M)%列出方程组
solutions=[solS,solL,solA,solM]
x1=double(solS)%x1*的值
x2=double(solL)%x2*的值
x3=double(solA)%x3*的值
x4=double(solM)%x4*的值
得:
从上述众多组解中,由于x*均大于0,则取最后一组解(并化为小数形式),即:
3.3结果仿真
程序:
clear;clc;
%--------------------------------------------------------------------
% 参数设置
%--------------------------------------------------------------------
a10=1.4;
a11=0.2;
a12=0.4;
a20=0.2;
a22=0.2;
a21=0.3;
a23=0.4;
a30=0.2;
a33=0.2;
a32=0.3;
a34=0.4;
a40=0.2;
a43=0.3;
%初值取平衡点E*
S=2.9197;
L=2.0417;
A=0.6667;
M=0.6979;
T=0:60;
%tao取值
t1=10;
t2=8;
t0=t1+t2;
for idx = 1:length(T)-1
if idx<=t0%第一阶段
S(idx+1)=S(idx)+S(idx)*(a10-a11*S(idx)-a12*L(idx));
L(idx+1)=L(idx)+L(idx)*(-a20-a22*L(idx)+a21*S(idx)-a23*A(idx));
A(idx+1)=A(idx)+A(idx)*(-a30-a33*A(idx)+a32*L(idx)-a34*M(idx));
M(idx+1)=M(idx)+M(idx)*(-a40+a43*A(idx));
elseif idx>t0%第二阶段,时滞开始起作用
S(idx+1)=S(idx)+S(idx)*(a10-a11*S(idx)-a12*L(idx-t1));
L(idx+1)=L(idx)+L(idx)*(-a20-a22*L(idx)+a21*S(idx-t2)-a23*A(idx-t1));
A(idx+1)=A(idx)+A(idx)*(-a30-a33*A(idx)+a32*L(idx-t2)-a34*M(idx-t1));
M(idx+1)=M(idx)+M(idx)*(-a40+a43*A(idx-t2));
end
end
%画出图像
figure
subplot(2,2,1),plot(T,S);
grid on;
xlabel('t');
ylabel('S(t)');
title('x1');
subplot(2,2,2),plot(T,L);
grid on;
xlabel('t');
ylabel('L(t)');
title('x2');
subplot(2,2,3),plot(T,A);
grid on;
xlabel('t');
ylabel('A(t)');
title('x3');
subplot(2,2,4),plot(T,M);
grid on;
xlabel('t');
ylabel('M(t)');
title('x4');
4结论
5参考文献
- 武利青,王晓云.一类具有双时滞四维捕食模型的定性分析.中北大学学报:自然科学版,2019(2):97-102.
- 李大虎,陈淑苹,童欢,朱文文.具有捕食者相互残杀项双时滞系统的Hopf分支.湖北师范学院学报:自然科学版,2011,31(3):76-80.
- 李颖. HR神经元模型的正平衡点稳定性分析.甘肃科技纵横,2017.01.022
- Yang, LX Yang, XF. Propagation Behavior of Virus Codes in the Situation That Infected Computers Are Connected to the Internet with Positive Probability.[J].Discrete Dynamics in Nature and Society,2012,(1):1-13
标签:x4,idx,食饵,笔记,matlab,x2,x3,x1,时滞 来源: https://blog.csdn.net/L1234X/article/details/120661320