/*仅当作学习笔记，若有纰漏欢迎友好交流指正，此外若能提供一点帮助将会十分荣幸*/

本系列谈论过单时滞，但还没提及过双时滞，本文将着重介绍一种双时滞系统并对其进行简单处理分析。

摘要：本文针对一个捕食网络，根据其特性建立起一个双时滞四维系统。首先分析其稳定性，找到其正平衡点，然后通过调试时滞量τ来模拟预测网络模型的发展前景。并根据所求得临界滞量τ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)-1if 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

双时滞四维捕食网络的分析【基于matlab的动力学模型学习笔记_6】相关推荐

带时滞的病毒模型计算模板【基于matlab的动力学模型学习笔记_1】
/*仅当作学习笔记,若有纰漏欢迎友好交流指正,此外若能提供一点帮助将会十分荣幸*/ 摘要:无论是生物病毒还是网络病毒,其内核的传播机理都有很多的相似之处.因此,本文在经典的SIR病毒模型基础上改造出 ...
一维离散动力学系统的混沌研究【基于matlab的动力学模型学习笔记_8】
摘要:混沌(Chaos)是指发生在确定系统中的貌似随机的不规则运动,本文将基于几种经典的一维动力学方程系统,根据其动力学方程研究其混沌产生过程以及相对应的MATLAB仿真. /*仅当作学习笔记,若有 ...
二维离散动力学系统的混沌研究【基于matlab的动力学模型学习笔记_9】
摘要:混沌(Chaos)是指发生在确定系统中的貌似随机的不规则运动,本文将基于经典的二维系统,然后根据动力学方程研究其混沌产生过程以及相对应的MATLAB仿真,再讨论Lyapunov指数以及正平衡点 ...
基于传染病模型中的再生数R0的讨论【基于matlab的动力学模型学习笔记_2】
/*仅当作学习笔记,若有纰漏欢迎友好交流指正,此外若能提供一点帮助将会十分荣幸*/ 在上一篇博文中介绍了病毒模型的基本计算思路方法,而本文将会重点讨论基本再生数R0-这个决定病毒是继续发展还是衰减的关 ...
基于Duffing系统的分数阶混沌研究【基于matlab的动力学模型学习笔记_5】
/*仅当作学习笔记,若有纰漏欢迎友好交流指正,此外若能提供一点帮助将会十分荣幸*/ 前面的几篇博文我们提到提到的都是整数阶模型,这里我们将对分数阶模型进行一个简单的研究. 摘要:与整数阶混沌相比,分数 ...
lyapunov指数求取时运用qr法与jacobi法之间的区别与联系【基于matlab的动力学模型学习笔记_10】
在进行lyapunov指数的求取时,需要知道离散动力学系统对应Jacobi矩阵的特征值,qr法与Jacobi法都可以求解矩阵特征值,其中qr法求解的是矩阵所有特征值,而Jacobi法求解的是矩阵的最大 ...
《对冲基金建模与分析基于MATLAB》简介及PDF下载
转 <对冲基金建模与分析--基于MATLAB>简介及PDF下载内容简介本书是关于用MATLAB对对冲基金进行建模和分析的入门读物.在对对冲基金的基本概念.分类.相关工具和指标系统介绍的 ...
计算机基础与程序设计（基于C语言）学习笔记
计算机基础与程序设计(基于C语言)学习笔记前言这是一个学习笔记课程导入在线学习工具:https://c.runoob.com/compile/11 为什么要学习程序设计 (1)存储程序和程序控 ...
【毕业设计之python系列】基于Flask的在线学习笔记的设计与实现
基于Flask的在线学习笔记的设计与实现摘要在线学习笔记系统是一种为学生和教师提供在线学习和教学的平台.本文基于Flask框架,设计并实现了一个在线学习笔记系统.该系统支持用户注册.登录.创建课程 ...

双时滞四维捕食网络的分析【基于matlab的动力学模型学习笔记_6】

0引言