红松鼠种群保护状态反馈脉冲控制第二类模型研究
Study on the Second Model of State Feedback Impulsive Model of Red Squirrel Protection
摘要: 英国本土红松鼠因外来灰松鼠的入侵而面临灭绝危机,为了保护红松鼠的生存量,英国政府提出既保护了英国本土红松鼠达到人们所希望的数量而灰松鼠也能保留人们所希望的种群规模,针对目前的情况,文章根据疱疹病毒对两种松鼠影响程度的差异,建立非Lota-Volterra竞争模型,详细分析了无脉冲时的正平衡点的情况,建立带有状态反馈的脉冲模型。利用后继函数法和微分方程的有关理论,研究了系统阶一周期解的存在性和稳定性。所得结论为红松鼠的保护工作提供了一定的理论依据。
Abstract: Native red squirrels in the UK are facing the danger of extinction due to the invasion of alien grey squirrels. In order to protect the stock of red squirrels, the British government proposed that the number of native red squirrels in the UK could be protected to the desired number while the population size of grey squirrels could be retained. In view of the present situation, according to the difference of herpes virus influence degree on two squirrels, the non-Lota-Volterra competition model is established, and the positive equilibrium point is analyzed in detail under the condition of no pulse, and the pulse model with state feedback is established. The existence and stability of first-order periodic solutions of the system are studied by using the method of successive functions and the theory of differential equations. The conclusions provided a theoretical basis for the conservation of red squirrels.
文章引用:王烁烁. 红松鼠种群保护状态反馈脉冲控制第二类模型研究[J]. 应用数学进展, 2021, 10(11): 3640-3649. https://doi.org/10.12677/AAM.2021.1011385

1. 背景

对红松鼠的保护刻不容缓 [1],利用数学模型研究红松鼠保护引起了学者的兴趣。近年来,灰松鼠种群增长比较快,而红松鼠受灰松鼠身上的疱疹病毒传染,红松鼠死亡率增高,而病毒对灰松鼠伤害较小 [2],因此灰松鼠种群增长比较快。英国本土红松鼠因外来灰松鼠的入侵而面临灭绝危机 [3],为了保护红松鼠的生存量,英国政府提出捕杀灰松鼠的方法限制灰松鼠的增长,这样就产生了一个问题:我们每年要捕杀多少灰松鼠,才能既保护了英国本土红松鼠达到人们所希望的数量而灰松鼠也能保留人们所希望的种群规模 [4]。

本文考虑疱疹病毒建立状态脉冲模型 [5]。

2. 模型建立

由于疱疹病毒对红松鼠影响较大,对灰松鼠影响较小,我们首先建立非Lota-Volterra竞争模型 [6]:

{ d x d t = x ( 1 x K 1 α y r y ) d y d t = y ( 1 y K 2 β x ) , y y (1)

{ Δ x = 0 Δ y = δ y , y = y

其中,x和y分别表示红松鼠和灰松鼠在t时刻的密度,系统(2)中微分方程组即为系统(1)两种群相互作用模型, k 1 k 2 分别表示红松鼠和灰松鼠的环境容纳量, α 表示每个灰松鼠占用红松鼠环境容纳量的数量, β 表示每个红松鼠占用灰松鼠环境容纳量的数量,r表示每个灰松鼠身上的病毒使红松鼠种群产生的死亡率。 y 表示允许灰松鼠保存的密度,也就是说在这个密度以内不会危及红松鼠的生存, δ 表示灰松鼠生存的密度达到 y 时,人们开始捕杀灰松鼠时的捕杀率,系统(1)中的: Δ x = x ( t + ) x ( t ) Δ y = y ( t + ) y ( t ) Δ x = 0 表示在捕杀灰松鼠时不要伤害红松鼠 [7]。

3. 无脉冲模型分析

模型(1)对应的无脉冲模型为:

{ d x d t = x ( 1 x K 1 α y r y ) d y d t = y ( 1 y K 2 β x ) (2)

其中,x和y分别表示红松鼠和灰松鼠在t时刻的密度, k 1 k 2 分别表示红松鼠和灰松鼠的环境容纳量, α 表示每个灰松鼠占用红松鼠环境容纳量的数量, β 表示每个红松鼠占用灰松鼠环境容纳量的数量,r表示每个灰松鼠身上的病毒使红松鼠种群产生的死亡率。 α y 表示灰松鼠t时刻占用红松鼠环境容纳量的数量。

3.1. 模型的向量场与一致有界性

定理1:系统(2)所有界最终一致有界。

证明:根据 d x d t = 0 d y d t = 0 得到水平等倾线为: y = K 2 β x ,为直线。

得到垂直等倾线为: x = K 1 α y 2 + ( r K 1 α ) y + K 1 ,为抛物线。

设直线 x = k 2 β 与等倾线交于一点A,G与x轴交于 A 1 ,过A作直线L,设L与y轴交于点B,直线L

记为 [8]:见图1

L = a x + y c = 0 , a > 0 , c > 0

直线段 A A 1 A B B O O A 1 ,围成一个闭区域G,见图2

d L d t = a d x d t + d y d t = a x [ 1 x K 1 α y ] + y [ 1 y K 2 β x ]

d L d t | L = 0 < 0

d x d t = x [ 1 x K 1 α y r y ] > 0

d y d t = y [ 1 y K 2 β x ] < 0

这说明轨线与线段L相交 [9],必进入区域G。

因此,我们有:系统(2)所有界最终一致有界。

3.2. 系统的平衡点分析

系统(2)存在四个平衡点: E 0 ( 0 , 0 ) E 1 ( 0 , K 2 ) E 2 ( K 1 , 0 ) E ( x , y )

系统(2)对应的雅各比矩阵:

J = ( 1 2 x K 1 α y r y α x 2 ( K 1 α y ) 2 r x β y 2 ( K 2 β x ) 2 1 2 x K 2 β x )

Figure 1. Vector filed

图1. 向量场

Figure 2. Vector filed

图2. 向量场

由于 { 1 x K 1 α y r y = 0 1 y K 2 β y = 0

得到方程:

α β 2 x 2 + ( α β + r β K 1 2 K 1 β α r 1 ) + α r K 2 2 α K 2 r K 1 K 2 + K 1 = 0

其中:

a = α r β 2

b = α β + r β K 1 2 K 1 β α r 1

c = α r K 2 2 α K 2 r K 1 K 2 + K 1 = ( α K 2 K 1 ) ( r K 2 1 )

Δ = ( α β + r β K 1 2 K 1 β α r 1 ) 2 4 α r β 2 ( α K 2 K 1 ) ( r K 2 1 )

1) 当 x 1 + x 2 = b a > 0 x 1 x 2 = c a > 0 Δ = 0

M1 α β + r β K 1 2 K 1 β α r 1 < 0 α β + r β K 1 < 2 K 1 β α r + 1

M2 ( α K 2 K 1 ) ( r K 2 1 ) > 0 α K 2 K 1 > 0 r K 2 1 > 0

M3 ( α K 2 K 1 ) ( r K 2 1 ) > 0 α K 2 K 1 < 0 r K 2 1 < 0

M4 ( α β + r β K 1 2 K 1 β α r 1 ) 2 = 4 α r β 2 ( α K 2 K 1 ) ( r K 2 1 )

系统(2)有一个正平衡点 E ( x , y ) y = K 2 β x * x * = 2 K 1 β α r r β K 1 α β + 1 2 α r β 2

此时系统的稳定性:

结论一: E 0 ( 0 , 0 ) 为不稳定的结点。

证明: E 0 ( 0 , 0 ) J ( 0 , 0 ) = ( 1 0 0 1 )

D = 1 , T = 2 , Δ = T 2 4 D = 0

Δ = T 2 4 D = 0 可知说明平衡点是结点。

所以 E 0 是不稳定的结点。

结论二: E 1 ( 0 , K 2 ) 为鞍点。

证明:平衡点, J ( 0 , K 2 ) = ( 1 r K 2 0 β 1 )

考虑M3成立

D = r K 2 1 < 0 , T = r K 2 , Δ = T 2 4 D = ( r K 2 2 ) 2

所以 E 2 是鞍点。

结论三: E 2 ( K 1 , 0 ) 为鞍点。

证明:平衡点 E 2 ( K 1 , 0 ) J ( K 1 , 0 ) = ( 1 α r K 1 0 1 2 K 1 K 2 β K 1 )

D = 1 < 0 , T = 0 , Δ = T 2 4 D = 4

所以 E 2 是鞍点。

结论四: E ( x , y ) 为稳定的结点。

证明:平衡点 E ( x , y ) J ( x , y ) = ( r y 1 α ( 1 r y ) 2 r x β 1 )

D = β ( α α r y + r x ) T = r y 2

Δ = T 2 4 D = ( r y 2 ) 2 4 β ( α α r y + r x )

考虑 M 1 , M 2 , M 3 , M 4 成立

D = 1 r y α β ( 1 r y ) 2 r β x = 1 r ( K 2 β x ) α β ( 1 r ( K 2 β x ) ) 2 r β x = 1 r ( K 2 β x ) α β ( 1 r K 2 + r β x ) 2 r β x = 1 r K 2 + r β x α β [ ( 1 r K 2 ) 2 + r 2 β 2 x 2 + 2 ( 1 r K 2 ) r β x ] r β x = 1 r K 2 + r β x α β ( 1 r K 2 ) 2 α r 2 β 3 x 2 2 α r β 2 ( 1 r K 2 ) x r β x

= 1 r K 2 α β ( 1 + r 2 K 2 2 2 r K 2 ) α r 2 β 3 x 2 2 α r β 2 ( 1 r K 2 ) x = 1 r K 2 α β α β r 2 K 2 2 + 2 α β r K 2 α r 2 β 3 x 2 2 α r β 2 ( 1 r K 2 ) x = 1 r K 2 + α β r K 2 α β r 2 K 2 2 + α β r K 2 α β α r 2 β 3 x 2 2 α r β 2 ( 1 r K 2 ) x = ( 1 r K 2 ) + α β r K 2 ( 1 r K 2 ) α β ( 1 r K 2 ) 2 α r β 2 x ( 1 r K 2 ) α r 2 β 3 x 2 = ( 1 r K 2 ) ( 1 + α β r K 2 α β 2 α r β 2 x ) α r 2 β 3 x 2

= ( 1 r K 2 ) ( 1 α β + α β r K 2 2 α r β 2 x ) α r 2 β 3 x 2 = ( 1 r K 2 ) ( 1 α β ( 1 r K 2 ) 2 α r β 2 x ) α r 2 β 3 x 2 = ( 1 r K 2 ) α β ( 1 r K 2 ) 2 2 α r β 2 x ( 1 r K 2 ) α r 2 β 3 x 2 = ( 1 r K 2 ) α β ( 1 r K 2 ) 2 2 α r β 2 2 K 1 β α r r β K 1 α β + 1 2 α r β 2 ( 1 r K 2 ) α r 2 β 3 ( 2 K 1 β α r r β K 1 α β + 1 2 α r β 2 ) 2 = ( 1 r K 2 ) α β ( 1 r K 2 ) 2 ( 2 K 1 β α r r β K 1 α β + 1 ) ( 1 r K 2 ) α r 2 β 3 ( 2 K 1 β α r r β K 1 α β + 1 2 α r β 2 ) 2

= ( 1 r K 2 ) α β ( 1 r K 2 ) 2 ( 2 K 1 β α r r β K 1 α β + 1 ) ( 1 r K 2 ) 1 4 α β ( 2 K 1 β α r r β K 1 α β + 1 ) 2 = ( 1 r K 2 ) [ 1 α β ( 1 r K 2 ) ( 2 K 1 β α r r β K 1 α β + 1 ) ] 1 4 α β ( 2 K 1 β α r r β K 1 α β + 1 ) 2 = ( 1 r K 2 ) ( 1 α β + α β r K 2 2 K 1 β α r + r β K 1 + α β 1 ) 1 4 α β ( 2 K 1 β α r r β K 1 α β + 1 ) 2 = ( 1 r K 2 ) ( α β r K 2 2 K 1 β α r + r β K 1 ) 1 4 α β ( 2 K 1 β α r r β K 1 α β + 1 ) 2 = ( 1 r K 2 ) ( α β r K 2 2 K 1 β α r + r β K 1 ) 4 α r β 2 ( α K 2 K 1 ) ( r K 1 1 ) 4 α β

= ( 1 r K 2 ) ( α β r K 2 2 K 1 β α r + r β K 1 ) r β ( α K 2 K 1 ) ( r K 1 1 ) = ( 1 r K 2 ) ( α β r K 2 2 K 1 β α r + r β K 1 ) + r β ( α K 2 K 1 ) ( 1 r K 1 ) = ( 1 r K 2 ) [ ( α β r K 2 2 K 1 β α r + r β K 1 ) + r β α K 2 r β K 1 ] = ( 1 r K 2 ) ( α β r K 2 2 K 1 β α r + r β α K 2 )

= ( 1 r K 2 ) ( 2 α β r K 2 2 K 1 β α r ) = 2 α β r ( 1 r K 2 ) ( K 2 K 1 ) = 2 α β r ( r K 2 1 ) ( K 1 K 2 ) > 0

所以 E ( x , y ) 是稳定的结点。

3.3. 系统极限环不存在性

定理2:系统(1)在第一象限不存在极限环。

证明:我们用DuLac极限环的存在性定理,记 B ( x , y ) = 1 x y

{ d x d t = x ( 1 x K 1 α y r y ) = f ( x , y ) d y d t = y ( 1 y K 2 β x ) = g ( x , y ) (2)

f 1 ( x , y ) = B ( x , y ) f ( x , y ) = 1 y x y ( K 1 α y ) r

g 1 ( x , y ) = B ( x , y ) g ( x , y ) = 1 x y x ( K 1 α y )

( B f ) x + ( B g ) y = 1 y ( K 1 α y ) 1 x ( K 2 β x ) < 0

因此,由DuLac定理得到结论,系统(1)在第一象限不存在极限环 [10]。

定理3:平衡点 E ( x , y ) 为全局渐近稳定。

证明:由结论4可知 E ( x , y ) 为稳定的焦点,而且 E ( x , y ) 是区域G中唯一的平衡点。又由定理2:系统(2)在第一象限不存在极限环。再由定理1:系统(1)所有界最终一致有界 [11]。

由此,我们就可以得出结论:平衡点 E ( x , y ) 为全局渐进稳定 [12]。

3.4. 全局稳定性

通过构造V函数证明无控制系统的全局稳定性 [13]。我们构造V函数:

V ( x , y ) = [ ( x x ) x ln ( x x ) ] + c [ ( y y ) y ln ( y y ) ]

d V d t = x x x d x d t + c y y y d y d t = ( x x x ) ( x x 2 K 1 α y r x y ) + c ( y y y ) ( y y 2 K 2 β x ) = ( x x ) ( 1 x K 1 β y r y ) + c ( y y ) [ 1 y K 2 β x ]

考虑到平衡点方程:

{ 1 x K 1 α y r y = 0 1 y K 2 β x = 0

d V d t = ( x x ) [ 1 x K 1 α y r y 1 + x K 1 α y + r y ] + c ( y y ) [ 1 y K 2 β x 1 + y K 2 β x ] = ( x x ) [ x K 1 α y r y + x K 1 α y + r y ] + c ( y y ) [ y K 2 β x + y K 2 β x ]

= ( x x ) [ r ( y y * ) + K 1 ( x x ) + x α ( y y ) ( K 1 α y ) ( K 1 α y ) ] + c ( y y ) [ K 2 ( y y ) + β ( x y x y ) ( K 2 β x ) ( K 2 β x ) ] = r ( y y * ) ( x x ) + ( x x ) [ K 1 ( x x ) + x α ( y y ) ( K 1 α y ) ( K 1 α y ) ] + c ( y y ) [ K 2 ( y y ) + β ( x y x y ) ( K 2 β x ) ( K 2 β x ) ]

= r ( y y * ) ( x x ) + K 1 ( x x ) ( x x ) ( K 1 α y ) ( K 1 α y ) + x α ( x x ) ( y y ) ( K 1 α y ) ( K 1 α y ) + [ c K 2 ( y y ) ( y y ) + c β ( y y ) ( x y x y ) ( K 2 β x ) ( K 2 β x ) ] = r ( y y * ) ( x x ) K 1 ( x x ) 2 ( K 1 α y ) ( K 1 α y ) x α ( x x ) ( y y ) ( K 1 α y ) ( K 1 α y ) + c K 2 ( y y ) ( y y ) ( K 2 β x ) ( K 2 β x ) + c β ( y y ) ( x y x y ) ( K 2 β x ) ( K 2 β x )

= r ( y y * ) ( x x ) K 1 ( x x ) 2 ( K 1 α y ) ( K 1 α y ) x α ( x x ) ( y y ) ( K 1 α y ) ( K 1 α y ) c K 2 ( y y ) 2 ( K 2 β x ) ( K 2 β x ) + c β ( y y ) ( x y x y ) ( K 2 β x ) ( K 2 β x ) = r ( y y * ) ( x x ) + c β ( y y * ) ( x y x y * ) ( K 2 β x ) ( K 2 β x ) K 1 ( x x ) 2 ( K 1 α y ) ( K 1 α y ) x α ( x x ) ( y y ) ( K 1 α y ) ( K 1 α y ) c K 2 ( y y ) 2 ( K 2 β x ) ( K 2 β x )

c = r ( x x ) ( K 2 β x ) ( K 2 β x ) β ( x y x y ) 时,

d V d t = [ K 1 ( x x ) 2 ( K 1 α y ) ( K 1 α y ) + x α ( x x ) ( y y ) ( K 1 α y ) ( K 1 α y ) + c K 2 ( y y ) 2 ( K 2 β x ) ( K 2 β x ) ] < 0

因为 d V d t 在第一象限内是定负的,所以在 E ( x , y ) 平衡点处是全局渐近稳定的 [14]。

4. 脉冲模型分析

4.1. 阶一周期解的存在性

定理4:系统(1)存在阶一周期解。

Figure 3. Successor function

图3. 后继函数

证明:脉冲集M: y = y ,相集N: y = ( 1 δ ) y ,直线 x = x 与相集交于一点A,过A点的轨线与脉冲集M交于一点 A 1 ,由向量场可以看出轨线 A A 1 的形状如图3所示,脉冲到相集上一点 A 0 ,这说明 A 0 为A的后继点 [15],同样可知从原点出来的一条轨线 B B 1 ,有 B 0 为B的后继点,在相集上建立坐标系,以N与Y轴的交点为坐标原点,分别设点A的坐标为a,点B的坐标为b,分别设点 A 0 的坐标为 a 0 B 0 点的坐标为 b 0 ,后继函数 [16]:

f ( A ) = a 0 a < 0 , f ( B ) = b 0 b > 0

这样,由文 [17] 的后继函数判别法证明了定理4,即脉冲状态反馈控制系统(1)存在阶一周期解 C C 1 。如图4

Figure 4. First order periodic solution

图4. 阶一周期解

4.2. 阶一周期解的稳定性

定理5:系统(1)的阶一周期解是稳定的。

证明:由定理2的证明知系统(2)中DuLac函数 B ( x , y ) 变换后得系统(3)。

{ d x d t = x ( 1 x K 1 α y r y ) = B ( x , y ) f ( x , y ) = f 1 ( x , y ) d y d t = y ( 1 y K 2 β x ) = B ( x , y ) g ( x , y ) = g 1 ( x , y ) (3)

由定理2的证明中的计算有:

( B f ) x + ( B g ) y = 1 y ( K 1 α y ) 1 x ( K 2 β x ) < 0

又系统(1)的轨线与系统(3)的轨线是完全相同的,由文献 [17] 即证明了阶一周期解 C C 1 是稳定的。

参考文献

[1] 陈兰荪. 数学生态学模型与研究方法(第二版), 生物数学丛书, 第19卷[M]. 北京: 科学出版社, 2017.
[2] 马知恩, 周义仓, 李承治. 常微分方程定性与稳定性方法[M]. 第2版. 北京: 科学出版社, 2015.
[3] Liu, Q., Zhang, M. and Chen, L.S. (2019) State Feedback Impulsive Therapy to SIS Model of Animal Infectious Disease. Physica A, 516, 222-232.
https://doi.org/10.1016/j.physa.2018.09.161
[4] Laksbmikantham, V., Bainov, D. and Simeonov, P. (1989) Theory of Impulsive Differential Equations. World Scientific, Singapore.
https://doi.org/10.1142/0906
[5] Wang, L.M., Liu, Z.J., et al. (2007) Impulsive Diffusion in Single Species Model. Chaos, Solitons & Fractals, 33, 1213-1219.
https://doi.org/10.1016/j.chaos.2006.01.102
[6] Liu, K.Y. and Chen, L.S. (2007) Harvesting Control for a Stage-Structured Predator-Prey Model with Ivlev’s Functional Response and Impulsive Stocking on Prey. Discrete Dynamics in Nature and Society, 2007, Article ID: 086482.
[7] 陈兰荪. 害虫治理与半连续动力系统几何理论[J]. 北华大学学报(自然科学版), 2011, 12(1): 1-9.
[8] 陈兰荪. “半连续动力系统”理论及应用[J]. 玉林师范学院学报, 2013, 34(2): 2-10.
[9] Wang, T.Y. and Chen, L.S. (2011) Nonlinear Analysis of a Microbial Pesticide Model with Impulsive State Feedback Control. Nonlinear Dynamics, 65, 1-10.
https://doi.org/10.1007/s11071-010-9828-x
[10] 陈兰荪. 半连续动力系统理论及应用[J]. 玉林师范学院学报, 2013, 34(2): 1-9.
[11] 桂占吉, 王凯华, 陈兰荪. 病虫害防治的数学理论与计算[M]. 北京: 科学出版社, 2014.
[12] 叶彦谦, 等. 极限环论[M]. 上海: 上海科学技术出版社, 1984.
[13] Zhang, M., Chen, L.S. and Li, Z.Y. (2019) Homoclinic Bifurcation of a State Feedback Impulsive Controlled Prey-Predator System with Holling-II Functional Response. Nonlinear Dynamics, 98, 929-942.
https://doi.org/10.1007/s11071-019-05235-8
[14] Fu, J.-B. and Chen, L.-S. (2018) Modelling and Qualitative Analysis of Water Hyacinth Ecological System with Two State-Dependent Impulse Controls. Complexity, 2018, Article ID: 4543976.
https://doi.org/10.1155/2018/4543976
[15] Liu, Q., Zhang, M. and Chen, L.S. (2019) State Feedback Impulsive Therapy to SIS Model of Animal Infectious Diseases. Physica A: Statistical Mechanics and its Applications, 516, 222-232.
https://doi.org/10.1016/j.physa.2018.09.161
[16] Zhang, M., Liu, K.Y., Chen, L.S. and Li, Z.Y. (2018) State Feedback Impulsive Control of Computer Worm and Virus with Saturated Incidence. Mathematical Biosciences & Engineering, 15, 1465-1478.
https://doi.org/10.3934/mbe.2018067
[17] Chen, L.S., Liang, X.Y. and Pei, Y.Z. (2018) The Periodic Solution of the Impulsive State Feedback Dynamical System. Communication in Mathematical Biology and Neuroscience, 2018, Article No. 14.