疫苗免疫与饱和发生率时滞系统的Hopf分支

Hopf Bifurcation in a Vaccine Immunity and Saturation Incidence Time Delay System

摘要

为探讨疫苗接种后病毒对人类影响的动态变化，构建含有非线性饱和发生率与疫苗接种率的SVIR传染病模型，并引入时滞\tau表征疫苗免疫作用下病毒的潜伏期。以时滞为变量，运用特征方程和分析方法，深入探讨模型在正平衡点处的稳定性、Hopf分支的产生条件及周期解特性，当时滞小于临界值{\tau }_{\mathrm{0}}时，系统稳定；超过临界值{\tau }_{\mathrm{0}}时，特征根穿越虚轴引发周期性震荡。通过MATLAB数值模拟设定疫苗保护率、潜伏期时长等参数，绘制不同时滞下各仓室人群数量的动态曲线，直观呈现系统从稳定态向周期振荡的转变。结果表明，时滞超临界值时模型输出呈显著周期性波动，与理论推导结论一致，验证了Hopf分支条件的正确性。

Abstract

To explore the dynamic changes in the impact of viruses on humans after vaccination， this paper constructed an SVIR infectious disease model that incorporated nonlinear saturation incidence and vaccination rates， introducing time delays τ to represent the incubation period of the virus under vaccine immunity. Using time delay as a variable， the paper applied characteristic equations and analytical methods to investigate the stability of the model at the positive equilibrium point， the conditions for Hopf bifurcation， and the characteristics of periodic solutions. The system is stable when the time delay is less than {\tau }_{\mathrm{0}}； when it exceeds {\tau }_{\mathrm{0}}， the characteristic roots cross the imaginary axis， leading to periodic oscillations. Through MATLAB numerical simulations， parameters such as vaccine efficacy and incubation period duration were set， and dynamic curves of the population quantities in different compartments were plotted under various time delays， visually depicting the transition of the system from a steady state to periodic oscillation. The results show that when the time delay exceeds the critical value， the model exhibits significant periodic fluctuations， consistent with the theoretical derivation and validating the correctness of the Hopf bifurcation conditions.

Graphical abstract

关键词

时滞 / 饱和发生率 / Hopf分支 / 传染病模型 / 稳定性

Key words

time delay / saturation incidence / Hopf bifurcation / infectious disease model / stability

在当今科技迅猛发展的时代，尽管人类在医学技术革新、公共卫生体系建设等领域成果斐然，但传染病对生命健康的威胁从未消退。以新冠肺炎疫情为例，世界卫生组织统计数据显示，全球累计确诊病例超6亿例，大量人口面临健康危机，经济秩序遭受冲击，充分暴露了传染病对社会运行的深远影响。制定传染病防控策略高度依赖理论模型支撑，模型精确度直接影响疾病控制效果。在流感传播模型中，通过模拟病毒在人群中的传播轨迹、不同场景下的感染概率等参数，可精准预测传播峰值，辅助公共卫生部门提前布局医疗资源、规划疫苗接种策略。学者对传染病模型的研究自古有之且从未停歇，从早期刻画易感者、感染者、康复者关系的SIR^［1］经典模型，到如今融合时滞效应、社交网络传播特性的复杂动力学模型，每一次理论突破都源于对传染病传播规律的深度探索。这些研究不仅推动流行病学理论迭代，更在埃博拉、登革热等疫情防控中发挥了关键作用。20世纪初，Kermack和Mckendrick^［2］开创性地提出了SIR仓室模型，并随着对现实情况理解的加深，不断融入各种影响因素，使得对传染病的研究更加精细。Sun等^［3］构建了一个考虑区域间旅行和历史疾病监测影响的两区域SIR模型，分析旅行率和监测数对病毒传播的影响。Das等^［4］研究带有饱和发生率的结核病SEIR模型，验证了在地方性平衡点的稳定性及周期解。Gao等^［5］研究了一个在时空异质环境中具有饱和发生率函数的扩散SIS传染病模型，深入探讨了复杂环境中传染病的传播机制。

Capasso和Serio^［6］提出了饱和发生率，其中\alpha是拥挤效应或心理因素的影响参数，表示感染传染病的感染人数，该发生率函数考虑了由于拥挤效应对传染病传播产生的抑制作用，并且当I增大到一定值时，发生率函数具有饱和状态。Shi等^［7］研究了具有饱和发生率和时间延迟的分数阶HIV模型，验证了模型在地方性平衡点和无病平衡点处的局部稳定性，使用时滞作为参数，研究一系列临界值时可能会出现Hopf分支，并验证其正确性。Mao等^［8］研究了一个考虑捕食者疾病和饱和发病率的延迟生态流行病学模型，利用迭代法推导出了无病平面平衡点和正平衡点全局稳定性的充分条件。Liu^［9］探讨了一个包含饱和治疗函数和非线性脉冲疫苗接种策略的SIR模型，利用分支理论，揭示了模型中跨临界分支与翻转分支的存在，并对传染病周期解的稳定性进行了详尽讨论。Ji等^［10］对随机SEIS传染病模型进行研究，该模型同样纳入了饱和发病率和时间延迟要素。他们证明了模型全局正解的存在性与唯一性，并通过构造Lyapunov函数，提出了确保地方性稳态随机稳定性的充分条件。

Ezekiel等^［11］对一类SIR传染病模型的动力学行为进行研究，分析了无病平衡点和正平衡点的局部稳定性，进一步讨论正平衡点处的Hopf分支的方向及稳定性。Attaullah等^［12］分析了其无病平衡点和正平衡点的局部稳定性以及全局稳定性，但在文中未讨论在正平衡点处产生的Hopf分支方向、周期解以及稳定性。因此，本文对此进行深入研究。

Attaullah等^［12］提出以下模型：

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\alpha (\mathrm{1}-\mu )\beta -(\alpha +\eta )S(t)-\frac{{\lambda }_{\mathrm{1}}S(t)I\left(t\right)}{\mathrm{1}+{\xi }_{\mathrm{1}}I\left(t\right)}+\gamma V(t),\\ \frac{\mathrm{d}V(t)}{\mathrm{d}t}=\alpha \beta \mu +\eta S(t)-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}V(t)I(t)-(\alpha +\gamma )V(t),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}=\frac{{\lambda }_{\mathrm{1}}S(t)I\left(t\right)}{\mathrm{1}+{\xi }_{\mathrm{1}}I\left(t\right)}+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}V(t)I(t)-(\alpha +\varphi +\phi )I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}=\varphi I(t)-\alpha R(t)\mathrm{。}\end{array}\right.（1）

依据微分方程解的存在性与非负性、各参数的生物学定义和真实的传染病特征数据确定参数的范围，使参数的选取具有实际意义。本文各参数含义以及参数取值范围见表1。

研究表明，传染病模型中引入时滞效应可有效刻画系统的动态演化特征，涵盖疫苗接种免疫滞后、感染潜伏期等延迟现象。汪雨琴等^［13］分析了一类具有时滞效应和Logistic增长的SIS传染病模型，基于Hopf分支理论，获得了地方病平衡点在不稳定时发生Hopf分支的条件。豆中丽^［14］对具有阶段结构的分数阶时滞生态流行病模型的稳定性进行研究，探讨在地方病平衡点附近存在Hopf分支的充分条件。因此，在模型中引入时滞机制不仅高度契合传染病传播的现实情境，具备显著的实践指导意义；更因其能够深入揭示系统复杂动态演化规律，而彰显出兼具理论探索深度与实际应用广度的研究价值。针对现有传染病模型未充分考量病毒传播时滞特性的局限，本研究在感染项中引入时滞参数，构建时滞动力学模型。疾病传播过程中，个体从感染到具备传染性存在时间间隔（如潜伏期），时滞\tau精准刻画了这一延迟，使模型更贴合“感染后需经历特定时段才具有传播能力”的真实传播流程。改进后的模型通过量化病毒传播周期与历史感染数据的动态关联，更精准地剖析潜伏期演化、免疫滞后等传播机制，进而揭示传播链的周期性波动规律。

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\alpha (\mathrm{1}-\mu )\beta -(\alpha +\eta )S(t)-\frac{{\lambda }_{\mathrm{1}}S(t)I\left(t\right)}{\mathrm{1}+{\xi }_{\mathrm{1}}I\left(t\right)}+\gamma V(t),\\ \frac{\mathrm{d}V(t)}{\mathrm{d}t}=\alpha \beta \mu +\eta S(t)-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}V(t)I(t-\tau )-(\alpha +\gamma )V(t),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}=\frac{{\lambda }_{\mathrm{1}}S(t)I\left(t\right)}{\mathrm{1}+{\xi }_{\mathrm{1}}I\left(t\right)}+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}V(t)I(t-\tau )-(\alpha +\varphi +\phi )I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}=\varphi I(t)-\alpha R(t)\mathrm{。}\end{array}\right.（2）

