一类具有非线性发生率的随机SIR模型的Melnikov混沌

Melnikov Analysis of Chaos in a Stochastic SIR Model with Nonlinear Incidence Rate

摘要

针对疾病的流行往往受到随机因素的干扰，本文通过对传染率施加有界噪声扰动，研究一类具有非线性发生率的SIR （Susceptible-Infectious-Recovered，易感-感染-康复） 模型的混沌动力学特性。基于同宿分岔理论，结合随机Melnikov方法，从数学理论的角度推导出模型可能发生混沌的充分条件，得到模型在有界噪声扰动下发生混沌的阈值，阐明噪声强度对SIR模型混沌动力学行为的影响。同时，利用数值模拟验证理论分析结果，数值模拟包括势能图、阈值曲线图、相图。研究成果有助于预测疾病传播态势，为公共卫生防控策略制定提供理论支撑，助力疫情科学防控。

Abstract

In view of the fact that the prevalence of diseases is often disturbed by random factors, this paper studies the chaotic dynamics of a class of SIR model with nonlinear incidence rate through the infection rate disturbed by bounded noise. Based on the homoclinic bifurcation and by using the stochastic Melnikov theory, from the perspective of mathematical theory, the sufficient condition for the possible occurrence of chaos in the model is derived, and the threshold of chaos in the model under bounded noise disturbance is obtained. Moreover, the influence of the noise amplitude on the chaotic dynamic behavior of the SIR model is elucidated. Numerical simulations, including potential, threshold curve graphs, and phase portraits, are used to verify the results of the theoretical analysis. The research results contribute to accurately predicting the disease transmission trends, provide theoretical support for formulating public health prevention and control strategies, and facilitate the scientific prevention and control of epidemics.

Graphical abstract

关键词

非线性发生率 / Melnikov理论 / 混沌 / 同宿分岔

Key words

nonlinear incidence rate / Melnikov theorem / chaos / homoclinic bifurcation

0　引言

传染病的流行严重危害人类健康，是影响社会安定的一个严峻问题。为有效防治传染病，众多学者^［1-3］通过建立数学模型来研究疾病的传播规律，从而预测疾病的流行趋势，为疾病预防与控制提供科学理论依据和应对措施。

本文考虑一类具有非线性发生率的SIR（Susceptible-Infectious-Recovered，易感-感染-康复）模型，

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=\mu (\mathrm{1}-S)-\beta (t)I^{\mathrm{2}}S,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\beta (t)I^{\mathrm{2}}S-(\gamma +\mu )I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}=\gamma I-\mu R,\end{array}\right.（1）

其中S(t),I(t),R(t)分别表示易感者，染病者和恢复者。\mu表示人口自然死亡率，\gamma是染病者的恢复率，两者均为正常数。假设人口数量恒定：S(t)+I(t)+R(t)=\mathrm{1}。非线性发生率\beta (t)I^{\mathrm{2}}S中\beta (t)表示染病者的传染率，这里\beta (t)>\mathrm{0}。

疾病的发病率表示单位时间内的新增病例数，研究具有非线性发生率的传染病模型，在流行病学的研究中起着十分重要的作用。很多学者^［4-8］对具有非线性发生率的传染病模型做了广泛的研究。如，Ruan 和Wang^［5］对一类具有生命动力学和饱和质量作用非线性发生率的传染病模型进行了全局定性分析和分岔分析；Jin、Wang和Xiao^［8］研究了一类具有非线性发生率的SIRS（Susceptible-Infectious-Recovered-Susceptible，易感-感染-康复-易感）模型的全局动力学行为。在疾病的传播过程中，由于各种形式的随机干扰无处不在，因此，为了更好地接近生物系统的实际情况，研究人员在模型中加入了噪声扰动^［9-12］，建立随机传染病模型，进而研究模型的动力学行为。如，Jiang等^［9］研究了一类具有双线性发生率的随机SIR模型的渐近性行为；Lahrouz和Omari^［10］研究了一类具有一般发生率的随机SIRS模型的灭绝性和平稳分布；曹虹等^［12］研究了一类具有非线性发生率的随机SIRS模型的全局正解唯一性和渐近行为。二十世纪八九十年代，相关数值研究结论^［13-14］表明：当周期激励传染率中的振幅充分大时，流行病模型可能存在混沌行为，这引起了研究人员的极大兴趣并促使许多研究人员在流行病学时间序列里寻找混沌的存在^［15-16］。Glendinning和Perry^［17］研究了周期激励传染率的SIR模型，随着周期振幅的变化，观察到了从倍周期分岔到混沌轨迹的产生，文中采用的非线性发生率为\beta (t)I^{\mathrm{2}}S，传染率为\beta (t)=\mathrm{16}\mu \left(\mathrm{1}+{\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}+{\epsilon }^{\mathrm{4}}{b}_{\mathrm{2}}+{\epsilon }^{\mathrm{5}}{b}_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}(\mu \epsilon \Omega t)\right)。 Diallo和Koné^［18］研究了同文献［17］的模型，从数学上证明了模型中混沌运动的存在，文中采用的非线性发生率为\beta (t)I^{p}S，传染率为\beta (t)=\frac{p^{\mathrm{2}p}}{{(p-\mathrm{1})}^{\mathrm{2}p-\mathrm{1}}}\mu \left(\mathrm{1}+{\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}+{\epsilon }^{\mathrm{4}}{b}_{\mathrm{2}}+{\epsilon }^{\mathrm{5}}{b}_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}(\mu \epsilon \Omega t)\right)。Diallo等^［19］研究了同文献［17］的模型的混沌，其非线性发生率同文献［18］，传染率为如下形式：

