红松鼠种群保护的状态反馈脉冲模型研究

Study on State Feedback Impulsive Model of Red Squirrel Protection

摘要: 英国本土红松鼠因外来灰松鼠的入侵而面临灭绝危机。针对目前的情况，文章根据疱疹病毒对两种松鼠影响程度的差异，建立竞争模型，并依据种群数量施加状态反馈脉冲，进而建立带有状态反馈的脉冲模型。利用后继函数法和微分方程的有关理论，研究了系统阶一周期解的存在性和稳定性。所得结论为红松鼠的保护工作提供了一定的理论依据。数值分析结果验证了文中定理的正确性。

Abstract: British native red squirrels are facing extinction crisis due to the invasion of foreign grey squirrels. In view of the current situation, the paper established a competitive model according to the difference of the impact of herpes virus on the two kinds of squirrels, and applied state feedback pulse according to the population size, and then established a pulse model with state feedback. The existence and stability of first-order periodic solutions of the system are studied by using the theory of subsequent functions and differential equations. The conclusion provides a theoretical basis for the protection of red squirrel. The numerical results verify the correctness of the theorem.

1. 背景

19世纪初期，灰松鼠从美国被引到了英国，灰松鼠比英国土生土长的红松鼠强壮，且繁殖能力强，在抢占食物、栖息地的竞争中，红松鼠总是处于劣势 [1]。灰松鼠还携带疱疹病毒(squirrelpox)，这种病毒会引发一种滤过性皮肤病，对灰松鼠自身影响较小，对红松鼠影响较大。在英国，有60%的红松鼠都携带这种病毒，患病后的红松鼠健康状况急剧下降，无法觅食，最快在15天内死亡。2008年的统计数据表明，红松鼠的数量比2007年锐减了60%，这个数据大大高于此前每年10%左右的下降率。

对红松鼠的保护刻不容缓，利用数学模型研究红松鼠保护引起了学者的兴趣。刘琼在文献 [1] 中，考虑了捕杀灰松鼠，投放药物治疗红松鼠自身疾病的措施，建立了依赖红松鼠密度施加状态反馈脉冲控制模型，并对模型进行了相应的分析。因为状态反馈脉冲微分方程可以描述某些运动状态在固定或不固定时刻加速变化或跳跃，使模型更切合实际。近年来，学者用它解决了很多问题 [2] [3] [4] [5] [6]。Jiang G R研究了带有阶段结构的害虫防治模型的脉冲效应，得到了系统存在周期解的充分条件 [7]。Jiang G R研究具Holling II类功能反应和状态反馈脉冲控制的食饵–捕食模型，得到了非平凡正解的存在性和稳定性的充分条件 [8]。

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

2. 模型建立

由于疱疹病毒对红松鼠影响较大，对灰松鼠影响较小，建立如下模型：

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=x_1\left(k_1-m_1{x}_1-{\beta }_{12}{x}_2-r_1{x}_1{x}_2\right)\\ \frac{\text{d}{x}_2}{\text{d}t}=x_2\left(k_2-m_2{x}_2-{\beta }_{21}{x}_1-r_2{x}_1{x}_2\right)\end{array}\}{x}_1>x_1^{\ast }\\ \begin{array}{l}\Delta x_1=\alpha x_1\\ \Delta x_2=-\beta \left(x_2-x_2^{\ast }\right)\end{array}\}{x}_1=x_1^{\ast }\end{array}$ (2.1)

其中， $x_1$ 和 $x_2$ 分别表示红松鼠和灰松鼠在t时刻的密度， $k_1$ 和 $k_2$ 分别表示红松鼠和灰松鼠的增长率， $m_1$ 和 $m_2$ 分别表示红松鼠和灰松鼠的种内竞争系数， ${\beta }_{12}$ 和 ${\beta }_{21}$ 分别表示红松鼠和灰松鼠的种间竞争系数， $r_1$ 和 $r_2$ 分别表示红松鼠和灰松鼠感染病毒的速率。

3. 无脉冲模型分析

模型(2.1)对应的无脉冲模型为：

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=x_1\left(k_1-m_1{x}_1-{\beta }_{12}{x}_2-r_1{x}_1{x}_2\right)\\ \frac{\text{d}{x}_2}{\text{d}t}=x_2\left(k_2-m_2{x}_2-{\beta }_{21}{x}_1-r_2{x}_1{x}_2\right)\end{array}$ (3.1)