1 平衡点与局部稳定性

根据所研究的问题，先给出一系列假设：

（H1）P_{\mathrm{1}}+P_{\mathrm{3}}>\mathrm{0},(P_{\mathrm{2}}+P_{\mathrm{4}})(P_{\mathrm{1}}+P_{\mathrm{3}})-P_{\mathrm{6}}>\mathrm{0}；

（H2） k_{\mathrm{0}}\ge \mathrm{0},\mathrm{\Delta }=k_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{k}_{\mathrm{1}}\le \mathrm{0}；

（H3） k_{\mathrm{0}}\ge \mathrm{0},\mathrm{\Delta }=k_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{k}_{\mathrm{1}}\ge \mathrm{0},z^{\mathrm{*}}>\mathrm{0},f(z)\le \mathrm{0}；

（H4）；

（H5） f^{\mathrm{\text{'}}}(z)\ne \mathrm{0}。

为研究系统的动力学行为，首先确定系统的平衡点。无病平衡点是指感染者数量为0的情况，即I=\mathrm{0}，可得R=\mathrm{0}，下面求解无病平衡点E_{\mathrm{0}}。

令方程

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=\alpha (\mathrm{1}-\mu )\beta -(\alpha +\eta )S+\gamma V=\mathrm{0},\\ \frac{\mathrm{d}V}{\mathrm{d}t}=\alpha \beta \mu +\eta S-(\alpha +\gamma )V=\mathrm{0}\mathrm{。}\end{array}\right.

即

\left\{\begin{array}{l}(\alpha +\eta )S-\gamma V=\alpha (\mathrm{1}-\mu )\beta ,\\ -\eta S+(\alpha +\gamma )V=\alpha \beta \mu \mathrm{。}\end{array}\right.

用矩阵表示为

\left(\begin{array}{cc}\alpha +\eta & -\gamma \\ -\eta & \alpha +\gamma \end{array}\right)\left(\begin{array}{c}S\\ V\end{array}\right)=\left(\begin{array}{c}\alpha (\mathrm{1}-\mu )\beta \\ \alpha \beta \mu \end{array}\right)，

计算行列式

\mathrm{\Delta }=(\alpha +\eta )(\alpha +\gamma )-\eta \gamma =\alpha (\alpha +\eta +\gamma )\mathrm{。}

利用克莱姆法则可得

\begin{array}{l}S=\frac{\left|\begin{array}{cc}\alpha (\mathrm{1}-\mu )\beta & -\gamma \\ \alpha \beta \mu & \alpha +\gamma \end{array}\right|}{\mathrm{\Delta }}=\frac{\alpha (\mathrm{1}-\mu )\beta (\alpha +\gamma )+\alpha \beta \mu \gamma }{\alpha (\alpha +\eta +\gamma )}=\frac{\beta (\alpha +\gamma -\mu \alpha )}{\alpha +\eta +\gamma },\\ V=\frac{\left|\begin{array}{cc}\alpha +\eta & \alpha (\mathrm{1}-\mu )\beta \\ -\eta & \alpha \beta \mu \end{array}\right|}{\mathrm{\Delta }}=\frac{\alpha \beta \mu (\alpha +\eta )+\alpha \eta (\mathrm{1}-\mu )\beta }{\alpha (\alpha +\eta +\gamma )}=\frac{\beta (\mu \alpha +\eta )}{\alpha +\eta +\gamma }\mathrm{。}\end{array}

由此得知，此系统的无病平衡点E_{\mathrm{0}}=(S_{\mathrm{0}},V_{\mathrm{0}},I_{\mathrm{0}},R_{\mathrm{0}})为

E_{\mathrm{0}}=(S_{\mathrm{0}},V_{\mathrm{0}},I_{\mathrm{0}},R_{\mathrm{0}})=(\frac{\beta (\alpha +\gamma -\mu \alpha )}{\alpha +\eta +\gamma },\frac{\beta (\mu \alpha +\eta )}{\alpha +\eta +\gamma },\mathrm{0,0})\mathrm{。}

基本再生数，即疾病在起始传播时，单一感染者于感染周期内，平均致使易感个体被传染的数量，是衡量疾病传播力与防控策略成效的关键指标。由文献［12］知该系统的基本再生数为

R_{\mathrm{0}}=\frac{\beta [{\lambda }_{\mathrm{1}}(-\alpha \mu +\alpha +\gamma )+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}(\alpha \mu +\eta )]}{(\alpha +\gamma +\eta )(\alpha +\varphi +\phi )}。

当R_{\mathrm{0}}\le \mathrm{1}时，意味着无病平衡点是稳定的，疾病会逐渐消失；当R_{\mathrm{0}}>\mathrm{1}时，系统还存在地方病平衡点E^{\mathrm{*}}=(S^{\mathrm{*}},V^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}})，疾病将会在社会中永久存在。下面求解地方病平衡点E^{\mathrm{*}}。

令方程\varphi I^{\mathrm{*}}-\alpha R^{\mathrm{*}}=\mathrm{0}推得R^{\mathrm{*}}=\frac{\varphi I^{\mathrm{*}}}{\alpha }。

将系统（1）各方程加和，代入总人口N^{\mathrm{*}}=S^{\mathrm{*}}+V^{\mathrm{*}}+I^{\mathrm{*}}+R^{\mathrm{*}}，求导得\frac{\mathrm{d}{N}^{\mathrm{*}}}{\mathrm{d}t}=\alpha \beta -\alpha N^{\mathrm{*}}-\phi I^{\mathrm{*}}=\mathrm{0},解得N^{\mathrm{*}}=\beta -\frac{\phi }{\alpha }{I}^{\mathrm{*}}，即N^{\mathrm{*}}=S^{\mathrm{*}}+V^{\mathrm{*}}+I^{\mathrm{*}}+R^{\mathrm{*}}=\beta -\frac{\phi }{\alpha }{I}^{\mathrm{*}}，S^{\mathrm{*}}+V^{\mathrm{*}}=\beta -\frac{\alpha +\phi -\varphi }{\alpha }{I}^{\mathrm{*}}。

对方程\frac{{\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}{I}^{\mathrm{*}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-(\alpha +\varphi +\phi )I^{\mathrm{*}}=\mathrm{0}，同时除以I^{\mathrm{*}}(I^{\mathrm{*}}\ne \mathrm{0})，得\frac{{\lambda }_{\mathrm{1}}S}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}=\alpha +\varphi +\phi。

联立方程\left\{\begin{array}{l}{S}^{\mathrm{*}}+V^{\mathrm{*}}=\beta -\frac{\alpha +\varphi +\phi }{\alpha }{I}^{\mathrm{*}},\\ \frac{{\lambda }_{\mathrm{1}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}{S}^{\mathrm{*}}+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}=\alpha +\varphi +\phi \mathrm{。}\end{array}\right.

解得

S^{\mathrm{*}}=\beta -\frac{\alpha +\varphi +\phi }{\alpha }{I}^{\mathrm{*}}-V^{\mathrm{*}}，

V^{\mathrm{*}}=\beta \cdot \frac{\eta +\alpha \mu -\frac{\eta k}{\beta }I}{{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}I+\alpha +\gamma +\eta }=\beta \cdot \frac{\eta +\alpha \mu -\eta \cdot \frac{\alpha +\varphi +\phi }{\alpha }I}{{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}I+\alpha +\gamma +\eta }\mathrm{。}

将S^{\mathrm{*}}\mathrm{、}{V}^{\mathrm{*}}表达式代入方程，可得I^{\mathrm{*}}满足二次方程a{I^{\mathrm{*}}}^{\mathrm{2}}+b{I}^{\mathrm{*}}+c=\mathrm{0}，其中

a=-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}\left[{\lambda }_{\mathrm{1}}k+\eta k{\xi }_{\mathrm{1}}+(\alpha +\varphi +\phi ){\xi }_{\mathrm{1}}\right]，

b={\lambda }_{\mathrm{1}}\left(\beta {\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}+k(\alpha +\gamma )\right)+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}\left(\beta (\eta +\alpha \mu ){\xi }_{\mathrm{1}}-\eta k\right)-(\alpha +\varphi +\phi )[{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}+(\alpha +\gamma +\eta ){\xi }_{\mathrm{1}}]，

c={\lambda }_{\mathrm{1}}\beta (\alpha +\gamma +\alpha \mu )+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}\beta (\eta +\alpha \mu )-(\alpha +\varphi +\phi )(\alpha +\gamma +\eta )。

因此，可得系统（1）在地方病平衡点E^{\mathrm{*}}=(S^{\mathrm{*}},V^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}})。