\beta (t)=\frac{p^{\mathrm{2}p}}{{(p-\mathrm{1})}^{\mathrm{2}p-\mathrm{1}}}\mu \left(\mathrm{1}+{\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}+{\epsilon }^{\mathrm{4}}{b}_{\mathrm{2}}+{\epsilon }^{\mathrm{5}}{b}_{\mathrm{3}}\left(\mathrm{s}\mathrm{i}\mathrm{n}(\mu \epsilon \Omega t)+\mathrm{s}\mathrm{i}\mathrm{n}(\sqrt[]{\mathrm{2}}\mu \epsilon \Omega t)\right)\right)。

在本文中，我们对模型（1）的传染率进行有界噪声干扰，研究模型（1）的混沌动力学行为。传染率为

\beta (t)=\mathrm{16}\mu \left(\mathrm{1}+{\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}+{\epsilon }^{\mathrm{4}}{b}_{\mathrm{2}}+{\epsilon }^{\mathrm{5}}{b}_{\mathrm{3}}\xi (t)\right)，（2）

其中b_{\mathrm{1}}\mathrm{、}{b}_{\mathrm{2}}\mathrm{、}{b}_{\mathrm{3}}均是一阶常数，满足：\mathrm{4}{b}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}}>\mathrm{0}\mathrm{。}\xi (t)属于有界噪声，其数学表达式为

\xi (t)=b_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{\mathrm{1}}t+\Psi ), \Psi =\sigma B(t)+\Gamma，

其中b_{\mathrm{3}}表示有界噪声的振幅，{\Omega }_{\mathrm{1}}和\sigma是正常数，{\Omega }_{\mathrm{1}}为有界噪声的平均频率，\sigma为噪声强度，B(t)为标准Wiener过程，\Gamma是在[\mathrm{0,2}\mathrm{\pi })内均匀分布的随机变量，表示随机相角。\xi (t)为广义上的宽带稳定有界随机过程^［20］，均值为零的非Gauss分布，其协方差函数为

C_{\xi }(\tau )=\frac{{b_{\mathrm{3}}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\sigma }^{\mathrm{2}}\left|\tau \right|}{\mathrm{2}}\right)\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{\mathrm{1}}\tau )\mathrm{。}

方差为C_{\xi }(\mathrm{0})={b_{\mathrm{3}}}^{\mathrm{2}}/\mathrm{2}，密度函数为p(\xi )=\mathrm{1}/\left(\mathrm{\pi }\sqrt[]{{b_{\mathrm{3}}}^{\mathrm{2}}-{\xi }^{\mathrm{2}}}\right)，谱密度函数为

