具有合作狩猎的食物链模型的Hopf分支.pdf

浙江大学学报（理学版）Journal of Zhejiang University（Science Edition）http：/ 50 卷第 5 期2023 年 9月Vol.50 No.5Sept.2023 具有合作狩猎的食物链模型的 Hopf分支韩婉琴，石垚*，包雄雄（长安大学 理学院，陕西 西安 710064）摘要：为研究捕食者间合作行为对食物链系统动力学行为的影响，在原有合作捕食模型基础上引入顶层捕食者（仅以捕食者种群为食），建立了具有合作狩猎的食物链模型。运用线性化方法讨论了平衡点的局部稳定性，并通过构造合适的 Lyapunov函数，给出了系统全局稳定的充分条件。利用中心流形约简定理导出了分支周期解稳定性的显式公式，并通过数值模拟验证了理论分析结果。结果显示，当捕食者间无合作时，正平衡点为稳定焦点，随着合作捕食参数的增大，系统出现稳定的极限环且随的增大不断胀大。说明如果捕食者种群密度过大，系统将产生持续的周期振荡，即食饵、中级捕食者、顶层捕食者要么以周期振荡的形式共存，要么种群数最终趋于稳定。因此，合作狩猎更有利于维护生态平衡。关键词：食物链；合作捕食；稳定性；Hopf分支；数值模拟中图分类号：O 175.26 文献标志码：A 文章编号：10089497（2023）0553912HAN Wanqin,SHI Yao,BAO Xiongxiong（College of Science，Changan University，Xi an 710064，China）Hopf bifurcation of a food chain model with cooperative hunting.Journal of Zhejiang University(Science Edition),2023,50(5):539550,557Abstract:In order to study the impact of cooperative behaviors between predators on the dynamic behavior of the food chain system,a top-level predator(which only feeds on the predator population)on the basis of the original cooperative hunting model is introduced,and a food chain model with cooperative hunting is established.The local stability of the equilibrium point is discussed by using the linearization method.By establishing an appropriate Lyapunov function,a sufficient condition for the global stability of the system at the equilibrium point is given.In addition,by using the central manifold reduction theorem,the explicit formulas of the branch periodic solution is studied.Numerical simulations are conducted to verify our theoretical analysis.The results show that when there is no cooperative relationship between predators,the positive equilibrium point is stable.With the increase of cooperative predation parameter,the system will have a stable limit cycle,and the limit cycle will continue to expand following the increase of cooperative predation parameter.Moreover,if the predator population density is too large,the system will produce continuous periodic oscillation,that is,the three species i.e.diet,mid-level predator,top-level predator,either coexist in the form of periodic oscillation,or the number of species eventually tends to be stable.Therefore,cooperative hunting is more conducive to maintaining ecological balance.Key Words:food chain;cooperative hunting;stability;Hopf bifurcation;numerical simulationDOI：10.3785/j.issn.1008-9497.2023.05.004收稿日期：20230215；修回日期：20230425；接受日期：20230505；出版日期：20230925.基金项目：国家自然科学基金资助项目（12201067，12101075）；陕西省自然科学基础研究计划项目（2021JQ-217，2022JQ-054）；中央高校基本科研业务费专项资金资助项目（300102122114，300102122104）.作者简介：韩婉琴（1999），ORCID：https：/orcid.org/0009-0009-4805-7958，女，硕士研究生，主要从事微分方程动力学及应用研究，E-mail：.*通信作者，ORCID：https：/orcid.org/0000-0003-1442-6646，E-mail：.浙 江 大 学 学 报（理学版）第 50 卷0引 言在自然界中，各种群不是孤立存在的，均与其他种群存在各种各样的关系，如竞争、合作、互惠等关系。各种群能够繁衍生息的关键在于彼此之间存在合作行为。研究表明，不少动物在捕食过程中为了提高攻击率和觅食便利化，会采用集体狩猎的方式，如狮子1、黑猩猩2、野犬3等。合作有利于物种的持续生存，保持生态系统的多样性。因此研究种群间的合作捕食强度对于促进生态系统的可持续发展有重要意义。本文基于 ALVES 等4提出的合作狩猎函数，研究三种群食物链模型：dudt=u()1-uk-(1+v)uv，dvdt=(1+v)uv-m1v-wv，dwdt=wv-m2w，（1）其中，u、v、w 分别表示食饵、捕食者、顶层捕食者的种群密度，u呈 Logistic增长，为食饵的出生率，k为食饵的环境容纳量，m1、m2分别为捕食者与顶层捕食者的自然死亡率，为合作捕食参数，可反映捕食者间的合作狩猎关系，=0表示捕食者间无合作关系，而越大（小），表示捕食者间的合作越强（弱）。所有参数均为正。ALVES等4研究了平衡点和分支的稳定性，结果表明，与无狩猎合作的情况相比，合作捕食有利于调节捕食者的存活率并产生振荡现象。ZHANG等5同时考虑了具有捕食者间狩猎合作和 Allee 效应的平面三次微分系统，发现合作捕食对降低捕食者灭绝的风险有一定的指导作用。对于具有狩猎合作和 Holling-II 型功能反应函数的反应扩散捕食-食饵系统的动力学问题，FU 等6讨论了非负平衡态的稳定性、Hopf分支、鞍结点分支和 B-T分支，结果表明，当合作捕食参数等于某阈值时，会出现超临界轨道渐近稳定的极限环。此外，狩猎合作会导致系统产生不稳定效应，出现双稳现象以及发生鞍结点分支等复杂动力学问题，与无合作7-12形成鲜明对比。本文在原有合作捕食模型的基础上引入顶层捕食者，使式（1）始终存在一个正平衡点，且在正平衡点出现 Hopf分支。此外，捕食者之间合作狩猎会影响种群的共存状态。当捕食者间无合作关系时，正平衡点为稳定的焦点，随着合作捕食参数的增大，系统出现一个稳定的极限环，且极限环随合作捕食参数的增大而胀大，这与无合作狩猎的情况形成了鲜明的对比。同时，如果捕食者的种群密度过大，系统将产生持续的周期振荡，即三个种群要么以周期振荡的形式共存，要么种群数量最终趋于稳定。因此，合作狩猎更有利于维护生态平衡。但在实际种群生态环境中，个体的随机运动是影响种群间动力学行为的重要因素，研究扩散对动力学行为的影响非常有意义。1平衡点的稳定性和存在性根据比较原理13，得到式（1）正解的先验估计。引 理 1 R3+=(u，v，w)|u 0，v 0，w 0是式（1）的正向不变集。证明 由式（1），可得 u(t)=u(0)exp0s ()1-uk-(1+v)v ds，v(t)=v(0)exp0s(1+v)u-m1-w ds，w(t)=w(0)exp0s(v-m2)ds，（2）所 以，若(u(0)，v(0)，w(0)0，则(u(t)，v()t，w()t)0，t0，+)。证毕。下面，在初始条件u 0，v 0，w 0下，得到系统解的最终有界性。引理 2 如果u 0，v 0，w 0，则式（1）的解非负且一致最终有界。证明 式（1）解的非负性可由引理 1 推证。下证式（1）解一致最终有界。令P(u(t)，v(t)，w(t)=u(t)+v(t)+w(t)，则dPdt=dudt+dvdt+dwdt=u()1-uk-m1v-m2w。设0 0，且是式（1）的正向不变集，吸引R3+中的所有正解，即式（1）的解一致最终有界，证毕。引 理 3 如 果(u，v，w)是 式（1）的 任 意 非 负解，则lim supt u(t)k，lim supt v(t)M，lim supt w(t)L。证 明 由 式（1）的 第 1 个 方 程，可 得dudtu()1-uk，进而lim supt u(t)k。故对任意的1 0，存在正实数T1 0，使得当tT1时，u(t)k+1。对任意的t T1，有d(u+v)dt=u()1-uk-m1v-wv(k+1)-m1v G-m1(u+v)，其中，G=(+m1)(k+1)，则lim supt u(t)+v(t)Gm1。因此存在一个充分大的正数M 0，使得lim supt v(t)M，所以对于任意的2 0，存在T2 T1，使得当t T2时，v(t)M+2。同理，对任意的t T2，有d(v+w)dt=(1+v)uv-m1v-m2w(k+1)1+(M+2)(M+2)-m2w H-m2(v+w)，其 中，H=(k+1)1+(M+2)+m2(M+2)，从而lim supt(v(t)+w(t)Hm2。因此存在充分大的正数L 0，使得lim supt w(t)L，所以对任意的3 0，存在T3 T2，使得当t T3时，w(t)L+3。证毕。下面讨论式（1）平衡点的存在性。引入无顶端捕食者系统，即w=0，模型为 ()1-uk-(1+v)v=0，(1+v)u-m1=0。（3）定理 1 设式（1）参数均大于 0，其中为食饵的出生率，则式（1）存在灭绝平衡点E0(0，0，0)和无捕食者存在的半平凡平衡点Ek(k，0，0)。当w 0时，式（1）存在唯一的共存平衡点E(u，v，w)。当w=0时，有 3种情况：（i）令1=(，k，m1)R4+：3 0，则式（3）在1中有唯一的正平衡态E*1(u*1，v*1，0)；（ii）令2=(，k，m1)R4+：20，则 式（3）在2中 有 2 个 内 部 平 衡 点，即E*i(u*i，v*i，0)(i=2，3)，或者无内部平衡点；（iii）令3=(，k，m1)R4+：20，30，则式（3）在3中无内部平衡点；其中，2=k(1-)，3=(m1-k)。证明 式（1）可整理为 dudt=u ()1-uk-(1+v)v，dvdt=v(1+v)u-m1-w，dwdt=w(v-m2)，（4）为讨论式（1）的平衡态，令式（4）3 个方程的右边等于 0：u ()1-uk-(1+v)v=0，v(1+v)u-m1-w=0，w(v-m2)=0，（5）显然，式（1）存在灭绝平衡点E0(0，0，0)，以及无捕食者存在的半平凡平衡点Ek(k，0，0)。下面分 2 种情况讨论：（1）w 0，即存在顶层捕食者；（2）w=0，即不存在顶层捕食者。当w 0时，式（1）存 在 唯 一 的 共 存 平 衡 点E(u，v，w)，其中u，v，w满足：()1-uk-(1+v)v=0，(1+v)u-m1-w=0，v-m2=0，（6）可得541浙 江 大 学 学 报（理学版）第 50 卷u=(-m2-m22)k，v=m2，w=k(-m2-m22)(1+m2)-m1。令(u*i，v*i)为式（3）的正解，由式（3）的第 2个方程，可得u*i=m11+v*i，代入式（3）的第 1个方程，可得0v*i3+1v*i2+2v*i+3=0，（7）其中，0=k20，1=2k0，2=k(1-)，3=(m1-k)。由于2，3的符号不确定，所以考虑 3个区域：1=(，k，m1)R4+：3 0，2=(，k，m1)R4+：2 0，3=(，k，m1)R4+：2 0，3 0。由笛卡尔符号法则14，可知：在区域1中，式（7）的符号 改 变 了 一 次，从 而 式（3）有 唯 一 的 正 平 衡 态E*1(u*1，v*1，0)；在区域2中，式（7）的符号改变了 2次，因此式（3）有 2 个内部平衡点，即E*i(u*i，v*i，0)(i=2，3)，或者无内部平衡点；在区域3中，式（7）的符号没有改变，从而式（3）无内部平衡点。证毕。定理 2 设式（1）的参数均大于 0，则（i）式（1）在平衡点E0(0，0，0)处是不稳定的；（ii）当k k时，式（1）在平衡点Ek(k，0，0)全局渐近稳定；（iii）令1=-()1-2uk+(1+v)v-u(1+2v)+m1，2=()1-2uk-(1+v)v u(1+2v)-m1+v(v-m2)+u(1+2v)(1+v)v，3=-()1-2uk-(1+v)vv(v-m2)，则当1 0，3 0且12 3时，式（1）在平衡点E(u，v，w)局部渐近稳定，当()1-uk(1+v)v，(1+v)um1+w，vm2，u(1-uk)m1v+m2w时，式（1）在平衡点E(u，v，w)全局渐近稳定。（iv）令c1=()1-2u*ik-(1+v*i)v*i+u*i(1+2v*i)-m1，c2=()1-2u*ik-(1+v*i)v*i u*i(1+2v*i)-m1+u*i(1+2v*i)(1+v*i)v*i，则 当c1 0时，式（3）在 平 衡 点E*i(u*i，v*i，0)(i=2，3)局部渐近稳定，当k4+(1+v)uv(1+v)v，(1+v)um1 时，式（3）在 平 衡 点E*i(u*i，v*i，0)(i=2，3)全 局 渐 近稳定。证明 首先，证明式（1）在平衡点E0(0，0，0)，Ek(k，0，0)，E(u，v，w)的稳定性。式（1）在平衡点E0(0，0，0)的 Jacobian矩阵为J(E0)=()000-m1000-m2，其特征值为，-m1，-m2，因为 Jacobian 矩阵的特征值 0，所以式（1）在平衡点E0(0，0，0)不稳定。式（1）在平衡点Ek(k，0，0)的 Jacobian矩阵为J(Ek)=()-k00k-m1000-m2，当k 0，3 0且12 3时，式（1）在 平 衡 点E(u，v，w)局 部 渐 近稳定。再证式（1）在平衡点Ek(k，0，0)和E(u，v，w)的全局稳定性。在平衡点Ek(k，0，0)，构造 Lyapunov函数L1=u+v+w，其沿式（1）轨线的全导数为dL1dt=dudt+dvdt+dwdt=u()1-uk-m1v-m2w，显 然，当u k时，对 任 意 的(u，v，w)0，均 有dL1dt 0，且dL1dt=0当且仅当(u，v，w)=(k，0，0)。由 Lyapunov-Lasalle 不变集原理16，当u k时，式（1）在平衡点Ek(k，0，0)全局渐近稳定。在平衡点E(u，v，w)，构造 Lyapunov函数L2=()u-u-ulnuu+()v-v-vlnvv+()w-w-wlnww，其沿式（1）轨线的全导数为dL2dt=()1-uududt+()1-vvdvdt+()1-wwdwdt=(u-u)()1-uk+(uv-uv)(1+v)-m1(v-v)-m2(w-w)-(wv-vw)=u()1-uk-u()1-uk+uv(1+v)-uv(1+v)-m1v+m1v-m2w+m2w-wv+vw=-()1-uk-(1+v)v u-(1+v)u-m1-w v-(v-m2)w+u()1-uk-m1v-m2w，显 然，当()1-uk(1+v)v，(1+v)u m1+w，v m2，u()1-uk m1v+m2w时，对 任 意 的(u，v，w)0，有dL2dt 0，且dL2dt=0当 且 仅 当(u，v，w)=(u，v，w)。根据 Lyapunov-Lasalle不变集原理16，式（1）在平衡点E(u，v，w)全局渐近稳定。式（3）在共存平衡点E*i的 Jacobian矩阵为J(E*i)=()1-2u*ik-(1+v*i)v*i-u*i(1+2v*i)(1+v*i)v*iu*i(1+2v*i)-m1，其对应的特征方程为2-c1+c2=0，（9）其中，c1=tr(J(E*i)=()1-2u*ik-(1+v*i)v*i+u*i(1+2v*i)-m1，c2=det(J(E*i)=()1-2u*ik-(1+v*i)v*i u*i(1+2v*i)-m1+u*i(1+2v*i)(1+v*i)v*i。由 Routh-Hurwitz 判据15，可得当c1 0时，式（3）在平衡点E*i(u*i，v*i，0)局部渐近稳定。式（6）在平衡点E*i(u*i，v*i，0)的全局稳定性：构造 Lyapunov函数L3=()u-u*i-u*ilnuu*i+()v-v*i-v*ilnvv*i，其沿式（6）轨线的全导数为dL3dt=()1-u*iududt+()1-v*ivdvdt=u()1-uk-u*i()1-uk+(1+v)u*iv+(1+v)uv-m1v-(1+v)uv*i+m1v*ik4+(1+v)uv-m1v-1-uk-(1+v)v u*i-(1+v)u-m1 v*i，显 然，如 果k4+(1+v)uv(1+v)v，(1+v)u m1，则dL3dt 0，当 且 仅 当(u，v，w)=(u*i，v*i，0)时 取 等 号。由 Lyapunov-Lasalle 不变集原理16，式（3）在平衡点E*i(u*i，v*i，0)全局渐近稳定。证毕。注 1 当w=0时，式（1）可约化为式（3），赵叶青等17讨论了式（3）平衡点的局部稳定性。基于上述工作，本文通过构造合适的 Lyapunov 函数，分析了式（3）在平衡点E*i(u*i，v*i，0)全局渐近稳定的充分条件。543浙 江 大 学 学 报（理学版）第 50 卷2Hopf 分支为深入探究合作捕食对式（1）的影响，选取合作捕食参数作为分支参数，进一步建立式（1）在正平衡点E(u，v，w)产生的 Hopf分支。定理 3 设j()=j()+i j()为式（1）在E(u，v，w)处的特征方程式（8）的特征值，当=h时，有j(h)=0，j(h)0，则当3=12时，如果djd=h 0，那么式（1）在正平衡点E(u，v，w)产生Hopf分支。证明 式（1）在E(u，v，w)处的特征方程为3+12+2+3=0。当3=12时，有(2+2)(+1)=0，（10）显然，式（10）有一对纯虚根1，2=2i和一个负根3=-1。设j()=j()+i j()为式（8）的特征值，且当=h时，有j(h)=0，j(h)0，将j()代入式（8），求导，分离实部和虚部，可得J()djd|=h-A()djd|=h+C()=0，A()djd|=h+J()djd|=h+D()=0，其中，J()=32j()+21()j()+2()-32j()，A()=6j()j()+21()j()，C()=2j()1()+2()j()+3()-1()2j()，D()=2j()j()1()+2()j()。下面验证横截条件djd|=h 0。因为A()D()+J()C()=1221()21()1()+62()1()21()+41()1()21()+21()2()21()+31()41()-31()21()21()+322()31()+33()21()+21()1()31()-21()1()1()21()+21()2()21()+21()3()1()+1()2()21()-1()2()21()+2()2()1()+2()3()-31()21()21()+31()41()-32()1()21()-33()，所以djd|=h=A()D()+J()C()J2()+A2()|=h=21(h)2(h)21(h)-1(h)2(h)21(h)+31(h)41(h)-33()，又因为1，2，3均不为零，所以djd|=h 0。证毕。注 2 由于参数h的符号不易确定，将通过数值模拟给出其算例。下面运用中心流形约简定理和正规型理论，分析 Hopf分支的方向和 Hopf分支周期解的稳定性。将式（1）的正平衡点E(u，v，w)平移至坐标原点。令 U=u-u，V=v-v，W=w-w，代入式（1），可得dUdt=(U+u)()1-U+uk-1+(V+v)(U+u)(V+v)，dVdt=1+(V+v)(U+u)(V+v)-m1(V+v)-(W+w)(V+v)，dWdt=(W+w)(V+v)-m2(W+w)。（11）为简便起见，用 u，v，w分别表示 U，V，W，可得dudt=(u+u)()1-u+uk-1+(v+v)(u+u)(v+v)：=f(u，v，w)dvdt=1+(v+v)(u+u)(v+v)-m1(v+v)-(w+w)(v+v)：=g(u，v，w)dwdt=(w+w)(v+v)-m2(w+w)：=h(u，v，w)，（12）则f(u，v，w)在点(0，0，0)的三阶泰勒级数为f(u，v，w)=f(0，0，0)+fu(0，0，0)u+fv(0，0，0)v+fw(0，0，0)w+12fuu(0，0，0)u2+12fvv(0，0，0)v2+12fww(0，0，0)w2+fuv(0，0，0)uv+fuw(0，0，0)uw+fvw(0，0，0)vw+16fuuu(0，0，0)u3+16fvvv(0，0，0)v3+16fwww(0，0，0)w3+12fuuv(0，0，0)u2v+544韩婉琴，等：具有合作狩猎的食物链模型的 Hopf分支第 5期12fuuw(0，0，0)u2w+fuvw(0，0，0)uvw+12fvvu(0，0，0)v2u+12fvvw(0，0，0)v2w+12fwwu(0，0，0)w2u+12fwwv(0，0，0)w2v。经计算，可得f(0，0，0)=0，fu(0，0，0)=()1-2uk-(1+v)v，fv(0，0，0)=-u(1+2v)，fw(0，0，0)=0，fuu(0，0，0)=-2k，fvv(0，0，0)=-2u，fww(0，0，0)=0，fuv(0，0，0)=-(1+2v)，fuw(0，0，0)=fvw(0，0，0)=0，fuuu(0，0，0)=fvvv(0，0，0)=fwww(0，0，0)=fuuv(0，0，0)=fuuw(0，0，0)=0，fuvw(0，0，0)=fvvu(0，0，0)=fvvw(0，0，0)=fvvw(0，0，0)=fwwv(0，0，0)=fwwu(0，0，0)=0。g(u，v，w)和h(u，v，w)在点(0，0，0)的泰勒展开与f(u，v，w)类似，此处不再叙述，其中，g(0，0，0)=0，gu(0，0，0)=(1+v)v，gv(0，0，0)=u(1+2v)-m1-w，gw(0，0，0)=-v，guu(0，0，0)=0，g(0，0，0)=0，gu(0，0，0)=(1+v)v，gv(0，0，0)=u(1+2v)-m1-w，gw(0，0，0)=-v，guu(0，0，0)=0，gvv(0，0，0)=2u，gww(0，0，0)=0，guv(0，0，0)=(1+2v)，guw(0，0，0)=0，gvw(0，0，0)=-1，guuu(0，0，0)=gvvv(0，0，0)=gwww(0，0，0)=0，guuv(0，0，0)=guuw(0，0，0)=guvw(0，0，0)=gvvu(0，0，0)=gvvw(0，0，0)=0，gwwv(0，0，0)=gwwu(0，0，0)=0，h(0，0，0)=hu(0，0，0)=huu(0，0，0)=hvv(0，0，0)=0，hww(0，0，0)=huv(0，0，0)=huw(0，0，0)=huuu(0，0，0)=hvvv(0，0，0)=gwww(0，0，0)=0，huuv(0，0，0)=huuw(0，0，0)=huvw(0，0，0)=hvvu(0，0，0)=hvvw(0，0，0)=hwwv(0，0，0)=0，hwwu(0，0，0)=0，hv(0，0，0)=w，hvw(0，0，0)=1，hw(0，0，0)=v-m2。因此式（12）可约化为Z=J(E)Z+F，其中，Z=(u，v，w)T，J(E)=()1-2uk-(1+v)v-u(1+2v)0(1+v)vu(1+2v)-m1-w-v0wv-m2=()fufv0gugvgw0hvhw，F=(g1(u，v，w)，g2(u，v，w)，g3(u，v，w)T，且g1(u，v，w)=-ku2-uv2-(1+2v)uv，g2(u，v，w)=uv2+(1+2v)uv-vw，g3(u，v，w)=vw。设1，3分别为1=2i和3=-1的特征向量，则1=-fv2+fufv2if2u+22i-hv2ihw，3=-fv11+fu1-hv11-hw。令D1=(Im(1)，Re(1)，3)T=fufv2f2u+2-fv2f2u+2-fv11+fu201-hv2hw0-hv11-hw：=p11p12p13p21p22p23p31p32p33，因为矩阵D1非退化，所以D-11=1|D1|q11q12q13q21q22q23q31q32q33，其中，|D1|=fvhv122(f2u+2)hw-fvhv122(f2u+2)(1-hw)，q11=0，q12=2hvhw12-hv212(1-hw)hw，545浙 江 大 学 学 报（理学版）第 50 卷q13=0，q21=hvfv12(f2u+2)(hw-1)，q22=-fufvhv12(f2u+2)(1-hw)-fvhv12(1+fu)hw，q23=fvhv22(f2u+2)hw，q31=-fv12f2u+2，q32=-fufv12f2u+2-fv121+fu，q33=fv22f2u+2。令Z=D1Y，则Y=D-11Z，又 因 为Y=(y1，y2，y3)T，所以Y.=(D-11J(E)D1)Y+U(Y)，（13）其中，U(Y)=D-11F(D1Y)=()m1m2m3=1|D1|q11g1+q12g2+q13g3q21g1+q22g2+q23g3q31g1+q32g2+q33g3，D1-1J(E)D1=0-2020000-1。因此式（13）可表达为y1y2y3=0-2020000-1y1y2y3+()m1m2m3，（14）可得 U1=B1U1+m(U1)，V1=C1V1+g(V1)，（15）其中，U1=(y1，y2)T R2+，V1=y3 R+，B1=()0-220，C1=(-1)，m=(m1，m2)T，g=(m3)，显然，C1的特征值具有负实部且矩阵B1的特征值的实部为零，m为由R2+到R2+的一阶连续可微映射，g为 由R+到R+的 一 阶 连 续 可 微 映 射，m(0)=m(0)=0，g(0)=g(0)=0，因此存在一个不变的中心流形18，可将其表示为Wc(0)=(U1，V1)R2+R|V1=d1(U1)，|U1|1，则d1(U1)=(U1)+O(U1q)。为了逼近中心流形，设y3=d1(y1，y2)=12(b11y21+2b12y1y2+b22y22)+O(y1，y2)，（17）代入式（15），得y3.=d1y1dy1dt+d1y2dy2dt=-1y3+m3。经计算，可得()2b12-12b11y21+()-2b12-12b22y22+(2b22-2b11-1b12)y1y2=Q1y21+Q2y1y2+Q3y22+O(y1，y2)，其中，Q1=1|D1|q31 -kp211-up221-(1+2v)p11p21+q32 up221+(1+2v)(p11p21-p21p31)+q33p21p31，Q2=1|D1|(q31-2kp11p12-2up21p22-(1+2v)(p11p22+p12p21)+q32 up21p22+(1+2v)(p12p21+p11p22)-(p22p31+p21p32)+q33(p22p31+)p21p32)，Q3=1|D1|q31 -2kp212-up222-(1+2v)p22p12 p21+q32 up222+(1+2v)(p12p22-p22p32)+q33p22p32。比较系数，可得546韩婉琴，等：具有合作狩猎的食物链模型的 Hopf分支第 5期Q1=2b12-12b11，Q2=2b22-2b11-1b12，Q3=-2b12-12b22，可解得b11=2(Q1+Q3)+1(2Q2+1Q1)/2-31/4-12，b12=21Q2/4-21(Q3-Q1)/2-31/2-12，b22=21Q3/2-12Q2/2+2(Q1+Q3)-31/4-12。定理 418 式（16）和式（15）具有相同的稳定性，即如果式（16）的零解是稳定（局部渐近稳定或者不稳定）的，那么式（15）的零解也是稳定（局部渐近稳定或者不稳定）的。由式（16），可得()y1y2=()0-22 0()y1y2+()m1m2，（18）其中，m1=1|D1|(q11g1+q12g2+q13g3)，m2=1|D1|(q21g1+q22g2+q23g3)，g1=-k p11y1+p12y2+12p13(b11y21+2b12y1y2+b22y22)2-u p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)2-(1+2v)p11y1+p12y2+12p13(b11y21+2b12y1y2+b22y22)p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)，g2=u p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)2+(1+2v)p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)p11y1+p12y2+12p13(b11y21+2b12y1y2+b22y22)-p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)p31y1+p32y2+12p33(b11y21+2b12y1y2+b22y22)，g3=p21y1+p22y2+12p23(b11y21+2b12y1y2+b22y22)p31y1+p32y2+12p33(b11y21+2b12y1y2+b22y22)。令mkij=mkyiyj|(0，0)，mkijl=mkyiyjyl|(0，0)，则分支周期解的方向和稳定性可由r=d(12-3)d|=h=12 m112(m111+m122)-m212(m211+m222)+m111m211+m122m222-(m1111+m2112+m1122+m2222)确定。综上，可得定理 5 定义 r=d(12-3)d|=h=12 m112(m111+m122)-m212(m211+m222)+m111m211+m122m222(m1111+m2112+m1122+m2222)，当r 0时，在正平衡点产生的 Hopf 分支为超临界分支，当r 0时，该 Hopf分支为亚临界分支。注 3 考 虑 计 算 的 复 杂 性，不 给 出m111，m112，m122，m222，m211，m212，m1111，m2112，m1122，m2222的具体表达式。3数值模拟采用数值模拟的方法验证平衡点的存在性以及边界平衡点、内部平衡点、共存平衡点的稳定性和Hopf 分支、混沌现象，并研究合作捕食强度对种群的影响，进而证实并拓展理论结果。首先，取（u，v，w）的初始值为(0.2，0.1，0.1)，公共参数取值如表 1所示，采用 Matlab软件计算，结果表明，当底层资源不够充足或灭绝时，所有物种均灭绝，见图1（a）。当食饵u以密度k持续生存，中级捕食者v和顶层捕食者w因无法获取足以生存的资源而547浙 江 大 学 学 报（理学版）第 50 卷灭绝，见图1（b）。当食饵u和中级捕食者v稳定共存，中级顶层捕食者w灭绝时，式（6）内部平衡点E*1=(0.049 3，0.029 4，0)，见图 1（c）。当食饵u和中级捕食者v稳定共存，顶层捕食者w灭绝时，式（6）有内部平 衡 点E*2=(0.323 0，0.285 9，0)，E*3=(0.451 4，0.109 8，0)，见图 1（d）。当食饵u、中级捕食者v和顶层捕食者w稳定共存时，经过充分长的时间，三者种群密度趋于共存平衡点E(u，v，w)，见图 1（e）。其次，模拟边界平衡点、内部平衡点及共存平衡点的稳定性，参数取值如表 1所示。在图 2（a）中，中级捕食者和顶层捕食者灭绝，食饵u最终将以密度k持续生存。在图 2（b）中，E*1是稳定的中心焦点，其外围被一族闭轨线族所环绕，说明在E*1附近食饵u与中级捕食者v呈周期变化。图 2（c）表明，E*2和E*3是稳定的结点或焦点，这与文献 18 中的结论一致。在图 2（d）中，食饵、图 1灭绝平衡点、边界平衡点、内部平衡点及共存平衡点的存在性Fig.1Existence of extinction equilibrium，boundary equilibrium，internal equilibrium and coexistence equilibrium points表 1公共参数及取值Table 1Public parameters and values参数取值k=0.5，=0.01，m1=0.01，m2=0.3，=10k=0.5，=6，m1=0.6，m2=0.3，=10k=0.5，=0.03，m1=0.05，m2=0，=0.5k=0.5，=1.5，m1=0.6，m2=0，=3k=0.5，=1，m1=0.1，m2=0.1，=1.2吸引子E0EkE*1E*i(i=2，3)E对应图序号图 1（a）图 1（b）图 1（c）图 1（d）图 1（e）548韩婉琴，等：具有合作狩猎的食物链模型的 Hopf分支第 5期中级捕食者和顶层捕食者三者稳定共存，此时E是稳定的焦点。再次，验证式（1）的 Hopf 分支。取=1，其他参数值同图 1（e），可分别得到u关于、v关于、w关于的分支图，如图 3（a）（c）所示，其中 HB表示Hopf 分 支 点，由 XXPAUT 软 件 得，式（1）在 点E(0.885 2，0.100 0，0.916 0)产生 Hopf 分支，且分支参数h为 1.477。在 1.477 的左侧，正平衡点E局部渐近稳定，在 1.477 的右侧，正平衡点E失去稳定性且出现稳定的周期解，此时由 Hopf分支产生的周期解的周期为 19.170，采用 Matlab 软件计算最大李雅普诺夫指数，结果如图 3（d）所示，发现当取某值时，最大李雅普诺夫指数大于零，系统从稳定态过渡至混沌态。