系统（2）在平衡点处的线性系统（3）为

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=-(\alpha +\eta )-\frac{{\lambda }_{\mathrm{1}}{I}^{\mathrm{*}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}S(t)+\gamma V(t)-\frac{{\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}}{(\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}{)}^{\mathrm{2}}}I(t),\\ \frac{\mathrm{d}V(t)}{\mathrm{d}t}=\eta S(t)+-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{I}^{\mathrm{*}}-(\alpha +\gamma )V(t)-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}I(t-\tau ),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}=\frac{{\lambda }_{\mathrm{1}}{I}^{\mathrm{*}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}S(t)+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{I}^{\mathrm{*}}V(t)+\frac{{\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}}{(\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}{)}^{\mathrm{2}}}-(\alpha +\varphi +\phi )I(t)+{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}I(t-\tau ),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}=\varphi I(t)-\alpha R(t)\mathrm{。}\end{array}\right.（3）

令

a_{\mathrm{11}}=-(\alpha +\eta )-\frac{{\lambda }_{\mathrm{1}}{I}^{\mathrm{*}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}，a_{\mathrm{12}}=\gamma，a_{\mathrm{21}}=\eta，a_{\mathrm{22}}=-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{I}^{\mathrm{*}}-(\alpha +\gamma )，

a_{\mathrm{31}}=\frac{{\lambda }_{\mathrm{1}}{I}^{\mathrm{*}}}{\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}}，a_{\mathrm{32}}={\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{I}^{\mathrm{*}}，a_{\mathrm{33}}=\frac{{\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}}{(\mathrm{1}+{\xi }_{\mathrm{1}}{I}^{\mathrm{*}}{)}^{\mathrm{2}}}-(\alpha +\varphi +\phi )，a_{\mathrm{43}}=\varphi，a_{\mathrm{44}}=-\alpha，

b_{\mathrm{23}}=-{\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}，b_{\mathrm{33}}={\lambda }_{\mathrm{2}}{\xi }_{\mathrm{2}}{V}^{\mathrm{*}}。

系统（3）在平衡点处的雅可比矩阵为

\mathit{J}(E^{\mathrm{*}})=\left(\begin{array}{cccc}{a}_{\mathrm{11}}& a_{\mathrm{12}}& a_{\mathrm{13}}& \mathrm{0}\\ a_{\mathrm{21}}& a_{\mathrm{22}}& b_{\mathrm{23}}& \mathrm{0}\\ a_{\mathrm{31}}& a_{\mathrm{32}}& a_{\mathrm{33}}+b_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& a_{\mathrm{43}}& a_{\mathrm{44}}\end{array}\right)，

特征方程为(\lambda -a_{\mathrm{44}})[{\lambda }^{\mathrm{3}}+P_{\mathrm{1}}{\lambda }^{\mathrm{2}}+P_{\mathrm{2}}\lambda +P_{\mathrm{3}}+(P_{\mathrm{4}}{\lambda }^{\mathrm{2}}+P_{\mathrm{5}}\lambda +P_{\mathrm{6}}){\mathrm{e}}^{-\lambda \tau }]=\mathrm{0}。

显而易见，方程的一个根为\lambda =-a_{\mathrm{44}}，即\lambda =\alpha，由此可研究方程

{\lambda }^{\mathrm{3}}+P_{\mathrm{1}}{\lambda }^{\mathrm{2}}+P_{\mathrm{2}}\lambda +P_{\mathrm{3}}+(P_{\mathrm{4}}{\lambda }^{\mathrm{2}}+P_{\mathrm{5}}\lambda +P_{\mathrm{6}}){\mathrm{e}}^{-\lambda \tau }=\mathrm{0}。（4）

其中，P_{\mathrm{1}}=-(a_{\mathrm{11}}+a_{\mathrm{22}}+a_{\mathrm{33}}),P_{\mathrm{2}}=a_{\mathrm{11}}{a}_{\mathrm{22}}+a_{\mathrm{11}}{a}_{\mathrm{33}}+a_{\mathrm{22}}{a}_{\mathrm{33}}-a_{\mathrm{12}}{a}_{\mathrm{21}}-a_{\mathrm{13}}{a}_{\mathrm{31}},P_{\mathrm{3}}=a_{\mathrm{12}}{a}_{\mathrm{21}}{a}_{\mathrm{33}}+a_{\mathrm{13}}{a}_{\mathrm{22}}{a}_{\mathrm{31}}-a_{\mathrm{13}}{a}_{\mathrm{21}}{a}_{\mathrm{32}}-a_{\mathrm{11}}{a}_{\mathrm{22}}{a}_{\mathrm{33}},P_{\mathrm{4}}=-b_{\mathrm{33}},P_{\mathrm{5}}=a_{\mathrm{11}}{b}_{\mathrm{33}}+a_{\mathrm{22}}{b}_{\mathrm{33}}-a_{\mathrm{32}}{b}_{\mathrm{23}},P_{\mathrm{6}}=-a_{\mathrm{11}}{a}_{\mathrm{22}}{b}_{\mathrm{33}}+a_{\mathrm{11}}{a}_{\mathrm{32}}{b}_{\mathrm{23}}+a_{\mathrm{12}}{a}_{\mathrm{21}}{b}_{\mathrm{33}}-a_{\mathrm{12}}{a}_{\mathrm{31}}{b}_{\mathrm{23}}。

当\tau =\mathrm{0}时，

{\lambda }^{\mathrm{3}}+(P_{\mathrm{1}}+P_{\mathrm{4}}){\lambda }^{\mathrm{2}}+(P_{\mathrm{2}}+P_{\mathrm{5}})\lambda +P_{\mathrm{3}}+P_{\mathrm{6}}=\mathrm{0}。（5）

依据Routh-Hurwitz判据，确定方程（5）的所有根都具有负实部的等价条件为（H1），故（H1）得到满足，根据稳定性理论，正平衡点是局部渐近稳定的。

当\tau >\mathrm{0}时，假设\lambda =\mathrm{i}\omega是方程（4）的纯虚根，将其代入方程可得

{(\mathrm{i}\omega )}^{\mathrm{3}}+P_{\mathrm{1}}{(\mathrm{i}\omega )}^{\mathrm{2}}+P_{\mathrm{2}}(\mathrm{i}\omega )+[P_{\mathrm{4}}{(\mathrm{i}\omega )}^{\mathrm{2}}+P_{\mathrm{5}}(\mathrm{i}\omega )+P_{\mathrm{6}}]{\mathrm{e}}^{-(\mathrm{i}\omega )\tau }+P_{\mathrm{3}}=\mathrm{0}。

分离虚实部得

\left\{\begin{array}{l}{P}_{\mathrm{5}}\omega \mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -(P_{\mathrm{4}}{\omega }^{\mathrm{2}}-P_{\mathrm{6}})\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau =P_{\mathrm{2}}\omega +{\omega }^{\mathrm{3}},\\ P_{\mathrm{5}}\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau +(P_{\mathrm{4}}{\omega }^{\mathrm{2}}-P_{\mathrm{6}})\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau =-P_{\mathrm{1}}{\omega }^{\mathrm{2}}+P_{\mathrm{3}}\mathrm{。}\end{array}\right.（6）

将式（6）两端同时平方相加，得

{\omega }^{\mathrm{6}}+k_{\mathrm{2}}{\omega }^{\mathrm{4}}+k_{\mathrm{1}}{\omega }^{\mathrm{2}}+k_{\mathrm{0}}=\mathrm{0}，（7）

其中k_{\mathrm{2}}=P_{\mathrm{1}}^{\mathrm{2}}-\mathrm{2}{P}_{\mathrm{2}}-P_{\mathrm{4}}^{\mathrm{2}},k_{\mathrm{1}}=-\mathrm{2}{P}_{\mathrm{1}}{P}_{\mathrm{3}}+P_{\mathrm{2}}^{\mathrm{2}}-P_{\mathrm{5}}^{\mathrm{2}}+\mathrm{2}{P}_{\mathrm{4}}{P}_{\mathrm{6}},k_{\mathrm{0}}=P_{\mathrm{3}}^{\mathrm{2}}-P_{\mathrm{6}}^{\mathrm{2}}。

令z={\omega }^{\mathrm{2}}，代入式（7）得

z^{\mathrm{3}}+k_{\mathrm{2}}{z}^{\mathrm{2}}+k_{\mathrm{1}}z+k_{\mathrm{0}}=\mathrm{0}。（8）

假设f(z)=z^{\mathrm{3}}+k_{\mathrm{2}}{z}^{\mathrm{2}}+k_{\mathrm{1}}z+k_{\mathrm{0}}，根据特征根的分布定理有以下结论产生：（1）若条件（H2）成立，无正实根；（2）若条件（H3）或（H4）成立，至少有一个正实根。

假设条件（H3）或（H4）总是得到满足，那么方程（8）存在三个正实根，分别为{\omega }_k(k=\mathrm{1,2},\mathrm{3})，其对应的始终临界值为{\tau }_k^j=\frac{\mathrm{1}}{{\omega }_k}(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{(-P_{\mathrm{5}}{\omega }_k)(P_{\mathrm{2}}{\omega }_k-{\omega }_k^{\mathrm{3}})+(-P_{\mathrm{1}}{\omega }_k^{\mathrm{2}}+P_{\mathrm{3}})(P_{\mathrm{4}}{\omega }_k^{\mathrm{2}}-P_{\mathrm{6}})}{(P_{\mathrm{4}}{\omega }_k^{\mathrm{2}}-P_{\mathrm{6}}{)}^{\mathrm{2}}+(-P_{\mathrm{5}}{\omega }_k{)}^{\mathrm{2}}}+\mathrm{2}j\mathrm{\pi })(k=\mathrm{1,2},\mathrm{3};j=\mathrm{0,1},\mathrm{2},\dots )，假设满足条件（H5），那么有横截性条件\mathrm{R}\mathrm{e}(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })\left|{}_{\lambda =\mathrm{i}\omega }\right.\ne \mathrm{0}成立。

下面对横截性条件进行证明，对于方程（5）左右两端进行求导，得

{(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })}^{-\mathrm{1}}=-\frac{\mathrm{3}{\lambda }^{\mathrm{2}}+\mathrm{2}{P}_{\mathrm{1}}\lambda +P_{\mathrm{2}}}{\lambda \left({\lambda }^{\mathrm{3}}+P_{\mathrm{1}}{\lambda }^{\mathrm{2}}+P_{\mathrm{2}}\lambda +P_{\mathrm{3}}\right)}+\frac{\mathrm{2}{P}_{\mathrm{4}}\lambda +P_{\mathrm{5}}}{\lambda \left(P_{\mathrm{4}}{\lambda }^{\mathrm{2}}+P_{\mathrm{5}}\lambda +P_{\mathrm{6}}\right)}-\frac{\tau }{\lambda }，