S_{\xi }(\omega )=\frac{(b_{\mathrm{3}}{\sigma )}^{\mathrm{2}}}{\mathrm{2}\mathrm{\pi }}\left(\frac{\mathrm{1}}{\mathrm{4}(\omega -{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}+\frac{\mathrm{1}}{\mathrm{4}(\omega +{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}\right)=\frac{{b_{\mathrm{3}}}^{\mathrm{2}}{\sigma }^{\mathrm{2}}}{\mathrm{4}\mathrm{\pi }}\frac{{\omega }^{\mathrm{2}}+{{\Omega }_{\mathrm{1}}}^{\mathrm{2}}+\frac{{\sigma }^{\mathrm{4}}}{\mathrm{4}}}{\left({\omega }^{\mathrm{2}}-{{\Omega }_{\mathrm{1}}}^{\mathrm{2}}-\frac{{\sigma }^{\mathrm{4}}}{\mathrm{4}}\right)+{\sigma }^{\mathrm{4}}{\omega }^{\mathrm{2}}}\mathrm{。}（3）

依据（3）式，图1画出了几组不同参数下有界噪声的谱密度。

不同于已有文献主要关注确定性条件下的模型，本文着眼于随机环境中模型的构建及混沌动力学分析，从数学理论的角度推导出模型可能发生混沌的充分条件，得到模型在有界噪声扰动下发生混沌的阈值。

1　随机Melnikov过程

考虑单自由度的Hamiltonian系统，该系统受到有界噪声的干扰，系统的方程形式如下，

\left\{\begin{array}{l}\dot{x}=\frac{\partial H}{\partial y},\\ \dot{y}=-\frac{\partial H}{\partial x}-\epsilon g(x, y)\frac{\partial H}{\partial y}+\epsilon f(x, y)\xi (t),\end{array}\right.（4）

其中x和y分别表示广义的位移和动量；H=H(x, y)是具有连续一阶导数的Hamiltonian量，\epsilon是任意小的正数，g(x, y)表示准线性阻尼系数，f(x, y)表示激励振幅。假设系统（4）有一个通过同宿轨道(x_{\mathrm{0}}(t), y_{\mathrm{0}}(t))连接到自身的双曲不动点。

如文献［21］所述，系统（4）的随机Melnikov过程为

M(t_{\mathrm{0}})={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\partial H}{\partial y}\left[-g(x_{\mathrm{0}}, y_{\mathrm{0}})\frac{\partial H}{\partial y}+f(x_{\mathrm{0}}, y_{\mathrm{0}})\xi (t+t_{\mathrm{0}})\right]\mathrm{d}t=M_d+Z(t_{\mathrm{0}}),（5）

其中M_d是Melnikov过程中因阻尼引起的分量，Z(t_{\mathrm{0}})是有界噪声干扰引起的分量。Melnikov过程的均值为

E[M(t_{\mathrm{0}})]=M_d=-{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}g(x_{\mathrm{0}}, y_{\mathrm{0}}){\left(\frac{\partial H}{\partial y}\right)}^{\mathrm{2}}\mathrm{d}t,（6）

其中E[•]表示期望算子。对于正阻尼，（6）式总是得出负值，这意味着系统（4）不可能有均值意义下的混沌。随机Melnikov过程（5）式中两个分量的均方值分别为

{M_d}^{\mathrm{2}}={\left({\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}g(x_{\mathrm{0}}, y_{\mathrm{0}}){\left(\frac{\partial H}{\partial y}\right)}^{\mathrm{2}}\mathrm{d}t\right)}^{\mathrm{2}},

{{\sigma }_Z}^{\mathrm{2}}=E[Z^{\mathrm{2}}(t_{\mathrm{0}})]=E\left[{\left({\int }_{-\infty }^{+\infty }\frac{\partial H}{\partial y}f(x_{\mathrm{0}}, y_{\mathrm{0}})\xi (t+t_{\mathrm{0}})\mathrm{d}t\right)}^{\mathrm{2}}\right]\mathrm{。}

为计算{{\sigma }_Z}^{\mathrm{2}}，（5）式中的Z(t_{\mathrm{0}})被改写为卷积积分，即

Z(t_{\mathrm{0}})={\int }_{-\infty }^{+\infty }(\partial H/\partial y)f(x_{\mathrm{0}}, y_{\mathrm{0}})\xi (t+t_{\mathrm{0}})\mathrm{d}t=h(t)\mathrm{*}\xi (t),

当\xi (t)是系统的输入时，h(t)=(\partial H/\partial y)f(x_{\mathrm{0}}, y_{\mathrm{0}})可视作时不变线性系统的脉冲响应函数。因此，可在频域中得到系统的方差，如下：

{{\sigma }_Z}^{\mathrm{2}}={{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left|H(\omega )\right|}^{\mathrm{2}}{S}_{\xi }(\omega )\mathrm{d}\omega ,

其中H(\omega )是系统的频率响应函数，它是脉冲响应函数h(t)的傅里叶变换。

综上，当{M_d}^{\mathrm{2}}={{\sigma }_Z}^{\mathrm{2}}时，随机Melnikov过程（5）在均方意义下存在简单零点。随机Melnikov过程通常应用于非线性振子系统的混沌检测^［22-23］。

2　模型的尺度变换

为了利用随机Melnikov理论对模型（1）的混沌进行分析，我们将对模型（1）进行一系列坐标变换，将模型化成扰动的Hamiltonian系统。由于S(t)=\mathrm{1}-I(t)-R(t)，模型（1）可以写成

\left\{\begin{array}{l}\frac{\mathrm{d}I}{\mathrm{d}t}=\beta (t)I^{\mathrm{2}}(\mathrm{1}-I-R)-(\gamma +\mu )I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}=\gamma I-\mu R\mathrm{。}\end{array}\right.

定义一个新的时间尺度t^{\text{'}}=\mu t， 那么

\left\{\begin{array}{l}\frac{\mathrm{d}I}{\mathrm{d}{t}^{\text{'}}}={\beta }^{\text{'}}(t^{\text{'}})I^{\mathrm{2}}(\mathrm{1}-I-R)-({\gamma }^{\text{'}}+\mathrm{1})I,\\ \frac{\mathrm{d}R}{\mathrm{d}{t}^{\text{'}}}={\gamma }^{\text{'}}I-R,\end{array}\right.（7）

这里，\beta =\mu {\beta }^{\text{'}}，\gamma =\mu {\gamma }^{\text{'}}。定义新变量z=\gamma I-R，因此，系统（7）变为

\left\{\begin{array}{l}\dot{R}=z=P_{\mathrm{1}}(R,z),\\ \dot{z}=\gamma \dot{I}-\dot{R}=P_{\mathrm{2}}(R,z)\mathrm{。}\end{array}\right.（8）

将系统（7）代入系统（8）的第二个方程，并结合I={\gamma }^{-\mathrm{1}}(R+z)， 则系统（8）变为

\left\{\begin{array}{l}\dot{R}=z=P_{\mathrm{1}}(R,z),\\ \begin{array}{l}\dot{z}=-(\mathrm{1}+\gamma )R+\beta (t)\frac{\mathrm{1}}{\gamma }{R}^{\mathrm{2}}-\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{1}+\frac{\mathrm{1}}{\gamma }\right)R^{\mathrm{3}}-(\mathrm{2}+\gamma )z+\beta (t)\frac{\mathrm{1}}{\gamma }{z}^{\mathrm{2}}-\beta (t)\frac{\mathrm{1}}{{\gamma }^{\mathrm{2}}}{z}^{\mathrm{3}}+\\ \mathrm{2}\beta (t)\frac{\mathrm{1}}{\gamma }Rz-\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{1}+\frac{\mathrm{3}}{\gamma }\right)R{z}^{\mathrm{2}}-\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{2}+\frac{\mathrm{3}}{\gamma }\right)R^{\mathrm{2}}z=P_{\mathrm{2}}(R,z)\mathrm{。}\end{array}\end{array}\right.（9）

系统（9）的驻点有z=\mathrm{0}且R=\mathrm{0}（对应无病平衡点），或z=\mathrm{0}且R满足：

\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{1}+\frac{\mathrm{1}}{\gamma }\right)R^{\mathrm{2}}-\beta (t)\frac{\mathrm{1}}{\gamma }R+(\mathrm{1}+\gamma )=\mathrm{0}, R\ne \mathrm{0}。（10）

若\mathrm{\Delta }>\mathrm{0}， （10）式有两个实根，即{\left(\beta (t)\mathrm{1}/\gamma \right)}^{\mathrm{2}}-\mathrm{4}\beta (t)(\mathrm{1}/\gamma )\left(\mathrm{1}+(\mathrm{1}/\gamma )\right)(\mathrm{1}+\gamma )>\mathrm{0}，化简得\beta (t)>\mathrm{4}{(\mathrm{1}+\gamma )}^{\mathrm{2}}。因此，如果\beta (t)={\beta }_{\mathrm{s}\mathrm{n}}(t)=\mathrm{4}{(\mathrm{1}+\gamma )}^{\mathrm{2}}，系统（9）有一个鞍结点分岔，此时地方病平衡点为

S_{\mathrm{s}\mathrm{n}}=\mathrm{1}/\mathrm{2}, I_{\mathrm{s}\mathrm{n}}=\mathrm{1}/\mathrm{2}(\mathrm{1}+\gamma ), R_{\mathrm{s}\mathrm{n}}=\gamma /\mathrm{2}(\mathrm{1}+\gamma )\mathrm{。}

系统（9）在z=\mathrm{0}处的Jacobian矩阵为\left(\begin{array}{cc}\partial P_{\mathrm{1}}/\partial R& \partial P_{\mathrm{1}}/\partial z\\ \partial P_{\mathrm{2}}/\partial R& \partial P_{\mathrm{2}}/\partial z\end{array}\right)=\left(\begin{array}{cc}\mathrm{0}& \mathrm{1}\\ J_{\mathrm{21}}& J_{\mathrm{22}}\end{array}\right)， 这里，

J_{\mathrm{21}}=-(\mathrm{1}+\gamma )+\mathrm{2}\beta (t)\frac{\mathrm{1}}{\gamma }R-\mathrm{3}\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{1}+\frac{\mathrm{1}}{\gamma }\right)R^{\mathrm{2}}，

J_{\mathrm{22}}=-(\mathrm{2}+\gamma )+\mathrm{2}\beta (t)\frac{\mathrm{1}}{\gamma }R-\beta (t)\frac{\mathrm{1}}{\gamma }\left(\mathrm{2}+\frac{\mathrm{3}}{\gamma }\right)R^{\mathrm{2}}。

当\beta (t)=\mathrm{4}{(\mathrm{1}+\gamma )}^{\mathrm{2}}时，矩阵行列式-J_{\mathrm{21}}=\mathrm{0}。如果迹J_{\mathrm{22}}=\mathrm{0}，我们会得到一种更退化的情形：Takens-Bogdanov点。如果要利用Melnikov方法，那么在该点附近可将系统按要求改写成近Hamiltonian形式^［24］。因此，若\gamma =\mathrm{1}, \beta =\mathrm{16}时，在鞍结点处Jacobian矩阵的迹等于零（此时\beta =\mathrm{4}{(\mathrm{1}+\gamma )}^{\mathrm{2}}，R=\gamma /\mathrm{2}(\mathrm{1}+\gamma )），此时退化平衡点为

S_{\mathrm{t}\mathrm{b}}=\mathrm{1}/\mathrm{2}, I_{\mathrm{t}\mathrm{b}}=\mathrm{1}/\mathrm{4}, R_{\mathrm{t}\mathrm{b}}=\mathrm{1}/\mathrm{4}。

为扰动该平衡点，通过以下方式定义新的变量和参数：

R=\frac{\mathrm{1}}{\mathrm{4}}+\theta , \gamma =\mathrm{1}+\vartheta , \beta (t)=\mathrm{16}\left(\mathrm{1}+b(t)\right),

其中\theta , \vartheta和b是任意小的正数。将上面变量代入系统（9），得到

\left\{\begin{array}{l}\dot{\theta }=z,\\ \dot{z}=\left(\frac{\mathrm{1}}{\mathrm{2}}b(t)-\frac{\mathrm{1}}{\mathrm{2}}\vartheta -\frac{\mathrm{1}}{\mathrm{4}}b(t)\vartheta \right)+\mathrm{2}b(t)\theta -\mathrm{8}{\theta }^{\mathrm{2}}+\left(\mathrm{3}b(t)-\vartheta \right)z-\mathrm{8}\theta z+O\left({\left(\left|\theta \right|+\left|z\right|+\left|\vartheta \right|+\left|b(t)\right|\right)}^{\mathrm{3}}\right)\mathrm{。}\end{array}\right.

设b(t)={\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}+{\epsilon }^{\mathrm{4}}{b}_{\mathrm{2}}+{\epsilon }^{\mathrm{5}}{b}_{\mathrm{3}}\xi (t), \vartheta ={\epsilon }^{\mathrm{2}}{b}_{\mathrm{1}}， 结合（2）式，得

\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}}b(t)-\frac{\mathrm{1}}{\mathrm{2}}\vartheta -\frac{\mathrm{1}}{\mathrm{4}}b(t)\vartheta =\left(\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{4}}{b_{\mathrm{1}}}^{\mathrm{2}}\right){\epsilon }^{\mathrm{4}}+\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{3}}{\epsilon }^{\mathrm{5}}\xi (t)+o({\epsilon }^{\mathrm{6}}),\\ \mathrm{2}b(t)=\mathrm{2}{b}_{\mathrm{1}}{\epsilon }^{\mathrm{2}}+o({\epsilon }^{\mathrm{4}}),\\ \mathrm{3}b(t)-\vartheta =\mathrm{2}{b}_{\mathrm{1}}{\epsilon }^{\mathrm{2}}+o({\epsilon }^{\mathrm{4}})\mathrm{。}\end{array}

定义新的坐标和时间：\theta ={\epsilon }^{\mathrm{2}}u, z={\epsilon }^{\mathrm{3}}v, \tau =\epsilon t,得到

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}\tau }=v,\\ \begin{array}{l}\frac{\mathrm{d}v}{\mathrm{d}\tau }=\left(\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{4}}{b_{\mathrm{1}}}^{\mathrm{2}}\right)+\mathrm{2}{b}_{\mathrm{1}}u-\mathrm{8}{u}^{\mathrm{2}}+\epsilon \left(\mathrm{2}{b}_{\mathrm{1}}v-\mathrm{8}uv+\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{3}}\xi (\tau )\right)+o({\epsilon }^{\mathrm{2}})=\\ \left(\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{8}}{b_{\mathrm{1}}}^{\mathrm{2}}\right)-\mathrm{8}{\left(u-\frac{\mathrm{1}}{\mathrm{8}}{b}_{\mathrm{1}}\right)}^{\mathrm{2}}+\epsilon \left(\mathrm{2}{b}_{\mathrm{1}}v-\mathrm{8}uv+\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{3}}\xi (\tau )\right)+o({\epsilon }^{\mathrm{2}})\mathrm{。}\end{array}\end{array}\right.