3.1. 系统的平衡点分析

系统(3.1)的平衡点为 $E_0\left(0,0\right)$， $E_1\left(0,\frac{k_2}{m_2}\right)$， $E_2\left(\frac{k_1}{m_1},0\right)$ 和 $E_3\left(x_1^{\ast },x_2^{\ast }\right)$，其中

$x_1^{\ast }=\frac{\rho \pm \sqrt{\left[{\rho }^2-4r_1^2\left({\beta }_{21}{r}_1-r_2{m}_1\right)\left(m_2{k}_1-{\beta }_{12}{k}_2\right)\right]}}{2r_1\left({\beta }_{21}{r}_1-m_1{r}_2\right)}$， $x_2^{\ast }=\frac{k_1-x_1^{\ast }}{{\beta }_{12}+r_1{x}_1^{*}}$，

$\rho =m_1{m}_2{r}_1-{\beta }_{21}{\beta }_{12}{r}_1$。当满足 $\frac{r_1}{r_2}<\frac{m_1}{{\beta }_{21}}$， $\frac{{\beta }_{12}}{m_2}<\frac{m_1}{{\beta }_{21}}$， $\frac{m_2}{{\beta }_{12}}>\frac{k_2}{k{}_1}$，存在一个正根 $E_3\left(x_1^{\ast },x_2^{\ast }\right)$。

容易知道 $E_0\left(0,0\right)$ 是不稳定的结点；当 $\frac{k_1}{k_2}>\frac{{\beta }_{12}}{m_2}$ 时， $E_1\left(0,\frac{k_2}{m_2}\right)$ 是鞍点，当 $\frac{k_1}{k_2}<\frac{{\beta }_{12}}{m_2}$ 时， $E_1\left(0,\frac{k_2}{m_2}\right)$ 是稳定的结点；当 $\frac{k_2}{k_1}>\frac{{\beta }_{21}}{m_1}$ 时， $E_2\left(\frac{k_1}{m_1},0\right)$ 是鞍点，当 $\frac{k_2}{k_1}<\frac{{\beta }_{21}}{m_1}$ 时， $E_2\left(\frac{k_1}{m_1},0\right)$ 是稳定的结点；当

$\left(m_1+r_1{x}_2^{\ast }\right)\left(m_2+r_1{x}_1^{\ast }\right)<\left({\beta }_{12}+r_1{x}_1^{\ast }\right)\left({\beta }_{21}+r_2{x}_2^{\ast }\right)$ 时， $E_3\left(x_1^{\ast },x_2^{\ast }\right)$ 是鞍点，当 $\left(m_1+r_1{x}_2^{\ast }\right)\left(m_2+r_1{x}_1^{\ast }\right)<\left({\beta }_{12}+r_1{x}_1^{\ast }\right)\left({\beta }_{21}+r_2{x}_2^{\ast }\right)$ 时， $E_3\left(x_1^{\ast },x_2^{\ast }\right)$ 是稳定的结点或焦点。

假设参数：

$k_1=0.05;k_2=0.03;m_1=0.02;m_2=0.3;{\beta }_{12}=0.008;{\beta }_{21}=0.006;r_1=0.00612;r_2=0.006$

如图1我们看出在 $E_3\left(x_1^{\ast },x_2^{\ast }\right)$ 点处是稳定的。

Figure 1. A track map of the neighborhood of E₃

图1. E₃附近的轨线图

3.2. 系统不存在极限环

我们用Bendixon-DuLac检验极限环的存在性，记

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=x_1\left(k_1-m_1{x}_1-{\beta }_{12}{x}_2-r_1{x}_1{x}_2\right)=f_1\left(x_1,x_2\right)\\ \frac{\text{d}{x}_2}{\text{d}t}=x_2\left(k_2-m_2{x}_2-{\beta }_{21}{x}_1-r_2{x}_1{x}_2\right)=g_1\left(x_1,x_2\right)\end{array}$ (3.2)

系统(3.2)中，函数 $f_1\left(x_1,x_2\right)$， $g_1\left(x_1,x_2\right)$ 在第一象限是光滑的，取 ${\rm B}\left(x_1,x_2\right)=\frac{1}{x_1{x}_2}$，则