将\lambda =\mathrm{i}\omega代入{(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })}^{-\mathrm{1}}表达式，并对其进行虚实部分离，有

\begin{array}{l}\mathrm{R}\mathrm{e}{(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })}^{-\mathrm{1}}\left|{}_{\lambda =\mathrm{i}\omega }\right.=\frac{\mathrm{3}{\omega }^{\mathrm{4}}+\mathrm{2}{P}_{\mathrm{1}}^{\mathrm{2}}{\omega }^{\mathrm{2}}-\mathrm{4}{P}_{\mathrm{2}}{\omega }^{\mathrm{2}}+P_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}{P}_{\mathrm{1}}{P}_{\mathrm{3}}}{(P_{\mathrm{1}}{\omega }^{\mathrm{2}}-P_{\mathrm{3}}{)}^{\mathrm{2}}+({\omega }^{\mathrm{2}}+P_{\mathrm{2}}{\omega )}^{\mathrm{2}}}-\frac{\mathrm{2}{P}_{\mathrm{4}}^{\mathrm{2}}{\omega }^{\mathrm{2}}+P_{\mathrm{5}}^{\mathrm{2}}-\mathrm{2}{P}_{\mathrm{6}}{P}_{\mathrm{4}}}{P_{\mathrm{4}}^{\mathrm{2}}{\omega }^{\mathrm{2}}+(P_{\mathrm{6}}-P_{\mathrm{4}}{\omega }^{\mathrm{2}}{)}^{\mathrm{2}}}=\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{\mathrm{3}{\omega }^{\mathrm{4}}+\mathrm{2}(P_{\mathrm{1}}^{\mathrm{2}}-P_{\mathrm{4}}^{\mathrm{2}}){\omega }^{\mathrm{2}}-\mathrm{4}{P}_{\mathrm{2}}{\omega }^{\mathrm{2}}+P_{\mathrm{2}}^{\mathrm{2}}-P_{\mathrm{5}}^{\mathrm{2}}+\mathrm{2}{P}_{\mathrm{6}}{P}_{\mathrm{4}}-\mathrm{2}{P}_{\mathrm{1}}{P}_{\mathrm{3}}}{(P_{\mathrm{1}}{\omega }^{\mathrm{2}}-P_{\mathrm{3}}{)}^{\mathrm{2}}+({\omega }^{\mathrm{2}}+P_{\mathrm{2}}{\omega )}^{\mathrm{2}}}\mathrm{。}\end{array}

令f^{\text{'}}(z)=\mathrm{3}{\omega }^{\mathrm{4}}+\mathrm{2}(P_{\mathrm{1}}^{\mathrm{2}}-P_{\mathrm{3}}^{\mathrm{2}}){\omega }^{\mathrm{2}}-\mathrm{4}{P}_{\mathrm{2}}{\omega }^{\mathrm{2}}+P_{\mathrm{2}}^{\mathrm{2}}-P_{\mathrm{4}}^{\mathrm{2}}+\mathrm{2}{P}_{\mathrm{5}}{P}_{\mathrm{3}}-\mathrm{2}{P}_{\mathrm{1}}{P}_{\mathrm{6}}，有\mathrm{R}\mathrm{e}{(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })}^{-\mathrm{1}}\left|{}_{\lambda =\mathrm{i}\omega }\right.=\frac{f^{\text{'}}(z)}{(P_{\mathrm{1}}{\omega }^{\mathrm{2}}-P_{\mathrm{3}}{)}^{\mathrm{2}}+({\omega }^{\mathrm{2}}+P_{\mathrm{2}}{\omega )}^{\mathrm{2}}}\ne \mathrm{0}成立。因为\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left\{\mathrm{R}\mathrm{e}{(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })}_{{}_{\lambda =\mathrm{i}\omega }}^{-\mathrm{1}}\right\}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left\{{(\frac{\mathrm{d}(\mathrm{R}\mathrm{e}\lambda )}{\mathrm{d}\tau })}_{\lambda =\mathrm{i}\omega }\right\}，所以\mathrm{R}\mathrm{e}(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau })\left|{}_{\lambda =\mathrm{i}\omega }\right.\ne \mathrm{0}，横截性条件证毕。

综合以上的结果，可得定理1。

定理1 （1）当\tau =\mathrm{0}时，若条件（H1）满足，地方病平衡点E^{\mathrm{*}}是局部渐近稳定的；（2）当\tau >\mathrm{0}时，若条件（H2）满足，地方病平衡点E^{\mathrm{*}}是局部渐近稳定的；（3）当\tau \in [\mathrm{0},{\tau }_{\mathrm{0}})时，若条件（H3）或（H4）满足，其地方病平衡点E^{\mathrm{*}}是局部渐近稳定的；（4）当\tau >{\tau }_{\mathrm{0}}时，若（H3）（H4）（H5）满足，地方病平衡点E^{\mathrm{*}}是不稳定的；（5）当\tau ={\tau }_{\mathrm{0}}时，若（H3）（H4）（H5）满足，其在地方病平衡点E^{\mathrm{*}}处产生Hopf分支。

2 Hopf分支方向及其稳定性

在\tau ={\tau }_{\mathrm{0}}的前提条件下，运用Hassard^［15］提到的中心流形定理和规范性方法，确定系统（2）的局部Hopf分支方向和周期解稳定性的表示形式。

令\mathit{u}(t)=(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t),u_{\mathrm{3}}(t),u_{\mathrm{4}}{(t))}^{\mathrm{T}}\in {\mathbb{R}}^{\mathrm{4}}，其中u_{\mathrm{1}}(t)=x_{\mathrm{1}}(\tau t)，u_{\mathrm{2}}(t)=x_{\mathrm{2}}(\tau t)，u_{\mathrm{3}}(t)=x_{\mathrm{3}}(\tau t)，u_{\mathrm{4}}(t)=x_{\mathrm{4}}(\tau t)，\tau ={\tau }_{\mathrm{0}}+\mu，\mu \in \mathbb{R}，则系统（2）在C=C([-\mathrm{1,0}],{\mathbb{R}}^{\mathrm{4}})上变为一般的泛函微分方程

\dot{\mathit{u}}(t)={\mathit{L}}_{\mu }({\mathit{u}}_t)+F(\mu ,{\mathit{u}}_t)。（9）

其中L_{\mu }:C\to {\mathbb{R}}^{\mathrm{4}}，F:\mathbb{R}\times C\to {\mathbb{R}}^{\mathrm{4}}分别由以下形式给出：

{\mathit{L}}_{\mu }(\mathit{\phi })=({\tau }_{\mathrm{0}}+\mu )\left(\begin{array}{cccc}{a}_{\mathrm{11}}& a_{\mathrm{12}}& a_{\mathrm{13}}& \mathrm{0}\\ a_{\mathrm{21}}& a_{\mathrm{22}}& \mathrm{0}& \mathrm{0}\\ a_{\mathrm{31}}& a_{\mathrm{32}}& a_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& a_{\mathrm{43}}& a_{\mathrm{44}}\end{array}\right)\left(\begin{array}{c}{\phi }_{\mathrm{1}}(\mathrm{0})\\ {\phi }_{\mathrm{2}}(\mathrm{0})\\ {\phi }_{\mathrm{3}}(\mathrm{0})\\ {\phi }_{\mathrm{4}}(\mathrm{0})\end{array}\right)+({\tau }_{\mathrm{0}}+\mu )\left(\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& b_{\mathrm{23}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& b_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}{\phi }_{\mathrm{1}}(-\mathrm{1})\\ {\phi }_{\mathrm{2}}(-\mathrm{1})\\ {\phi }_{\mathrm{3}}(-\mathrm{1})\\ {\phi }_{\mathrm{4}}(-\mathrm{1})\end{array}\right)，

