具有猎物避难所和Holling II 型功能反应函数的离散捕食者-食饵模型的分岔分析*

张佩雪，牛利娟，庞茹一

（西安工程大学理学院，陕西 西安 710048）

0 引言

捕食者和食饵之间动态关系的普遍性和重要性,使得捕食者-食饵模型长期以来一直是生态学和数学生态学的主要课题之一[1].如果捕食者或者天敌的数量过多,食饵有可能灭绝.避难所对捕食者-食饵的相互作用具有稳定作用.Feller[2]建立了具有猎物避难所的捕食者-食饵模型.Bick 等[3]研究发现避难所的存在可以保护食饵不被灭绝.近年来,离散捕食者-食饵模型受到许多学者的关注[4].物种间相互作用的微分方程模型是数学在生物学中的经典应用之一.马兆芝等[5]研究发现,避难所大小是影响动力学行为的一个关键性因素.Rodriguez等[6]研究了在避难区,捕食者无法捕食猎物.Teng[7]研究了捕食者-食饵模型的复杂动力学行为.Fu 等[8]讨论了Holling II 型捕食者-食饵模型的动力学性质.Hu 等[9]研究了捕食者-食饵模型的Flip 分岔和Neimark-Sacker 分岔.程利芳[10]讨论了分岔理论在生态模型中的应用.Khan 等[11]借助中心流形定理和分岔理论研究了离散捕食者-食饵模型的分岔.

介绍离散模型之前,我们引入具有猎物避难所和Holling II 型功能反应函数的连续捕食者-食饵模型[12]:

其中:x,y分别表示在t时刻的食饵和捕食者种群密度,a,k,α,β,γ,c都是正常数.这里α代表食饵的内在生长速率;k代表食饵的携带能力;γ是捕食者的死亡率;β/α是单位时间内每个捕食者可吃掉的最大食饵数量; 1/a是达到该比率一半所需的食饵密度;c是表示每个捕获的食饵的新生捕食者数量的转换因子;m表示食饵避难率,βx/(1+ax)表示Holling II 型[13]功能反应函数.

下面我们用欧拉向前方法离散化模型(1):

本文讨论了模型(2)平衡点的存在性,利用中心流形定理和分岔理论给出模型(2)的Flip 分岔和Neimark-Sacker 分岔的存在条件.通过数值模拟发现模型(2)出现周期解、极限环、稳定周期解、稳定极限环、不稳定极限环.当参数在某个特定值时,模型(2)的解会出现混沌现象.螨类捕食者-食饵的相互作用通常表现出空间避难所,这为食饵提供了一定程度的保护,使其免受捕食,减少了因捕食而灭绝的机会,此模型对实际生活和理论指导有一定的参考价值.

1 平衡点分析

很容易看出E0(0,0)和E1(k,0)是模型(2)的两个平衡点.

此外,当满足条件cβ-γa＞0 和0≤m≤1-γ/(k(cβ-γa))时,E2(x2,y2)存在正平衡点,x2=γ/k(cβ-γa)(1-m),y2=ac(k(cβ-γa)(1-m))/k(cβ-γa)(1-m)2.

下面讨论正平衡点E2(x2,y2)是否存在Flip 分岔和Neimark-Sacker 分岔.记在E2(x2,y2)的雅可比矩阵为

J(E2)的相应特征方程可以写成F(λ)=λ2-trJ(E2)λ+detJ(E2).其中,

显然F(1)＞0,当α=α1=[4k(1+a(1-m)x2)2+2kaβ(1-m)2x2y2+kγβ(1-m)y2]/2x2(1+a(1-m)x2)2,此时F(-1)=0,trJ(E2)≠0.当trJ(E2)=0 时,有α=α*=[4k(1+a(1-m)x2)2+kaβ(1-m)2x2y2]/2x2(1+a(1-m)x2)2.则在E2可能存在Flip 分岔,定义M={(a,m,k,α,β) :α=α1,α≠α*,a＞0,m＞0,k＞0,γ＞0,β＞0}.当detJ(E2)=1 时,存在共轭复根,且|λ1|=|λ2|=1.c=c1为(cβ-γa)2/c(k(cβ-γa)(1-m)-γ)-[x2aβ(1-m)+γβ]/x2(1-m)(1+a(1-m)x2)2=0 的解.则在E2可能存在Neimark-Sacker 分岔,定义N={(a,m,k,α,β):c=c1,a＞0,m＞0,k＞0,γ＞0,β＞0}.