因此，设\eta =u-b_{\mathrm{1}}/\mathrm{8}， 得到

\left\{\begin{array}{l}\frac{\mathrm{d}\eta }{\mathrm{d}\tau }=v,\\ \frac{\mathrm{d}v}{\mathrm{d}\tau }=\left(\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{8}}{b_{\mathrm{1}}}^{\mathrm{2}}\right)-\mathrm{8}{\eta }^{\mathrm{2}}+\epsilon \left(b_{\mathrm{1}}v-\mathrm{8}\eta v+\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{3}}\xi (\tau )\right)+o({\epsilon }^{\mathrm{2}})\mathrm{。}\end{array}\right.（11）

接下来，进一步缩放，将系统（11）化为标准形式。设x=-\mathrm{8}\eta , y=-\mathrm{8}v，则系统（11）化为

\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}\tau }=y,\\ \frac{\mathrm{d}y}{\mathrm{d}\tau }={b_{\mathrm{1}}}^{\mathrm{2}}-\mathrm{4}{b}_{\mathrm{2}}+x^{\mathrm{2}}+\epsilon \left(b_{\mathrm{1}}y+xy-\mathrm{4}{b}_{\mathrm{3}}\xi (\tau )\right)+o({\epsilon }^{\mathrm{2}})\mathrm{。}\end{array}\right.（12）