\mathit{F}(\mu ,\mathit{\phi })=({\tau }_{\mathrm{0}}+\mu )(F_{\mathrm{1}}(\mu ,\mathit{\phi }),F_{\mathrm{2}}(\mu ,\mathit{\phi }),\hspace{1em}{F}_{\mathrm{3}}(\mu ,\mathit{\phi }),\hspace{1em}{F}_{\mathrm{4}}{(\mu ,\mathit{\phi }))}^{\mathrm{T}}。（10）

其中

\mathit{\phi }=({\phi }_{\mathrm{1}},{\phi }_{\mathrm{2}},{\phi }_{\mathrm{3}},{\phi }_{\mathrm{4}})\in C([-\mathrm{1,0}],{\mathbb{R}}^{\mathrm{4}})，

F_{\mathrm{1}}(\mu ,\mathit{\phi })=b_{\mathrm{11}}{\phi }_{\mathrm{1}}(\mathrm{0}){\phi }_{\mathrm{3}}(\mathrm{0})，

F_{\mathrm{2}}(\mu ,\mathit{\phi })=c_{\mathrm{11}}{\phi }_{\mathrm{2}}(\mathrm{0}){\phi }_{\mathrm{3}}(-\mathrm{1})，

F_{\mathrm{3}}(\mu ,\mathit{\phi })=d_{\mathrm{11}}{\phi }_{\mathrm{3}}^{\mathrm{2}}(\mathrm{0})+d_{\mathrm{12}}{\phi }_{\mathrm{3}}(\mathrm{0}){\phi }_{\mathrm{1}}(\mathrm{0})+d_{\mathrm{13}}{\phi }_{\mathrm{3}}(-\mathrm{1}){\phi }_{\mathrm{2}}(\mathrm{0})。

故应用Riesz表示定理，对于\mathit{\phi }\in \left[-\mathrm{1,0}\right]存在有界变分函数\mathit{\eta }(\theta ,\mu )，使得

{\mathit{L}}_{\mu }(\phi )={\int }_{-\mathrm{1}}^{\mathrm{0}}d\mathit{\eta }(\theta ,\mu )\mathit{\phi }(\theta ),\hspace{1em}\phi \in C\left(\left[-\mathrm{1,0}\right],{\mathbb{R}}^{\mathrm{4}}\right)，

其中

\mathit{\eta }(\theta ,\mu )=({\tau }_{\mathrm{0}}+\mu )\left(\begin{array}{cccc}{a}_{\mathrm{11}}& a_{\mathrm{12}}& a_{\mathrm{13}}& \mathrm{0}\\ a_{\mathrm{21}}& a_{\mathrm{22}}& \mathrm{0}& \mathrm{0}\\ a_{\mathrm{31}}& a_{\mathrm{32}}& a_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& a_{\mathrm{43}}& a_{\mathrm{44}}\end{array}\right)\delta (\theta )+({\tau }_{\mathrm{0}}+\mu )\left(\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& b_{\mathrm{23}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& b_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\delta (\theta +\mathrm{1})，

\delta (\theta )是Dirac-delta函数。

对于\mathit{\phi }\in C\left(\left[-\mathrm{1,0}\right],{\mathbb{R}}^{\mathrm{4}}\right)，定义

\mathit{A}(\mu )\mathit{\phi }(\theta )=\left\{\begin{array}{l}\dot{\mathit{\phi }}(\theta ),\theta \in [-\mathrm{1,0}),\\ {\int }_{-\mathrm{1}}^{\mathrm{0}}\mathrm{d}\mathit{\eta }(\theta ,\mu )\mathit{\phi }(\theta ),\theta =\mathrm{0}\mathrm{。}\end{array}\right.

\mathit{R}(\mu )\mathit{\phi }(\theta )=\left\{\begin{array}{l}{(\mathrm{0,0})}^{\mathrm{T}},\theta \in [-\mathrm{1,0}),\\ \mathit{F}(\mu ,\theta ),\hspace{1em}\theta =\mathrm{0}\mathrm{。}\end{array}\right.

于是系统（10）可以改写为抽象的微分方程

\dot{\mathit{u}}(t)=\mathit{A}(\mu ){\mathit{u}}_t+\mathit{R}(\mu ){\mathit{u}}_t，（11）

其中，u=(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}},u_{\mathrm{4}});u_t(\theta )=u(t+\theta ),\theta \in [-\mathrm{1,0}]。

对于\mathit{\psi }\in C\left(\left[-\mathrm{1,0}\right],({\mathbb{R}}^{\mathrm{4}}{)}^{\mathrm{*}}\right)，定义A=A(\mathrm{0})的伴随算子A^{\mathrm{*}}：

A^{\mathrm{*}}\mathit{\psi }(s)=\left\{\begin{array}{l}-\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}s},s\in (\mathrm{0,1}],\\ {\int }_{-\mathrm{1}}^{\mathrm{0}}\mathit{\psi }(-s)\mathrm{d}\mathit{\eta }(s,\mu ),s=\mathrm{0}\mathrm{。}\end{array}\right.

对于\mathit{\phi }\in C\left(\left[-\mathrm{1,0}\right],{\mathbb{R}}^{\mathrm{4}}\right)，\mathit{\psi }\in C\left(\left[-\mathrm{1,0}\right],\hspace{1em}({\mathbb{R}}^{\mathrm{4}}{)}^{\mathrm{*}}\right)，为了使算子A和伴随算子A^{\mathrm{*}}的特征向量规范化，定义双线性形式：

〈\mathit{\psi }(s),\mathit{\phi }(\theta )〉={\overline{\mathit{\psi }}}^{\mathit{T}}(\mathrm{0})\mathit{\phi }(\mathrm{0})-{\int }_{\theta =-\mathrm{1}}^{\mathrm{0}}{\int }_{\xi =\mathrm{0}}^{\theta }{\overline{\mathit{\psi }}}^{\mathrm{T}}(\xi -\theta )\mathrm{d}\mathit{\eta }(\theta )\mathit{\phi }(\xi )\mathrm{d}\xi。

令\mathit{q}(\theta )和{\mathit{q}}^{\mathrm{*}}(\theta )分别为A和A^{\mathrm{*}}对于\lambda (\mathrm{0})=\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}和\overline{\lambda }(\mathrm{0})=-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}的特征向量，则有\mathit{A}(\mathrm{0})\mathit{q}(\theta )=\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\mathit{q}(\theta ),{\mathit{A}}^{\mathrm{*}}(\mathrm{0}){\mathit{q}}^{\mathrm{*}}(s)=-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}{\mathit{q}}^{\mathrm{*}}(s)，通过计算可得：

\mathit{q}(\theta )=\left(\begin{array}{l}\mathrm{1}\\ {\alpha }_{\mathrm{2}}\\ {\alpha }_{\mathrm{3}}\\ {\alpha }_{\mathrm{4}}\end{array}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta },\theta \in [-\mathrm{1,0})，{\mathit{q}}^{\mathrm{*}}(s)=B\left(\begin{array}{l}\mathrm{1}\\ {\beta }_{\mathrm{2}}\\ {\beta }_{\mathrm{3}}\\ {\beta }_{\mathrm{4}}\end{array}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}s},s\in (\mathrm{0,1}]，

其中

{\alpha }_{\mathrm{2}}=\frac{\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}}-a_{\mathrm{13}}{\alpha }_{\mathrm{3}}}{a_{\mathrm{12}}},{\alpha }_{\mathrm{3}}=\frac{a_{\mathrm{12}}{a}_{\mathrm{21}}-\mathrm{i}{\omega }_{\mathrm{0}}(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}})-a_{\mathrm{22}}(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}})}{-(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{22}})a_{\mathrm{31}}-a_{\mathrm{12}}{b}_{\mathrm{23}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}},{\alpha }_{\mathrm{4}}=\frac{a_{\mathrm{43}}{\alpha }_{\mathrm{3}}}{\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{44}}}，

{\beta }_{\mathrm{2}}=\frac{-\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}}-a_{\mathrm{31}}{\beta }_{\mathrm{3}}}{a_{\mathrm{21}}},{\beta }_{\mathrm{3}}=\frac{a_{\mathrm{12}}{a}_{\mathrm{21}}-\mathrm{i}{\omega }_{\mathrm{0}}(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}})-a_{\mathrm{22}}(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}})}{-(\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{22}})a_{\mathrm{31}}-a_{\mathrm{21}}{a}_{\mathrm{32}}},{\beta }_{\mathrm{4}}=\mathrm{0}，