2 Flip 分岔

本节应用中心流形定理[14]研究了Flip 分岔的存在性.选择α作为分岔参数,给参数α一个小扰动，并且令un=xn-x2,vn=yn-y2,模型(2)变为

将模型(3)在(un,vn,)=(0,0,0)时进行泰勒级数展开至二阶

令

由于|λ1|=-1,|λ2|≠-1,通过计算,得到矩阵J(E2)的特征值分别为λ1=-1,λ2=3-αx2/k+αβ(1-m)2x2y2/(1+a(1-m)x2)2.

设矩阵

使用翻转

则模型(3)变成以下形式

对模型(5)应用中心流形定理.假设Wc(0,0,0)表示模型(5)在=0 的小邻域内的中心流形,则Wc(0,0,0)计算为

此外,对于限制于中心流形Wc(0,0,0)的映射,我们考虑映射G*

在(un,vn,)=(0,0,0)时定义两个判别量l1≠0 和l2≠0.

因此,根据上述分析和文献[14]中的定理3.1,我们给出了在正平衡点Flip 分岔存在的条件.

定理1如果l2≠0,模型(2)在E2(x2,y2)处存在Flip 分岔.如果l2＞0,则从E2(x2,y2)分岔的周期2 点是稳定的; 如果l2＜0,则从E2(x2,y2) 分岔的周期2 点是不稳定的.

3 Neimark-Sacker 分岔

本节选择了c作为分岔参数来研究在E2(x2,y2)的Neimark-Sacker 分岔.令c=c1+,un=xn-x2,vn=yn-y2,其中||≪1,模型(2)变为

将模型(7)在(un,vn)=(0,0)时进行泰勒级数展开至三阶

其中c11,c12,···,c15,c21,c22,···,c25同模型(4)中的a11,a12,···,a15,a21,a22,···,a25,

模型(8)在原点(0,0)处线性化的特征方程为λ2-P()λ+Q(),其中,

使用翻转

则模型(7)变成以下形式

其中,

定义第一Lyapunov 指数如下所示[14]

其中,

因此,根据上述分析和文献[15]中的Neimark-Sacker 分岔存在理论,我们得到定理2,表明Neimark-Sacker分岔存在和方向的参数条件.

定理2如果满足非退化线性条件且L≠0,模型(2)在E2(x2,y2)处存在Neimark-Sacker 分岔.若L＜0(或者L＞0),则当＞0(或者＜0) 时,吸引(或者排斥) 不变闭合曲线从E2(x2,y2)分岔.

4 数值模拟

例1选择参数a=0.1,c=4,k=0.4,m=0.5,β=0.3,γ=0.2,通过数值计算得平衡点E2(0.338 983,1.269 07)和α=2.454 24,λ1=-1,λ2=0.926 375.模型(2)具有平衡点E2(x2,y2)和(a,m,k,α,β)∈M.我们只改变参数α以查看模型(2)的动力学行为变化.其中α∈[2.2,3.4].如图1(a)所示,对于平衡点E2(x2,y2)，在α＜2.454 24 时是稳定的,当α=2.454 24 时失去了稳定性,当α＞2.454 24 时,存在Flip 分岔.此外,随着α的增大,出现了一个混沌集合.图2(a)是参数α∈[2.2,3.4] 时模型(2)的相图.

图1 模型(2)的分岔图

图2 模型(2)不同α 值、c 值的相图

例2选择参数a=5,k=0.4,α=3,m=0.1,β=0.4,γ=1.5,通过数值计算得出平衡点E2(0.222 222,0.740 741)和c=3.75,λ1,2=0.499 999±0.866 026i.模型(2)具有平衡点E2(x2,y2)和(a,m,k,α,β)∈N.我们只改变参数c以查看模型(2)的动力学行为变化.其中c∈[3.5,4.2].如图1(b)所示,对于平衡点E2(x2,y2)，在c＜3.75 时是稳定的,当c=3.75 时失去了稳定性,当c＞3.75 时,存在Neimark-Sacker 分岔.此外,随着c的增大,出现了一个混沌集合.图2(b)是参数c∈[3.5,4.2]时模型(2)的相图,有周期解和极限环出现.