设\mathrm{4}{b}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}}=c^{\mathrm{2}}，由已知条件知其为严格正常数（恒正常数）。假设c>\mathrm{0}，同时假设b_{\mathrm{1}}>\mathrm{0}（时动力学无意义），且b_{\mathrm{3}}=\mathrm{0}（系统（12）此时为Takens-Bogdanov分岔的标准范式）。因而系统（12）得到两个稳定点：S_{\pm }(\pm c, \mathrm{0})，S_{+}为鞍点且：

（i） 若\mathrm{2}{b_{\mathrm{1}}}^{\mathrm{2}}>\mathrm{4}{b}_{\mathrm{2}}>{b_{\mathrm{1}}}^{\mathrm{2}}，则S_{-}是不稳定焦点；

（ii）若l(b_{\mathrm{1}})>\mathrm{4}{b}_{\mathrm{2}}>\mathrm{2}{b_{\mathrm{1}}}^{\mathrm{2}}，这里l(b_{\mathrm{1}})\approx (\mathrm{74}/\mathrm{25}){b_{\mathrm{1}}}^{\mathrm{2}}，则S_{-}是稳定焦点，且被不稳定的周期轨道包围；

（iii）若\mathrm{4}{b}_{\mathrm{2}}>l(b_{\mathrm{1}})，则S_{-}是稳定焦点，且没有围绕它的周期轨道。

{b_{\mathrm{1}}}^{\mathrm{2}}=\mathrm{2}{b}_{\mathrm{2}}处的分岔是标准（亚临界）Hopf分岔，在l(b_{\mathrm{1}})=\mathrm{4}{b}_{\mathrm{2}}处破坏了周期轨道的分岔是同宿（全局）分岔。图2（a）—图2（c）展示的是b_{\mathrm{3}}=\mathrm{0}时的系统（12）的相图。图2（a）中，S_{+}(\sqrt[]{\mathrm{2}}, \mathrm{0})是鞍点，S_{-}(-\sqrt[]{\mathrm{2}}, \mathrm{0})是不稳定焦点。图2（b）中，S_{+}(\sqrt[]{\mathrm{6}}, \mathrm{0})是鞍点，S_{-}(-\sqrt[]{\mathrm{6}}, \mathrm{0})是稳定焦点，被不稳定的周期轨道所环绕。图2（c）中，S_{+}(\mathrm{2}\sqrt[]{\mathrm{2}}, \mathrm{0})是鞍点，S_{-}(-\mathrm{2}\sqrt[]{\mathrm{2}}, \mathrm{0})是稳定焦点，周围没有周期轨道所环绕。

3　Melnikov混沌

当\epsilon足够小时，我们可以将Melnikov方法应用于形式为（12）的SIR模型。若\epsilon =\mathrm{0}， 则系统（12）是未受干扰的系统，即

\left\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}\tau }=y,\\ \frac{\mathrm{d}y}{\mathrm{d}\tau }=-(\mathrm{4}{b}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}})+x^{\mathrm{2}}\mathrm{。}\end{array}\right.（13）

系统（13）是Hamiltonian系统，其Hamiltonian函数为

H(x,y)=\frac{\mathrm{1}}{\mathrm{2}}{y}^{\mathrm{2}}+(\mathrm{4}{b}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}})x-\frac{\mathrm{1}}{\mathrm{3}}{x}^{\mathrm{3}}\mathrm{。}

势能函数为V(x)=(\mathrm{4}{b}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}})x-\frac{\mathrm{1}}{\mathrm{3}}{x}^{\mathrm{3}}。 Jacobian矩阵为\left(\begin{array}{cc}\mathrm{0}& \mathrm{1}\\ \mathrm{2}x& \mathrm{0}\end{array}\right)， 特征多项式为{\lambda }^{\mathrm{2}}-\mathrm{2}x=\mathrm{0}。系统（13）有两个稳定点：S(c, \mathrm{0})是鞍点，C(-c, \mathrm{0})是中心，其被一组连续的周期轨道\left(x_T(\tau ), y_T(\tau )\right)包围，并且存在有界的同宿轨道\left(x_h(\tau ), y_h(\tau )\right)满足：