$\frac{\partial \left({\rm B}{f}_1\right)}{\partial x_1}+\frac{\partial \left({\rm B}{g}_1\right)}{\partial x_2}=-\frac{m_1}{x_2}-r_1-\frac{m_2}{x_1}-r_2<0$。

因此，系统(3.2)在第一象限不存在极限环。

4. 脉冲模型分析

4.1. 阶一周期解的存在性

定理4.1系统(2.1)存在阶一周期解

证明：根据该平面系统上的轨线走势可知 [10]，我们假设A点坐标是 $A\left(x_{1A},x_{2A}\right)$， $A_1$ 的坐标是 $\left(x_{1A_1},x_{2A_1}\right)$， $A^{\ast }$ 的坐标是 $A^{\ast }\left(x_{1A^{\ast }},x_{2A^{\ast }}\right)$，假设B的坐标是 $B\left(x_{1B},x_{2B}\right)$， $B_1$ 的坐标是 $\left(x_{1B_1},x_{2B_1}\right)$， $B^{\ast }$ 的坐标是 $B^{\ast }\left(x_{1B^{\ast }},x_{2B^{\ast }}\right)$， $p_0$ 的坐标是 $p_0\left(x_{1p_0},x_{2p_0}\right)$，从 $A\in N$ 点出发的轨线与脉冲集M交于点 $A^{\ast }$，再由脉冲函数映射到 $A_1\in N$，此时 $A_1$ 在A的下方，必定有 $x_2^{\ast }<x_{2A_1}<x_{2A^{\ast }}<x_{2A}$，则 $F\left(A\right)=x_{2A_1}-x_{2A}<0$，从 $B\in N$ 点出发的轨线与脉冲集M交于点 $B^{\ast }$，再由脉冲函数映射到 $B_1\in N$， $B_1$ 在B的上边， $x_{2B}<x_{2B_1}<x_{2B^{\ast }}<x_2^{\ast }$，则 $F\left(B\right)=x_{2B_1}-x_{2B}>0$ 使得 $p_0\in N,p_1\in M,\phi :p_1\to p_2\in N$ 有 $F\left(p\right)=x_{2p_2}-x_{2p_1}=0$， $p_0$ 一定在AB两点之间 [11]。所以系统(2.1)必存在阶一周期解 $\Gamma \left(p_0,T\right)$，如图2所示：

Figure 2. Phase diagram

图2. 相图

4.2. 阶一周期解的稳定性

定理4.2 假设 $\Gamma$ 是某状态反馈脉冲控制系统的阶一周期解。如果对于任何 $\epsilon >0$，在相集上存在点p的 $\delta$ 的领域 $U\left(p,\delta \right)$， $\delta >0$ 对任意点 $p_1\in U\left(p,\delta \right)$ 和以为 $p_1$ 初始点的状态反馈脉冲系统的轨线 $f\left(p_1,t\right)$，存在T，当 $t>T$ 时有 $\rho \left(f\left(p_1,t\right),\Gamma \right)<\epsilon$，则称 $\Gamma$ 为是轨道稳定的 [12]。

我们知道对任意 $\sigma \left(0<\sigma <1\right)$ 系统(1.1)，有个唯一的周期解 $\Gamma$。我们在相集N上确定一点A (如图3)，ABC轨线 + 线段CA = 阶一周期解 [13]。我们把它定义为 $\Gamma$，然后，另一接近A点的点 $S_0$，我们得到一系列的点 $\left\{S_0,S_1,\cdots ,S_k,S_{k+1},\cdots \right\}$，这里的 $S_{i+1}$ 是 $S_i$ 的后继点， $i=0,1,\cdots ,k,\cdots$ 在相集N上建立这些点，假设A的坐标是0， $s_0,s_1,\cdots ,s_k,s_{k+1},\cdots$ 是 $S_0,S_1,\cdots ,S_k,S_{k+1},\cdots$ 的坐标。我们有：