B=\frac{\mathrm{1}}{\mathrm{1}+{\overline{\beta }}_{\mathrm{2}}{\alpha }_{\mathrm{2}}+{\overline{\beta }}_{\mathrm{3}}{\alpha }_{\mathrm{3}}-{\tau }_{\mathrm{0}}(b_{\mathrm{23}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\overline{\beta }}_{\mathrm{2}}{\alpha }_{\mathrm{3}}+b_{\mathrm{23}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\overline{\beta }}_{\mathrm{3}}{\alpha }_{\mathrm{3}}){\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}}，

则有\mathit{q}(\theta )和{\mathit{q}}^{\mathrm{*}}(\theta )满足〈{\mathit{q}}^{\mathrm{*}}(s),\mathit{q}(\theta )〉=\mathrm{1},〈{\mathit{q}}^{*}(s),\overline{\mathit{q}}(\theta )〉=\mathrm{0}。

令X_t为方程（12）在\mu =\mathrm{0}时的解，定义

z(t)=〈q^{*},X_t〉，

\mathit{W}(t,\theta )={\mathit{X}}_t(\theta )-\mathrm{2}\mathrm{R}\mathrm{e}\left\{z(t)\mathit{q}(\theta )\right\}={\mathit{X}}_t(\theta )-z(t)\mathit{q}(\theta )-\overline{z}(t)\overline{\mathit{q}}(\theta )。（12）

在中心流形C_{\mathrm{0}}上，有

\mathit{W}(t,\theta )=\mathit{W}(z(t),\overline{z}(t),\theta )={\mathit{W}}_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\mathit{W}}_{\mathrm{11}}(\theta )z\overline{z}+{\mathit{W}}_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\cdots，（13）

其中，z、\overline{z}为中心流形C_{\mathrm{0}}在{\mathit{q}}^{\mathrm{*}}和{\overline{\mathit{q}}}^{\mathrm{*}}方向上的状态坐标。考虑X_t和W的实解，在\mu =\mathrm{0}时，有

\dot{z}(t)=\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}z(t)+{{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})\left(\begin{array}{l}{F}_{\mathrm{1}}(\mu ,\mathit{\phi })\\ F_{\mathrm{2}}(\mu ,\mathit{\phi })\\ F_{\mathrm{3}}(\mu ,\mathit{\phi })\\ F_{\mathrm{4}}(\mu ,\phi )\end{array}\right)=\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}z(t)+g(z,\overline{z})，

其中，\dot{z}(t)=\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}z(t)+g(z,\overline{z})。根据中心流形理论和规范型变换，可将中心流形上的g(z,\overline{z})展开为关于z和\overline{z}的幂级数形式：

g(z,\overline{z})=g_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{11}}(\theta )z\overline{z}+g_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{21}}(\theta )\frac{z^{\mathrm{2}}\overline{z}}{\mathrm{2}}+\cdots。（14）

结合中心流形解的定义W(t,\theta )以及g(z,\overline{z})可得：

{\mathit{X}}_t(\theta )=\mathit{W}(t,\theta )+\mathrm{2}\mathrm{R}\mathrm{e}\left\{z(t)\mathit{q}(\theta )\right\}，

即

{\mathit{X}}_t(\theta )={\mathit{W}}_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\mathit{W}}_{\mathrm{11}}(\theta )z\overline{z}+{\mathit{W}}_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+(\mathrm{1},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}{)}^{\mathrm{T}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }z+(\mathrm{1},{\overline{\alpha }}_{\mathrm{2}},{\overline{\alpha }}_{\mathrm{3}}{)}^{\mathrm{T}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }\overline{z}+\cdots。

结合式（11）可得g(z,\overline{z})={{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})f_{\mathrm{0}}(z,\overline{z})，展开为

\begin{array}{l}g(z,\overline{z})={\tau }_{\mathrm{0}}\overline{B}(r_{\mathrm{11}}+{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{21}}+{\overline{\beta }}_{\mathrm{3}}{r}_{\mathrm{31}})\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\tau }_{\mathrm{0}}\overline{B}(r_{\mathrm{12}}+{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{22}}+r_{\mathrm{32}}\overline{\beta })\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\\ {\tau }_{\mathrm{0}}\overline{B}({\overline{r}}_{\mathrm{11}}+{\overline{\beta }}_{\mathrm{2}}{\overline{r}}_{\mathrm{21}}+{\overline{r}}_{\mathrm{31}}{\overline{\beta }}_{\mathrm{3}})z\overline{z}+\mathrm{2}{\tau }_{\mathrm{0}}\overline{B}(r_{\mathrm{13}}+{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{23}}+r_{\mathrm{33}}{\overline{\beta }}_{\mathrm{3}})\frac{z^{\mathrm{2}}\overline{z}}{\mathrm{2}}+\cdots \mathrm{。}\end{array}

其中，

r_{\mathrm{11}}=\mathrm{2}{b}_{\mathrm{11}}，r_{\mathrm{12}}=\mathrm{2}{b}_{\mathrm{11}}，r_{\mathrm{13}}=b_{\mathrm{11}}[W_{\mathrm{11}}^{\mathrm{3}}(\mathrm{0})+W_{\mathrm{20}}^{\mathrm{3}}(\mathrm{0})+W_{\mathrm{11}}^{\mathrm{1}}(\mathrm{0})+W_{\mathrm{20}}^{\mathrm{1}}(\mathrm{0})]，