\underset{\tau \to \pm \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\left(x_h(\tau ), y_h(\tau )\right)=(c, \mathrm{0})\mathrm{。}

同宿轨道表达式^［24］为

\left\{\begin{array}{l}{x}_h(\tau )=c-\mathrm{3}c\mathrm{c}\mathrm{o}\mathrm{s}{\mathrm{h}}^{-\mathrm{2}}\left(\tau \sqrt[]{\frac{c}{\mathrm{2}}}\right),\\ y_h(\tau )=\mathrm{3}c\sqrt[]{\mathrm{2}c}\mathrm{c}\mathrm{o}\mathrm{s}{\mathrm{h}}^{-\mathrm{2}}\left(\tau \sqrt[]{\frac{c}{\mathrm{2}}}\right)\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\tau \sqrt[]{\frac{c}{\mathrm{2}}}\right)\mathrm{。}\end{array}\right.

图3（a）和图3（b）分别显示了当b_{\mathrm{1}}=\mathrm{0}和b_{\mathrm{2}}=\mathrm{1}时的相图和势能图。图3（a）中，C(-\mathrm{2}, \mathrm{0})是中心，S(\mathrm{2}, \mathrm{0})是鞍点，且通过同宿轨道{\mathrm{\Gamma }}_{\mathrm{h}\mathrm{o}\mathrm{m}}连接自身。

有界噪声干扰传染率的随机Melnikov函数为

\begin{array}{c}M(t_{\mathrm{0}})={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{y}_h(t)\left(b_{\mathrm{1}}{y}_h(t)+x_h(t)y_h(t)-\mathrm{4}{b}_{\mathrm{3}}\xi (t+t_{\mathrm{0}})\right)\mathrm{d}t=\\ \frac{\mathrm{6}}{\mathrm{5}}{(\mathrm{2}c)}^{\mathrm{5}/\mathrm{2}}{b}_{\mathrm{1}}-\frac{\mathrm{3}}{\mathrm{7}}{(\mathrm{2}c)}^{\mathrm{7}/\mathrm{2}}-{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{4}{b}_{\mathrm{3}}{y}_h(t)\xi (t+t_{\mathrm{0}})\mathrm{d}t=I+Z(t_{\mathrm{0}}),\end{array}（14）

其中I=\frac{\mathrm{6}}{\mathrm{5}}{(\mathrm{2}c)}^{\mathrm{5}/\mathrm{2}}{b}_{\mathrm{1}}-\frac{\mathrm{3}}{\mathrm{7}}{(\mathrm{2}c)}^{\mathrm{7}/\mathrm{2}}, Z(t_{\mathrm{0}})=-{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{4}{b}_{\mathrm{3}}{y}_h(t)\xi (t+t_{\mathrm{0}})\mathrm{d}t。这里t_{\mathrm{0}}表示Poincaré映射的横截面时间，可以看作是有界噪声扰动的初始时间。如果Melnikov函数存在简单零点，那么Poincaré映射的稳定和不稳定流形便横截相交，此时便会出现Smale马蹄形变换意义下的混沌。（14）式中前两个积分表示由阻尼产生的Melnikov过程的均值，最后一个积分表示由有界噪声产生的Melnikov过程的随机部分。我们可以从概率统计角度分析M(t_{\mathrm{0}})的简单零点。假设随机过程E[Z(t_{\mathrm{0}})]的均值为零，则E[M(t_{\mathrm{0}})]=I，若是一个非零常数，则在平均意义下系统（12）永远不会发生混沌。若它等于零，则在平均意义下系统（12）产生混沌的条件是：\frac{\mathrm{6}}{\mathrm{5}}{(\mathrm{2}c)}^{\mathrm{5}/\mathrm{2}}{b}_{\mathrm{1}}-\frac{\mathrm{3}}{\mathrm{7}}{(\mathrm{2}c)}^{\mathrm{7}/\mathrm{2}}=\mathrm{0}。图4展示了在(b_{\mathrm{1}}, b_{\mathrm{2}} , t_{\mathrm{0}})空间中I=\mathrm{0}表示的曲面。

接下来考虑Melnikov过程（14）在均方意义下是否存在简单零点。脉冲响应函数为

h(t)=y_h(t)=\mathrm{3}c\sqrt[]{\mathrm{2}c}\mathrm{c}\mathrm{o}\mathrm{s}{\mathrm{h}}^{-\mathrm{2}}(t\sqrt[]{c/\mathrm{2}})\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(t\sqrt[]{c/\mathrm{2}}),

频率响应函数为

\begin{array}{c}H(\omega )={\int }_{-\infty }^{+\infty }h(t)exp(-j\omega t)dt={\int }_{-\infty }^{+\infty }\mathrm{3}c\sqrt[]{\mathrm{2}c}cos{h}^{-\mathrm{2}}(t\sqrt[]{c/\mathrm{2}})tanh(t\sqrt[]{c/\mathrm{2}})exp(-j\omega t)dt=\\ \frac{\mathrm{j}{\mathrm{\pi }}^{\mathrm{2}}\omega }{\mathrm{2}\left(-\frac{\mathrm{1}}{\mathrm{3}\sqrt[]{\mathrm{2}c}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)},\end{array}

还需要考虑{\left[M(t_{\mathrm{0}})-E\left(M(t_{\mathrm{0}})\right)\right]}^{\mathrm{2}}和{\left[\mathrm{0}-E\left(M(t_{\mathrm{0}})\right)\right]}^{\mathrm{2}}的均值是否相等。因为

\begin{array}{c}E{\left[M(t_{\mathrm{0}})-E\left(M(t_{\mathrm{0}})\right)\right]}^{\mathrm{2}}=E{\left[Z(t_{\mathrm{0}})\right]}^{\mathrm{2}}={{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left|H(\omega )\right|}^{\mathrm{2}}{S}_{\xi }(\omega )\mathrm{d}\omega =\\ {{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left|H(\omega )\right|}^{\mathrm{2}}\frac{(b_{\mathrm{3}}{\sigma )}^{\mathrm{2}}}{\mathrm{2}\mathrm{\pi }}\left(\frac{\mathrm{1}}{\mathrm{4}(\omega -{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}+\frac{\mathrm{1}}{\mathrm{4}(\omega +{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}\right)\mathrm{d}\omega ={{\sigma }_Z}^{\mathrm{2}}\end{array}

且E{\left[\mathrm{0}-E\left(M(t_{\mathrm{0}})\right)\right]}^{\mathrm{2}}=E{\left[M(t_{\mathrm{0}})\right]}^{\mathrm{2}}=I^{\mathrm{2}}， 因而若满足

I^{\mathrm{2}}={{\sigma }_Z}^{\mathrm{2}}，（15）

则Melnikov过程（14）在均方意义下存在简单零点。

定理 1 基于Melnikov过程，均方意义下若满足I^{\mathrm{2}}={{\sigma }_Z}^{\mathrm{2}}，即

\begin{array}{c}{\left[\frac{\mathrm{6}}{\mathrm{5}}{(\mathrm{2}c)}^{\mathrm{5}/\mathrm{2}}{b}_{\mathrm{1}}-\frac{\mathrm{3}}{\mathrm{7}}{(\mathrm{2}c)}^{\mathrm{7}/\mathrm{2}}\right]}^{\mathrm{2}}=\\ {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\left(\frac{\mathrm{j}\mathrm{\pi }{\omega }^{\mathrm{2}}}{\mathrm{2}\left(-\frac{\mathrm{1}}{\mathrm{3}\sqrt[]{\mathrm{2}c}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)}\right)}^{\mathrm{2}}\times \frac{(b_{\mathrm{3}}{\sigma )}^{\mathrm{2}}}{\mathrm{2}\mathrm{\pi }}\left(\frac{\mathrm{1}}{\mathrm{4}(\omega -{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}+\frac{\mathrm{1}}{\mathrm{4}(\omega +{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}\right)\mathrm{d}\omega =\\ \frac{\mathrm{9}\mathrm{\pi }{b_{\mathrm{3}}}^{\mathrm{2}}c}{\mathrm{4}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{{\sigma }^{\mathrm{2}}{\omega }^{\mathrm{4}}}{\mathrm{s}\mathrm{i}\mathrm{n}{\mathrm{h}}^{\mathrm{2}}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\mathrm{h}}^{\mathrm{2}}\left(\frac{\mathrm{\pi }\omega }{\mathrm{4}\sqrt[]{c/\mathrm{2}}}\right)}\left(\frac{\mathrm{1}}{\mathrm{4}(\omega -{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}+\frac{\mathrm{1}}{\mathrm{4}(\omega +{\Omega }_{\mathrm{1}}{)}^{\mathrm{2}}+{\sigma }^{\mathrm{4}}}\right)\mathrm{d}\omega ,\end{array}（16）

则模型（1）可能存在混沌。

（16）式中的积分可以通过数值方法计算。令b_{\mathrm{1}}=\mathrm{2}, b_{\mathrm{2}}=\mathrm{3}/\mathrm{2}, {\mathrm{\Omega }}_{\mathrm{1}}=\mathrm{2}。（15）式中的均方准则通过数值计算得到系统（12）中混沌发生的有界噪声的振幅阈值，如图5（a）和图5（b），图5（a）展示了\sigma \in [\mathrm{0}, \mathrm{3}]上振幅阈值曲线，图5（b）是图5（a）在\sigma \in [\mathrm{0}, \mathrm{1}]上放大的振幅阈值曲线，在这些曲线的上方可能发生混沌。图6（a）和图6（b）分别显示了系统（12）在b_{\mathrm{3}}=\mathrm{7.5}, \sigma =\mathrm{0.5}和b_{\mathrm{3}}=\mathrm{10}, \sigma =\mathrm{2}时的相图，我们可以看到混沌轨道。

4　结论

本文研究了有界噪声对一类具有非线性发生率的SIR模型混沌运动的影响。基于同宿分岔，利用随机Melnikov理论得到了可能发生混沌运动的充分条件，并建立了有界噪声强度使混沌发生的阈值。结果表明，噪声的存在会增大阈值幅值，在参数空间中使可能出现的混沌区域缩小。对原系统相图的数值计算也验证了噪声强度越大，混沌发生的阈值的幅值会随着噪声强度的增加而上移。结果揭示了混沌吸引子受到有界噪声的扩散，并且噪声强度越大，吸引子越扩散。

参考文献

[1] 马知恩, 周义仓, 王稳地, 等. 传染病动力学的数学建模与研究[M]. 北京: 科学出版社, 2004.

[2] MA Z E, ZHOU Y C, WANG W D, et al. Mathematical Modeling and Research on Dynamics of Infectious Diseases[M]. Beijing: Science Press, 2004.

[3] 靳祯, 孙桂全, 刘茂省. 网络传染病动力学建模与分析[M]. 北京: 科学出版社, 2014.

[4] JIN Z, SUN G Q, LIU M S. Dynamic Modeling and Analysis of Network Infectious Diseases[M]. Beijing: Science Press, 2014.

[5] 季春燕, 蒋达清. 随机生物模型和传染病模型[M]. 北京: 科学出版社, 2018.

[6] JI C Y, JIANG D Q. Stochastic Biological Model and Infectious Disease Model[M]. Beijing: Science Press, 2018.

[7] DERRICK W R, VAN DEN DRIESSCHE P. Homoclinic Orbits in a Disease Transmission Model with Nonlinear Incidence and Nonconstant Population[J]. Discrete Contin Dyn Syst B, 2003, 3(2): 299-309. DOI: 10.3934/dcdsb.2003.3.299 .

[8] RUAN S G, WANG W D. Dynamical Behavior of an Epidemic Model with a Nonlinear Incidence Rate[J]. J Differ Equ, 2003, 188(1): 135-163. DOI: 10.1016/S0022-0396(02)00089-X .

[9] VAN DEN DRIESSCHE P, WATMOUGH J. A Simple SIS Epidemic Model with a Backward Bifurcation[J]. J Math Biol, 2000, 40(6): 525-540. DOI: 10.1007/s002850000032 .

[10] ALEXANDER M E, MOGHADAS S M. Periodicity in an Epidemic Model with a Generalized Non-linear Incidence[J]. Math Biosci, 2004, 189(1): 75-96. DOI: 10.1016/j.mbs.2004.01.003 .

[11] JIN Y, WANG W D, XIAO S W. An SIRS Model with a Nonlinear Incidence Rate[J]. Chaos Solitons Fractals, 2007, 34(5): 1482-1497. DOI: 10.1016/j.chaos.2006.04.022 .

[12] JIANG D Q, YU J J, JI C Y, et al. Asymptotic Behavior of Global Positive Solution to a Stochastic SIR Model[J]. Math Comput Model, 2011, 54(1/2): 221-232. DOI: 10.1016/j.mcm.2011.02.004 .

[13] LAHROUZ A, OMARI L. Extinction and Stationary Distribution of a Stochastic SIRS Epidemic Model with Non-linear Incidence[J]. Stat Probab Lett, 2013, 83(4): 960-968. DOI: 10.1016/j.spl.2012.12.021 .

[14] CAI Y L, JIAO J J, GUI Z J, et al. Environmental Variability in a Stochastic Epidemic Model[J]. Appl Math Comput, 2018, 329: 210-226. DOI: 10.1016/j.amc.2018.02.009 .

[15] 曹虹, 樊小琳, 聂麟飞. 一类具有非线性发生率的随机SIRS传染病模型研究[J]. 数学的实践与认识, 2024, 54(10): 144-151.

[16] CAO H, FAN X L, NIE L F. A Study on a Stochastic SIRS Epidemic Model with Nonlinear Incidence[J]. Math Pract Theory, 2024, 54(10): 144-151.

[17] ARON J L, SCHWARTZ I B. Seasonality and Period-doubling Bifurcations in an Epidemic Model[J]. J Theor Biol, 1984, 110(4): 665-679. DOI: 10.1016/S0022-5193(84)80150-2 .

[18] KUZNETSOV Y A, PICCARDI C. Bifurcation Analysis Pf Periodic SEIR and SIR Epidemic Models[J]. J Math Biol, 1994, 32(2): 109-121. DOI: 10.1007/BF00163027 .

[19] SCHAFFER W M, KOT M. Nearly One Dimensional Dynamics in an Epidemic[J]. J Theor Biol, 1985, 112(2): 403-427. DOI: 10.1016/S0022-5193(85)80294-0 .

[20] ELLNER S, TURCHIN P. Chaos in a Noisy World: New Methods and Evidence from Time-series Analysis[J]. Am Nat, 1995, 145(3): 343-375. DOI: 10.1086/285744 .

[21] GLENDINNING P, PERRY L P. Melnikov Analysis of Chaos in a Simple Epidemiological Model[J]. J Math Biol, 1997, 35(3): 359-373. DOI: 10.1007/s002850050056 .

[22] DIALLO O, KONÉ Y. Melnikov Analysis of Chaos in a General Epidemiological Model[J]. Nonlinear Anal Real World Appl, 2007, 8(1): 20-26. DOI: 10.1016/j.nonrwa.2005.03.032 .

[23] DIALLO O, KONÉ Y, MAIGA A. Melnikov Analysis of Chaos in an Epidemiological Model with Almost Periodic Incidence Rates[J]. Appl Math Sci, 2008, 2(28): 1377-1386.

[24] WIGGINS S. Introduction to Applied Nonlinear Dynamical Systems and Chaos[M]. Cham: Springer New York, 2003. DOI: 10.1007/b97481 .

[25] LIU W Y, ZHU W Q, HUANG Z L. Effect of Bounded Noise on Chaotic Motion of Duffing Oscillator Under Parametric Excitation[J]. Chaos Solitons Fractals, 2001, 12(3): 527-537. DOI: 10.1016/S0960-0779(00)00002-3 .

[26] 朱培, 任兴民, 秦卫阳, 等. 非线性双稳态压电俘能系统在有界噪声与谐和激励下的Melnikov混沌运动阈值分析[J]. 西北工业大学学报, 2017, 35(4): 669-675. DOI: 10.3969/j.issn.1000-2758.2017.04.017 .

[27] ZHU P, REN X M, QIN W Y, et al. The Analysis of Melnikov Chaotic Motion Threshold of Nonlinear Bistable Piezoelectric Energy Harvester System Subjected Bounded Noise and Harmonic Excitation[J]. J Northwest Polytech Univ, 2017, 35(4): 669-675. DOI: 10.3969/j.issn.1000-2758.2017.04.017 .

[28] 孔令辉, 刘丁杨, 蹇开林. 非均质物理双摆的混沌特性研究[J]. 重庆大学学报, 2024, 47(2): 106-118. DOI: 10.11835/j.issn.1000-582X.2023.255 .

[29] KONG L H, LIU D Y, JIAN K L. Study on Chaotic Characteristics of Heterogeneous Physical Double Pendulum[J]. J Chongqing Univ, 2024, 47(2): 106-118. DOI: 10.11835/j.issn.1000-582X.2023.255 .

[30] GUCKENHEIMER J, HOLMES P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[M]. Cham: Springer New York, 1983.

基金资助

国家自然科学基金(12571537)

山西省自然科学基金(201901D111041)

山西省自然科学基金(201601D202002)