$s_k=\{\begin{array}{l}-d_k,S_k在点A的左侧\\ d_k,S_k在点A的右侧\end{array}$

这里的 $d_k$ 是 $S_k$ 和A之间的距离。

Figure 3. $s_0,s_1,\cdots ,s_k,s_{k+1},\cdots$ are the subsequent points of $S_0,S_1,\cdots ,S_k,S_{k+1},\cdots$ respectively

图3. $s_0,s_1,\cdots ,s_k,s_{k+1},\cdots$ 是 $S_0,S_1,\cdots ,S_k,S_{k+1},\cdots$ 的坐标

引理4.2 f来表示抽象脉冲微分动力系统(2.1)的映射 [14] [15]。当函数

${\left|\frac{\text{d}\bar{s}}{\text{d}s}\right|}_{s=0}<1$

是稳定的，否则是不稳定的。

一般来说，我们假设 $P\left(x,y\right)$ 和 $Q\left(x,y\right)$ 是系统(2.1)中的函数，具有任意阶偏导，我们假设在一定的轨道上，从A到C时间是T，周期闭合轨道 $\Gamma$ 也是T。

我们建立一个新的坐标，然而第二个点n作为法向量的个数(这里我们假设是正方向上的) (如图4所示)。由第一个坐标，我们可以建立曲线ABC的参数方程 [16] [17] [18] [19]。

Figure 4. Establish coordinate system $\left(s,n\right)$ on point A

图4. 存在点系统 $\left(s,n\right)$ 在点A上

$x=\varphi \left(s\right),y=\psi \left(s\right)$，这里的s是参数。我们表示坐标 $\left(x,y\right)$ 与 $\left(s,n\right)$ 之间的直线表达形式：

$x=\varphi \left(s\right)-n{\psi }^{\prime }\left(s\right),y=\psi \left(s\right)+n{\varphi }^{\prime }(\; s\; )$

根据(2.1)和(4.2)我们有：

$\frac{\text{d}n}{\text{d}s}=\frac{Q{\varphi }^{\prime }-P{\psi }^{\prime }-n\left(P{\varphi }^{″}+Q{\psi }^{″}\right)}{P{\varphi }^{\prime }+Q{\psi }^{\prime }}\triangleq F\left(s,n\right)$ (4.3)

很明显 $n=0$ 是上述方程(4.3)的解，一阶偏导数 $F\left(s,n\right)$。关于n是连续的。然后上述方程可以简化为：

$\frac{\text{d}n}{\text{d}s}={{F^{\prime }}_n\left(s,n\right)|}_{n=0}+{\rm O}(\; n\; )$

我们可以计算出：

${{F^{\prime }}_n\left(s,n\right)|}_{n=0}=\frac{P_0^2{Q}_{y_0}-P_0{Q}_0\left(P_{y_0}+Q_{x_0}\right)+Q_0^2{P}_{x_0}}{{\left(P_0^2+Q_0^2\right)}^{\frac{3}{2}}}\triangleq H\left(s\right)$ (4.4)

这里的 $P_0$ 和 $Q_0$ 代表的是P和Q在A的坐标。 $P_{x_0}$， $P_{y_0}$ 以及 $Q_{y_0}$， $Q_{y_0}$ 是 $n=0$ 时P和Q的偏导数，关于 $x,y$，很明显 $P_0{\varphi }^{″}+Q_0{\psi }^{″}=0$ 然后(4.4)的直线方程是：

$\frac{\text{d}n}{\text{d}s}=H\left(s\right)n$ (4.5)

解是： $n=n_0{\text{e}}^{{\displaystyle {\int }_0^{s}H\left(s^{\prime }\right)\text{d}{s}^{\prime }}},n_0=n(\; 0\; )$

定理4.1 假设h是轨线 $\stackrel{⌢}{ABC}$ 的长度，是系统(2.1)阶一周期解 $\Gamma$ 的一部分，如果把 $\stackrel{⌢}{ABC}$ 集成到一起，用 $H\left(s\right)$ 来表示是负的，说明 $\Gamma$ 是稳定的，否则是不稳定的。

定理4.2 假设该区域的阶一周期解 $\Gamma$ 是凸闭的。如果满足 ${\displaystyle {\int }_0^{{\rm T}}{P}_{x_0}+Q_{y_0}}<0$，说明 $\Gamma$ 是稳定的，否则是

不稳定的。

然后我们直接得出推论：

推论4.1 确定闭区域阶一周期解 $\Gamma$ 是凸的，如果 $P_{x_0}+Q_{y_0}<0$ 成立，证明 $\Gamma$ 是稳定的，否则是不稳定的。