具有常数捕获猎物和Allee效应的Leslie-Gower模型的动力学研究
The Dynamics Study of the Leslie-Gower Model with Constant Prey Capture and Effect
DOI: 10.12677/aam.2024.137328, PDF, HTML, XML, 下载: 5  浏览: 11 
作者: 孙 营, 孙福芹:天津职业技术师范大学理学院,天津
关键词: Allee效应有界性稳定性平衡点Allee Effect Boundedness Stability Equilibrium Point
摘要: 文章考虑具Allee效应的Leslie-Gower模型。在常数捕获的情况下,研究了捕食者和食饵的动力学行为,分析了模型的有界性、平衡点的稳定性和类型以及Hopf分支现象。
Abstract: In this paper, the dynamic behavior of predator and prey is studied by adding the restriction of the Allee effect of predator to the Leslie Gower model in the case of constant predator capture. In this paper, the boundedness of the model, the stability and type of equilibrium point and Hopf branching phenomenon are analyzed.
文章引用:孙营, 孙福芹. 具有常数捕获猎物和Allee效应的Leslie-Gower模型的动力学研究[J]. 应用数学进展, 2024, 13(7): 3425-3432. https://doi.org/10.12677/aam.2024.137328

1. 引言

捕食–食饵模型是生态动力学中经典的数学模型[1],同时也是生态学和生物数学一个重要的研究课题。它描述了捕食者、食饵之间相互作用的动态关系,其中一个经典的模型之一是Leslie-Gower模型[2],当捕食者种群数量与其喜爱的食物种群数量成比例时,一般可假设捕食者的最大环境容纳量与食饵种群数量成正比。Leslie-Gower模型同时考虑捕食者和猎物两个种群的数量变化,其中捕食者数量变化受猎物数量的影响,猎物数量变化也受捕食者数量的影响,模型中猎物增长模型为线性方程,捕食者受环境容纳量限制。为了研究捕食者–食饵模型的动力学,其中从生物学的角度研究不同的分支是很重要的。捕食者–食饵模型的动力学变化可以通过模型参数的分支来实现。如Gupta等人研究了一个具有非线性捕食的Leslie-Gower型食饵–捕食者模型的稳定性和分支[3]

经典的Leslie-Gower模型为:

{ dx dt =rx( 1 x k )axy dy dt =ρy( 1 y bx ) (1.1)

其中xy分别表t时刻被捕食者和捕食者的种群密度,r表示被捕食者的内在增长率,ρ表示捕食者的内在增长率,k为被捕食者的环境容纳量,ax表示单位时间内被捕食者捕获率,bx为捕食者的承载能力, y/ ( bx ) 被称作Leslie-Gower项[4],其中上述参数都为正数。

2. 模型组成和初步结果

2.1. 模型介绍

本文在经典的Leslie-Gower模型被捕食者的Allee效应[5],则被捕食者在t时刻的增长模型为:

dx dt =rx( 1 x k )( xA )axy (1.2)

此外,本文还考虑被捕食者的商业价值,对被捕食者进行常数捕获的Holling-I型模型,则被捕食者在t时刻的增长模型为:

dx dt =rx( 1 x k )( xA )axyh (1.3)

其中0 < A < 1,且A为被捕食者中Allee效应的阈值,Allee效应曾在研究鱼缸试验中被提出,并且大量的研究[6]表明了当物种种群密度很低的时候会发生Allee效应,同时Allee效应对一些面临数量威胁的物种有着十分重要的研究和指导意义。

目前对于Leslie-Gower捕食者被捕食者模型的研究篇幅较多,例如2023年Yong Yao [7]分析了具有群体行为和不断捕获猎物的Leslie-Gower捕食者–食饵模型的动力学问题;为了研究不同捕捞方法物种之间的相互作用,May [8]提出了对捕食者和被捕食者都进行常数捕获的Leslie-Gower的Holling-I型的捕食者与被捕食者模型,还有一些学者研究了Leslie-Gower的Holling-II型的捕食者与被捕食者模型[9]

基于以上的一些研究,本文提出的改进的Leslie-Gower模型为:

{ dx dt = r 1 x( 1 x k )( xA )axyh dy dt = r 2 y( 1 y x ) (1.4)

对模型系统(1.4)利用无量纲化的方法减少参数,以达到化简模型系统的目的。首先作如下变换:

x ˜ = x k y ˜ = y k t ˜ = r 1 kt     a ˜ = a r 1 h ˜ = h 1 r 1 k 2 A ˜ = A k

为了方便我们将变换后变量、参数做如下记法:

x ˜ ~x y ˜ ~y t ˜ ~t     a ˜ ~a h ˜ ~h A ˜ ~A

我们使用符号“ ”表示“记为”。

则模型简化为:

{ dx dt =x( 1x )( xA )axyh dy dt =sy( 1 y x ) (1.5)

其中 0<A<1 ,其余参数均为正。

2.2. 有界性

定理1 系统(1.5)的解是有界的。

证明:当 x0 时,

dx dt x( 1x )( xA )

xB e 0 t ( 1s )( sA )ds

由Bellman-Gronwall不等式可知:

xB e t 3 3 + ( 1+A ) 2 t 2 At

因此 lim t supx( t )=0 ,令 Z( t )=x( t )+y( t ) ,对t两边求导得:

dZ(t) dt = dx dt + dy dt =x( 1x )( xA )axyh+sy( 1 y x ) =x[ ( 1x )( xA )ay h x ]+y( s sy x ) x( 1x )( xA )+y( s sy x ) x 2 ( 1x )+y( s sy x )

4 27 Z( t )+x+y+y( s sy y ) 4 27 Z( t )+0+y( s sy x +1 ) Z( t )+ 4 27 +( s y 2 +sy+y )

f( y )=s y 2 +sy+y ,则:

f max = ( s+1 ) 2 4s dZ( t ) dt z( t )+ 4 27 + ( s+1 ) 2 4s

再令 M= 4 27 + ( s+1 ) 2 4s ,可得 dZ( t ) dt z( t )+M

由Gronwall不等式得: Z( t )M+( z( 0 )M ) e t ,故 lim t z( t )=x( t )+y( t )=M

综上,系统(1.5)的解是有界的。

2.3. 平衡点分析

2.3.1. 平衡点的存在性

对于系统(1.5)在(0, 0)处无定义,故(0, 0)不是平衡点。

dx/ dt =0 dy/ dt =0 ,则考虑关于xy的方程组为:

{ x( 1x )( xA )axyh=0 sy( 1 y x )=0 (1.6)

定理2 系统(1.5)的平衡点为:

E 1 =( x 2 ,0 )  E 2 =( x 3 ,0 )  E 3 =( x 4 ,0 )  E 4 =( x 6 , x 6 )  E 5 =( x 7 , x 7 )  E 6 =( x 8 , x 8 )

其中

x 1 = A+1 A 2 A+1 3 x 2 = A+1+ A 2 A+1 3 x 5 = ( Aa+1 ) ( Aa+1 ) 2 3A 3 x 6 = ( Aa+1 )+ ( Aa+1 ) 2 3A 3

对于模型系统正平衡点的讨论在本节占主要部分,这是后续对主要结论验证的基础。

考虑当y=0时系统(1.5)正平衡点个数,相当于研究下面一元三次方程的正根的分布情况:对于方程 x 3 x 2 A x 2 +Ax+h=0 ,我们令 f( x )= x 3 ( A+1 ) x 2 +Ax+h ,下面分析 f( x ) 根的情况:

f ( x )=3 x 2 2( A+1 )x+A ,由二次函数相关知识可得 Δ>0

f ( x )=0 的解为: x 1 = A+1 A 2 A+1 3 x 2 = A+1+ A 2 A+1 3 ,分析导函数与原函数关系知:

f( x )= x 3 ( A+1 ) x 2 +Ax+h ( 0, x 1 ) ( x 2 , ) 上单调递增,在 ( x 1 , x 2 ) 上单调递减,下面讨论 f( x ) 根的三种情况:

1):当 f( x 2 )>0 时, f( x )=0 x>0 没有正解,故此种情况无平衡点;

2):当 f( x 2 )=0 时, f( x )=0 x>0 有一个正解,故此种情况有一个平衡点,设为 E 1 =( x 2 ,0 )

3):当 f( x 2 )<0 时, f( x )=0 x>0 有两个正解,设为 x 3 x 4 ,其中 0< x 1 < x 3 < x 2 < x 4 ,故此种情况有两个平衡点,设为 E 2 =( x 3 ,0 ) E 3 =( x 4 ,0 )

同理考虑当 x=y 时系统(1.5)正平衡点个数,系统转化为:

x 3 ( Aa+1 ) x 2 +Ax+h=0 ,令 g( x )= x 3 ( Aa+1 ) x 2 +Ax+h ,由二次函数相关知识得:

Δ=4 ( Aa+1 ) 2 12A=4[ ( Aa+1 ) 2 3A ]

由于 Δ<0 Δ=0 以及 Δ>0 x= 1 3 ( Aa+1 )0 g( x )=0 无正解,

故下面只讨论当 Δ>0 且对称轴 x= 1 3 ( Aa+1 )0 时解得分布情况,即 0<a<A+1 3A ,设 g ( x )=0 的解为 x 5 x 6 ,分析得知 x 5 x 6 同号且大于零, x 5 = ( Aa+1 ) ( Aa+1 ) 2 3A 3      x 6 = ( Aa+1 )+ ( Aa+1 ) 2 3A 3

分析导函数与原函数关系知:

x 3 ( Aa+1 ) x 2 +Ax+h=0 ( 0, x 5 ) ( x 6 , ) 上单调递增,在 ( x 5 , x 6 ) 上单调递减,下面讨论 g( x ) 根的三种情况:

1) 当 g( x 6 )>0 时, g( x )=0 无正解,故无平衡点;