r_{\mathrm{21}}=\mathrm{2}{c}_{\mathrm{11}}{\alpha }_{\mathrm{3}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}，r_{\mathrm{23}}=(c_{\mathrm{11}}(W_{\mathrm{11}}^{\mathrm{3}}(-\mathrm{1})+W_{\mathrm{20}}^{\mathrm{3}}(-\mathrm{1})+W_{\mathrm{11}}^{\mathrm{2}}(\mathrm{0}){\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\alpha }_{\mathrm{3}}+W_{\mathrm{20}}^{\mathrm{2}}(-\mathrm{1}){\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\overline{\alpha }}_{\mathrm{3}}，

r_{\mathrm{31}}=\frac{d_{\mathrm{11}}}{\mathrm{2}}+d_{\mathrm{12}}+d_{\mathrm{13}}{\alpha }_{\mathrm{3}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}，r_{\mathrm{32}}=\frac{d_{\mathrm{11}}}{\mathrm{2}}+\mathrm{2}{d}_{\mathrm{12}}+d_{\mathrm{13}}({\overline{\alpha }}_{\mathrm{3}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}+{\alpha }_{\mathrm{3}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}})，

r_{\mathrm{33}}=d_{\mathrm{12}}[W_{\mathrm{11}}^{\mathrm{3}}(\mathrm{0})+W_{\mathrm{20}}^{\mathrm{3}}(\mathrm{0})+W_{\mathrm{11}}^{\mathrm{1}}(\mathrm{0})+W_{\mathrm{20}}^{\mathrm{1}}(\mathrm{0})]\mathrm{。}

基于文献［14］方法，可得如下系数：

g_{\mathrm{20}}={\tau }_{\mathrm{0}}\overline{B}[r_{\mathrm{11}}+{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{21}}+{\overline{\beta }}_{\mathrm{3}}{r}_{\mathrm{31}}]，

g_{\mathrm{11}}={\tau }_{\mathrm{0}}\overline{B}(r_{\mathrm{12}}+{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{22}}+r_{\mathrm{32}}{\overline{\beta }}_{\mathrm{3}})，

g_{\mathrm{02}}={\tau }_{\mathrm{0}}\overline{B}[{\overline{r}}_{\mathrm{11}}+{\overline{\beta }}_{\mathrm{2}}{\overline{r}}_{\mathrm{21}}+r_{\mathrm{33}}{\overline{\beta }}_{\mathrm{3}}]，

\begin{array}{l}{g}_{\mathrm{21}}={\tau }_{\mathrm{0}}\overline{B}\{(b_{\mathrm{11}}+d_{\mathrm{12}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{20}}^{\mathrm{1}}(\mathrm{0})+(b_{\mathrm{11}}+d_{\mathrm{12}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{11}}^{\mathrm{1}}(\mathrm{0})+(c_{\mathrm{11}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\overline{\alpha }}_{\mathrm{3}}{\overline{\beta }}_{\mathrm{3}}+d_{\mathrm{12}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\overline{\alpha }}_{\mathrm{3}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{20}}^{\mathrm{2}}(\mathrm{0})+\\ (c_{\mathrm{11}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\alpha }_{\mathrm{3}}{\overline{\beta }}_{\mathrm{2}}+d_{\mathrm{13}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\alpha }_{\mathrm{3}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{11}}^{\mathrm{2}}(\mathrm{0})+(b_{\mathrm{11}}+\frac{d_{\mathrm{11}}}{\mathrm{2}}{\overline{\beta }}_{\mathrm{3}}+d_{\mathrm{12}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{11}}^{\mathrm{3}}(\mathrm{0})+\\ (b_{\mathrm{11}}+\frac{d_{\mathrm{11}}}{\mathrm{2}}{\overline{\beta }}_{\mathrm{3}}+d_{\mathrm{12}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{20}}^{\mathrm{3}}(\mathrm{0})+(c_{\mathrm{11}}{\overline{\beta }}_{\mathrm{2}}+d_{\mathrm{11}}{\overline{\beta }}_{\mathrm{3}})W_{\mathrm{11}}^{\mathrm{3}}(-\mathrm{1})+c_{\mathrm{11}}{\overline{\beta }}_{\mathrm{2}}{W}_{\mathrm{20}}^{\mathrm{3}}(-\mathrm{1})\mathrm{。}\end{array}

为了确定g_{\mathrm{21}}，还需计算出W_{\mathrm{20}}(\theta )和W_{\mathrm{11}}(\theta )，由式（11） （12）可得

\dot{\mathit{W}}=\left\{\begin{array}{l}\mathbf{A}\mathbf{W}-\mathrm{2}\mathrm{R}\mathrm{e}\left\{{{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})\mathit{F}(\mu ,\mathit{\phi })\mathit{q}(\theta )\right\},\theta \in [-\mathrm{1,0}),\\ \mathbf{A}\mathbf{W}-\mathrm{2}\mathrm{R}\mathrm{e}\left\{{{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})\mathit{F}(\mu ,\mathit{\phi })\mathit{q}(\theta )\right\}+\mathit{F}((\mu ,\mathit{\phi }),\theta =\mathrm{0}\mathrm{。}\end{array}\right.

即

\dot{W}=\mathbf{A}\mathbf{W}+\mathit{H}(z,\overline{z},\theta )，（15）

其中

\mathit{H}(z,\overline{z},\theta )={\mathit{H}}_{\mathrm{20}}(\theta )\frac{z^{\mathrm{2}}}{\mathrm{2}}+{\mathit{H}}_{\mathrm{11}}(\theta )z\overline{z}+{\mathit{H}}_{\mathrm{02}}(\theta )\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+\cdots。（16）

组合式（12） （13），再将式（17）代入式（16）可得

\left\{\begin{array}{l}(\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}-\mathit{A}){\mathit{W}}_{\mathrm{20}}(\theta )={\mathit{H}}_{\mathrm{20}}(\theta ),\\ -\mathit{A}{\mathit{W}}_{\mathrm{11}}(\theta )={\mathit{H}}_{\mathrm{11}}(\theta )\mathrm{。}\end{array}\right.

由式（16）可知，当\theta \in [-\mathrm{1,0})时，

\begin{array}{l}\mathit{H}(z,\overline{z},\theta )=-\mathrm{2}\mathrm{R}\mathrm{e}\left\{{{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})F(\mu ,\phi )\mathit{q}(\theta )\right\}=\\ -{{\overline{\mathit{q}}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0})F_{\mathrm{0}}\mathit{q}(\theta )-{{\mathit{q}}^{\mathrm{*}}}^{\mathrm{T}}(\mathrm{0}){\overline{F}}_{\mathrm{0}}\overline{\mathit{q}}(\theta )=\\ -g(z,\overline{z})\mathit{q}(\theta )-\overline{g}(z,\overline{z})\overline{\mathit{q}}(\theta )=\\ -\left(g_{\mathrm{20}}\frac{z^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{11}}z\overline{z}+g_{\mathrm{02}}\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{21}}\frac{z^{\mathrm{2}}\overline{z}}{\mathrm{2}}\right)\mathit{q}(\theta )-\left({\overline{g}}_{\mathrm{20}}\frac{{\overline{z}}^{\mathrm{2}}}{\mathrm{2}}+{\overline{g}}_{\mathrm{11}}z\overline{z}+{\overline{g}}_{\mathrm{02}}\frac{z^{\mathrm{2}}}{\mathrm{2}}+g_{\mathrm{21}}\frac{{\overline{z}}^{\mathrm{2}}z}{\mathrm{2}}\right)\overline{\mathit{q}}(\theta )\mathrm{。}\end{array}（17）

通过比较式（16） （17）的系数，有

\left\{\begin{array}{l}{\mathit{H}}_{\mathrm{20}}(\theta )=-g_{\mathrm{20}}\mathit{q}(\theta )-{\overline{g}}_{\mathrm{20}}\overline{\mathit{q}}(\theta ),\\ {\mathit{H}}_{\mathrm{11}}(\theta )=-g_{\mathrm{11}}\mathit{q}(\theta )-{\overline{g}}_{\mathrm{11}}\overline{\mathit{q}}(\theta )\mathrm{。}\end{array}\right.

由\mathit{A}(\mu )\mathit{\phi }(\theta )的定义以及(\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}-\mathit{A}){\mathit{W}}_{\mathrm{20}}(\theta )={\mathit{H}}_{\mathrm{20}}(\theta )可得

\begin{array}{l}{\dot{\mathit{W}}}_{\mathrm{20}}(\theta )=\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}{\mathit{W}}_{\mathrm{20}}(\theta )-{\mathit{H}}_{\mathrm{20}}(\theta )=\\ \mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}{\mathit{W}}_{\mathrm{20}}(\theta )+g_{\mathrm{20}}\mathit{q}(\theta )+{\overline{g}}_{\mathrm{20}}\overline{\mathit{q}}(\theta )=\\ \mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}{\mathit{W}}_{\mathrm{20}}(\theta )+g_{\mathrm{20}}\mathit{q}(\mathrm{0}){\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }+{\overline{g}}_{\mathrm{20}}\overline{\mathit{q}}(\mathrm{0}){\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }\mathrm{。}\end{array}

同理可算出{\dot{W}}_{\mathrm{11}}(\theta )，其中{\mathit{W}}_{\mathrm{20}}(\theta )，{\mathit{W}}_{\mathrm{11}}(\theta )结果如下：

{\mathit{W}}_{\mathrm{20}}(\theta )=\frac{\mathrm{i}{g}_{\mathrm{20}}\mathit{q}(\mathrm{0})}{{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }+\frac{\mathrm{i}{\overline{g}}_{\mathrm{20}}\overline{\mathit{q}}(\mathrm{0})}{\mathrm{3}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }+{\mathit{E}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }，

{\mathit{W}}_{\mathrm{11}}(\theta )=-\frac{\mathrm{i}{g}_{\mathrm{11}}\mathit{q}(\mathrm{0})}{{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }+\frac{\mathrm{i}{\overline{g}}_{\mathrm{11}}\overline{\mathit{q}}(\mathrm{0})}{{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}\theta }+{\mathit{E}}_{\mathrm{2}}。

这里，{\mathit{E}}_{\mathrm{1}}=\left(E_{\mathrm{1}}^{(\mathrm{1})},E_{\mathrm{1}}^{(\mathrm{2})},E_{\mathrm{1}}^{(\mathrm{3})},E_{\mathrm{1}}^{(\mathrm{4})}\right),{\mathit{E}}_{\mathrm{2}}=\left(E_{\mathrm{2}}^{(\mathrm{1})},E_{\mathrm{2}}^{(\mathrm{2})},E_{\mathrm{2}}^{(\mathrm{3})},E_{\mathrm{2}}^{(\mathrm{4})}\right)\in {\mathbb{R}}^{\mathrm{4}}是常量，在\theta =\mathrm{0}时，根据\mathit{A}(\mu )\mathit{\phi }(\theta )的定义得到{\mathit{E}}_{\mathrm{1}}和{\mathit{E}}_{\mathrm{2}}的解为

{\mathit{E}}_{\mathrm{1}}={\left(\begin{array}{cccc}\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{11}}& -a_{\mathrm{12}}& -a_{\mathrm{13}}& \mathrm{0}\\ -a_{\mathrm{21}}& \mathrm{2}i{\omega }_{\mathrm{0}}-a_{\mathrm{22}}& -b_{\mathrm{23}}{\mathrm{e}}^{-\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}& \mathrm{0}\\ -a_{\mathrm{31}}& -a_{\mathrm{32}}& \mathrm{2}i{\omega }_{\mathrm{0}}-a_{\mathrm{33}}-b_{\mathrm{33}}{\mathrm{e}}^{-\mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -a_{\mathrm{43}}& \mathrm{2}\mathrm{i}{\omega }_{\mathrm{0}}-a_{\mathrm{44}}\end{array}\right)}^{-\mathrm{1}}\left(\begin{array}{c}{r}_{\mathrm{11}}\\ \begin{array}{l}{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{21}}\\ {\overline{\beta }}_{\mathrm{3}}{r}_{\mathrm{31}}\\ \mathrm{0}\end{array}\end{array}\right)，

{\mathit{E}}_{\mathrm{2}}=-{\left(\begin{array}{cccc}{a}_{\mathrm{11}}& a_{\mathrm{12}}& a_{\mathrm{13}}& \mathrm{0}\\ a_{\mathrm{21}}& a_{\mathrm{22}}& b_{\mathrm{23}}& \mathrm{0}\\ a_{\mathrm{31}}& a_{\mathrm{32}}& a_{\mathrm{33}}+b_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& a_{\mathrm{43}}& a_{\mathrm{44}}\end{array}\right)}^{-\mathrm{1}}\left(\begin{array}{c}{r}_{\mathrm{12}}\\ \begin{array}{l}{\overline{\beta }}_{\mathrm{2}}{r}_{\mathrm{22}}\\ {\overline{\beta }}_{\mathrm{3}}{r}_{\mathrm{32}}\\ \mathrm{0}\end{array}\end{array}\right)。

于是有

C_{\mathrm{1}}(\mathrm{0})=\frac{\mathrm{i}}{\mathrm{2}{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}(g_{\mathrm{20}}{g}_{\mathrm{11}}-\mathrm{2}{\left|g_{\mathrm{11}}\right|}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{3}}{\left|g_{\mathrm{02}}\right|}^{\mathrm{2}})+\frac{\mathrm{1}}{\mathrm{2}}{g}_{\mathrm{21}}，（18）

{\mu }_{\mathrm{2}}=\frac{\mathrm{R}\mathrm{e}\left\{C_{\mathrm{1}}(\mathrm{0})\right\}}{\mathrm{R}\mathrm{e}\left\{{\lambda }^{\text{'}}({\tau }_{\mathrm{0}})\right\}},{\beta }_{\mathrm{2}}=\mathrm{2}\mathrm{R}\mathrm{e}{C}_{\mathrm{1}}(\mathrm{0}),T_{\mathrm{2}}=\frac{\mathrm{I}\mathrm{m}\left\{C_{\mathrm{1}}(\mathrm{0})\right\}+{\mu }_{\mathrm{2}}\mathrm{I}\mathrm{m}\left\{{\lambda }^{\text{'}}({\tau }_{\mathrm{0}})\right\}}{{\omega }_{\mathrm{0}}{\tau }_{\mathrm{0}}}\mathrm{。}（19）

定理2 对于系统（2），由\tau ={\tau }_{\mathrm{0}}得到的Hopf分支的周期解可由式（19）描述，在\mu =\mathrm{0}邻域中，有

（1）{\mu }_{\mathrm{2}}确定Hopf分支的方向，如果，则Hopf分支是超临界的（次临界的）；

（2）{\beta }_{\mathrm{2}}确定分支周期解的稳定性，如果，则分支周期解是稳定的（不稳定的）；

（3）T_{\mathrm{2}}确定分支周期解的周期，如果，则分支周期解的周期增大（减小）。

3 数值模拟

对系统（2）进行数值模拟，给出一个局部Hopf分支存在性的数值结果。考虑系统

\left\{\begin{array}{l}\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\mathrm{0.001}\times (\mathrm{1}-\mathrm{0.7})\times \mathrm{10}-(\mathrm{0.001}+\mathrm{0.05})S(t)-\frac{\mathrm{0.9}S(t)I(t)}{\mathrm{1}+\mathrm{0.9}I(t)}+\mathrm{0.01}V(t),\\ \frac{\mathrm{d}V(t)}{\mathrm{d}t}=\mathrm{0.01}\times \mathrm{0.7}+\mathrm{0.05}S(t)-\mathrm{0.4}\times \mathrm{10}V(t)I(t)-(\mathrm{0.001}+\mathrm{0.01})V(t),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}=\frac{\mathrm{0.9}S(t)I(t)}{\mathrm{1}+\mathrm{0.9}I(t)}+\mathrm{0.4}\times \mathrm{10}V(t)I(t)-(\mathrm{0.001}+\mathrm{0.2}+\mathrm{0.01})I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}=\mathrm{0.2}I(t)-\mathrm{0.001}R(t)\mathrm{。}\end{array}\right.（20）

通过计算可得知系统（20）存在正平衡点E^{\mathrm{*}}=(\mathrm{0.036,0.044,0.057,9.401})，且{\omega }_{\mathrm{0}}\approx \mathrm{0.082},{\tau }_{\mathrm{0}}\approx \mathrm{13.979}。

系统（20）在时滞参数为和\tau >{\tau }_{\mathrm{0}}时的数值模拟结果对比见图1和图2，当时，正平衡点E^{\mathrm{*}}是渐近稳定的；当\tau =\mathrm{15}>{\tau }_{\mathrm{0}}\approx \mathrm{13.979}时，系统（20）的正平衡点失稳，产生了一个小振幅的极限环（周期解），并在正平衡点附近振动。

4 结论

本研究针对潜伏期时滞效应下的SVIR传染病动力学模型展开系统性分析，揭示了时滞参数对疾病传播稳定性的关键影响机制。通过构建具有潜伏期特征的时滞微分方程，结合Hopf分支理论与 中心流形定理，从理论上证明了系统在正平衡点处存在局部渐近稳定性阈值。再结合数值仿真结果，发现与理论推理证明相同。表明当潜伏期时滞低于临界值{\tau }_{\mathrm{0}}时，系统的稳定性特征可支持现有防控措施（如症状监测、病例隔离及密接追踪）的有效实施，这为早期疫情控制提供了重要依据。当潜伏期超过临界阈值时，系统的Hopf分支现象将导致周期性振荡产生，表明疾病可能呈现波动式传播态势。结合病毒传播的真实环境，这种情况下病毒传播会难以控制，极大危害公众健康。

参考文献

[1] ESSOUFI E H， ZAFRAR A. Boundary optimal control of time-space SIR model with nonlinear Robin boundary condition［J］. International Journal of Dynamics and Control， 2022， 10（4）： 1279-1290.

[2] KERMACK W O， MCKENDRICK À. A contribution to the mathematical theory of epidemics［J］. Proc. A 1 August 1927， 115 （772）： 700-721.

[3] SUN G W， JIN Z， MAI A L. Dynamics of a two-patch SIR model with disease surveillance mediated infection force［J］. Communications in Nonlinear Science and Numerical Simulation， 2024， 132： 107872.

[4] DAS S， SRIVASTAVA P K， BISWAS P. Exploring Hopf-bifurcations and endemic bubbles in a tuberculosis model with behavioral changes and treatment saturation［J］. Chaos， 2024， 34（1）： 013126.

[5] GAO D Z， LEI C X， PENG R， et al. A diffusive SIS epidemic model with saturated incidence function in a heterogeneous environment［J］. Nonlinearity， 2024， 37（2）： 025002.

[6] CAPASSO V， SERIO G. A generalization of the Kermack-McKendrick deterministic epidemic model［J］. Mathematical Biosciences， 1978， 42（1-2）： 43-61.

[7] SHI R Q， ZHANG Y H. Stability analysis and Hopf bifurcation of a fractional order HIV model with saturated incidence rate and time delay［J］. Alexandria Engineering Journal， 2024， 108： 70-88.

[8] MAO S X， XU R， LI Z， et al. Global stability of an eco-epidemiological model with time delay and saturation incidence［J］. Discrete Dynamics in Nature and Society， 2011， 2011： 730783.

[9] LIU X S， DAI B X. Qualitative and bifurcation analysis of an SIR epidemic model with saturated treatment function and nonlinear pulse vaccination［J］. Discrete Dynamics in Nature and Society， 2016， 2016： 9146481.

[10] JI X， YUAN S， LI J. Stability of a stochastic seis model with saturation incidence and latent period［J］. Journal of Applied Analysis and Computation， 2017， 7（4）： 1652-1673.

[11] EZEKIEL I D， IYASE S A， ANAKE T. Stability and Hopf bifurcation analysis of an infectious disease delay model［J］. Mathematics and Computers in Simulation， 2023， 213： 78-114.

[12] ATTAULLAH A， KHURSHAID A， ZEESHAN Z， et al. Computational framework of the SVIR epidemic model with a non-linear saturation incidence rate［J］. Axioms， 2022， 11（11）： 651.

[13] 汪雨琴， 谢景力， 赵佳璐， 等. 一类具有时滞效应的SIS传染病模型［J］. 吉首大学学报（自然科学版）， 2025， 46（1）： 12-17.

[14] 豆中丽. 阶段结构的分数阶时滞生态流行病模型的动力学分析［J］. 云南大学学报（自然科学版）， 2025， 47（1）： 10-20.

[15] HASSARD B D， KAZARINOFF N D， WAN Y H. Theory and applications of Hopf bifurcation［M］. Cambridge： Cambridge University Press， 1981.

基金资助

吉林省教育厅科学技术研究资助项目“乙类乙管时期新型冠状病毒肺炎传播的建模分析与风险预测”(JJKH20240891KJ)