2) 当 g( x 6 )=0 时, g( x )=0 有一个正解,故有一个平衡点,设为 E 4 =( x 6 , x 6 )

3) 当 g( x 6 )<0 时, g( x )=0 有两个正解,故有两个平衡点,设为 E 5 =( x 7 , x 7 ) E 6 =( x 8 , x 8 )

2.3.2. 平衡点的类型

下面我们考虑系统(1.5)在平衡点的类型,系统(1.5)的Jacobian矩阵为 J( E ) ,我们令

J( E )=( a 11 a 12 a 21 a 22 )=( 3 x 2 +2( 1+A )xayA ax s y 2 x 2 s 2sy x )

其中 a 21 >0 ,且

det( J )= λ 1 λ 2 =[ 3 x 2 +2( 1+A )xayA ]( s 2sy x )+ as y 2 x

Tr( J )= λ 1 + λ 2 =3 x 2 +2( 1+A )xA+s

1) 下面考虑当 y=0 时,Jacobian矩 J( E ) 列式和迹分别为:

det( J )= λ 1 λ 2 =[ 3 x 2 +2( 1+A )xA ]s

Tr( J )= λ 1 + λ 2 =3 x 2 +2( 1+A )xA+s

在平衡点 E 1 =( x 2 ,0 ) det( J( E 1 ) )=0 Tr( J( E 1 ) )=s>0 ,因此雅可比矩阵 J( E 1 ) 的两个特征根 λ 1 =0 λ 2 >0 λ 1 >0 λ 2 =0 ,故 E 1 =( x 2 ,0 ) 为高阶平衡点[10] [11]

同理可得:

det( J( E 2 ) )=[ 3 x 4 2 2( 1+A ) x 4 +A ]s<0

Tr( J( E 2 ) )=[ 3 x 3 2 2( A+1 )x+A ]+s>0

det( J( E 3 ) )=[ 3 x 4 2 2( 1+A ) x 4 +A ]s<0

因此雅可比矩阵 J( E 2 ) 的两个特征根 λ 1 >0 λ 2 >0 ,故 E 2 =( x 3 ,0 ) 为结点,且为不稳定结点。

雅可比矩阵 J( E 3 ) 的两个特征根 λ 1 >0 λ 2 <0 λ 1 <0 λ 2 >0 ,故 E 3 =( x 4 ,0 ) 为鞍点。

2) 考虑 y=x 时系统(1.5)平衡点类型,则此时Jacobian矩阵 J( E ) 的行列式和迹分别为:

det( J )= λ 1 λ 2 =s[ 3 x 2 2( Aa+1 )x+A ]

Tr( J )= λ 1 + λ 2 =3 x 2 +2( 1+A a 2 )xAs

在平衡点 E 4 =( x 6 , x 6 ) det( J( E 4 ) )=0 因此雅可比矩阵 J 有为零的特征根,

E 4 =( x 6 , x 6 ) 为高阶平衡点。

同理可得:

det( J( E 5 ) )=s[ 3 x 7 2 2( Aa+1 ) x 7 +A ]<0

det( J( E 6 ) )=s[ 3 x 8 2 2( Aa+1 ) x 8 +A ]>0

Tr( J( E 6 ) )=3 x 8 2 +2( A a 2 +1 ) x 8 As

因此雅可比矩阵 J( E 5 ) 的特征根 λ 1 >0 λ 2 <0 λ 1 <0 λ 2 >0 ,故 E 5 =( x 7 , x 7 ) 为鞍点。由于 Tr( J( E 6 ) ) 的正负性未知,故下面讨论 Tr( J( E 6 ) ) 的正负性:

0<a<A+1 3A x 8 > ( A a 2 +1+ 3( A+s ) ( A a 2 +1 ) 2 )/3 时,平衡点 E 6,1 =( x 8 , x 8 ) 的秩 Tr( J( E 6,1 ) )>0 ,因此雅可比矩阵 J( E 6,1 ) 的特征根 λ 1 <0 λ 2 <0 ,故为稳定结点。

0<a<A+1 3A

Aa+1+ ( Aa+1 ) 2 3A 3 < x 8 < A a 2 +1+ 3( A+s ) ( A a 2 +1 ) 2 3 时,

平衡点 E 6,2 =( x 8 , x 8 ) 的秩 Tr( J( E 6,2 ) )>0 ,因此雅可比矩阵 J( E 6,2 ) 的特征根 λ 1 >0 λ 2 >0 ,故 E 6,2 =( x 8 , x 8 ) 为不稳定结点。

3. Hopf分岔

定理3 0<a<A+1 3A 时,若

x 8 > A a 2 +1+ 3( A+s ) ( A a 2 +1 ) 2 3 .

则系统(1.5)在 E 6,1 =( x 8 , x 8 ) 处发生Hopf分岔。

证:模型系统(1.5)的Jacobian矩阵为:

J( E )=( a 11 a 12 a 21 a 22 )=( 3 x 2 +2( 1+A )xayA ax s y 2 x 2 s 2sy x ) det( J( E ) )= λ 1 λ 2 =[ 3 x 2 +2( 1+A )xayA ]( s 2sy x )+ as y 2 x

Tr( J( E ) )= λ 1 + λ 2 =3 x 2 +2( 1+A )xA+s

设特征方程为 P( λ )= λ 2 Tr( J( E ) )λ+det( J( E ) ) ,当 λ 1 + λ 2 =0 时,会出现一对虚根: λ 1 =i det( J ) λ 2 =i det( J ) ,我们选取s为分支参数,故当 y=x 时,

J( E )=( a 11 a 12 s s ) Tr( J )= a 11 s Det( J )=s( a 11 + a 12 )

则需满足 a 11 >0 ,则 P( λ )= λ 2 ( a 11 s )λs( a 11 + a 12 ) ,且 λ( s )=α( s )+iβ( s )

α( s )= a 11 s 2 β( s )= 4s( a 11 + a 12 ) ( a 11 s ) 2

α( a 11 )=0 Re( dλ ds )= 1 2 0 ,证毕。

下面计算第一Lyapunov系数 σ 来研究Hopf分支的稳定性[12]

作变换 X=x x 8 Y=y y 8 ,这样将平衡点移到原点处,然后在原点处对系统泰勒展开到三阶,则有:

dX dt = a 10 X+ a 01 Y+ a 11 XY+ a 20 X 2 + a 02 Y 2 + a 30 X 3 + a 03 Y 3 + a 21 X 2 Y+ a 12 X Y 2 +P( X,Y )

dY dt = b 10 X+ b 01 Y+ b 11 XY+ b 20 X 2 + b 02 Y 2 + b 30 X 3 + b 03 Y 3 + b 21 X 2 Y+ b 12 X Y 2 +Q( X,Y )

其中:

a 10 =3 x 8 2 +2( 1+A ) x 8 a y 8 A a 01 =a x 8 a 11 =a a 20 =3 x 8 +A+1 a 02 =0 a 03 =3

a 03 =0 a 21 =0 a 12 =0 b 10 = s y 8 2 x 8 2 b 01 =s 2s y 8 y 8 b 11 = 2s y 8 x 8 2 b 20 = 2s y 8 2 x 8 3 b 02 = 2s x 8

b 30 = 6s y 8 2 x 8 4 b 03 =0 b 21 = 4s y 8 x 8 3 b 21 = 4s y 8 x 8 3 b 12 = 2s x 8 2

P( X,Y )= i+j=4 + a ij x i y j Q( X,Y )= i+j=4 + b ij x i y j

则由文献[13]中公式知系统的第一Lyapunov系数 σ 为:

σ= 3π 2 a 01 Δ 3 2 { [ a 10 b 10 ( a 11 b 02 + a 02 b 11 + a 11 2 )+ a 10 a 01 ( a 20 b 11 + a 11 b 02 + b 11 2 ) + b 10 2 ( a 11 a 02 +2 a 02 b 02 ) 2 a 10 b 10 ( b 02 2 a 20 a 02 )2 a 10 a 01 ( a 20 2 b 20 b 02 ) a 01 2 ( 2 a 20 b 20 + b 11 b 20 )+( a 10 b 01 2 a 10 2 )( b 11 b 02 a 11 a 20 ) ( a 10 2 + a 01 b 10 )[ 3( b 10 b 03 a 01 a 30 )+2 a 10 ( a 21 + b 12 )+( b 10 a 21 a 01 b 21 ) ] }

其中 Δ= a 10 b 01 a 01 b 10 =s[ 3 x 8 2 +2( 1+A a 2 )xA ]+as x 8 =sTr( J )+as x 8

因为 Tr( J )<0 ,所以 Δ= a 10 b 01 a 01 b 10 >0 ,则由上述分析可得

σ>0 时: E 6,1 =( x 8 , x 8 ) 产生亚临界Hopf分支;

σ<0 时: E 6,1 =( x 8 , x 8 ) 产生超临界Hopf分支。

4. 结语

本文研究了具有不断收获猎物和Allee效应的Leslie-Gower捕食者–被捕食者模型,分析了该系统的有界性、平衡点稳定性和类型以及Hopf分支。通过Gronwall不等式分析系统的有界性,计算雅可比矩阵分析平衡点的稳定性和类型,以及计算第一Lyapunov系数分析产生亚临界分叉和超临界分叉的条件。

参考文献

[1] 周艳, 张存华. 具有集群行为的捕食者-食饵反应扩散系统的稳定性和Turing不稳定性[J]. 山东大学学报(理学版), 2021, 56(7): 73-81.
[2] 李梦婷, 周文. 具有群体性行为的捕食者-食饵系统稳定性分析[J]. 安庆师范大学学报: 自然科学版, 2023, 29(2): 103-109.
[3] Cai, L., Chen, G. and Xiao, D. (2012) Multiparametric Bifurcations of an Epidemiological Model with Strong Allee Effect. Journal of Mathematical Biology, 67, 185-215.
https://doi.org/10.1007/s00285-012-0546-5
[4] 李彦, 陈博. 一类扩散Leslie型捕食-食饵模型的动力学性质[J]. 淮阴师范学院学报(自然科学版), 2023, 22(4): 283-288.
[5] González-Olivares, E., Mena-Lorca, J., Rojas-Palma, A. and Flores, J.D. (2011) Dynamical Complexities in the Leslie-Gower Predator-Prey Model as Consequences of the Allee Effect on Prey. Applied Mathematical Modelling, 35, 366-381.
https://doi.org/10.1016/j.apm.2010.07.001
[6] 方侃, 曾怀杰, 陈晓英. 具有捕食者Allee效应的Leslie-Gower模型的动力学分析[J]. 延边大学学报(自然科学版), 2022, 48(1): 25-29, 40.
[7] Yao, Y. (2023) Dynamics of a Leslie-Gower Type Predator-Prey System with Herd Behavior and Constant Harvesting in Prey. arXiv: 2307.12117.
[8] May, R.M., Beddington, J.R., Clark, C.W., Holt, S.J. and Laws, R.M. (1979) Management of Multispecies Fisheries. Science, 205, 267-277.
https://doi.org/10.1126/science.205.4403.267
[9] Aziz-Alaoui, M.A. and Daher Okiye, M. (2003) Boundedness and Global Stability for a Predator-Prey Model with Modified Leslie-Gower and Holling-Type II Schemes. Applied Mathematics Letters, 16, 1069-1075.
https://doi.org/10.1016/s0893-9659(03)90096-6
[10] Brauer, F. (1979) Decay Rates for Solutions of a Class of Differential-Difference Equations. SIAM Journal on Mathematical Analysis, 10, 783-788.
https://doi.org/10.1137/0510074
[11] Harcourt, D.G. (1971) Population Dynamics of Leptinotarsa decemlineata (Say) in Eastern Ontario: III. Major Population Processes. The Canadian Entomologist, 103, 1049-1061.
https://doi.org/10.4039/ent1031049-7
[12] González-Olivares, E., Mena-Lorca, J., Rojas-Palma, A. and Flores, J.D. (2011) Dynamical Complexities in the Leslie-Gower Predator-Prey Model as Consequences of the Allee Effect on Prey. Applied Mathematical Modelling, 35, 366-381.
https://doi.org/10.1016/j.apm.2010.07.001
[13] Perko, L. (2001) Differential Equations and Dynamical Systems. Springer.