三维Chen系统的动力学行为分析

Analysis of Dynamical Behavior in a Three-Dimensional Chen System

摘要: 本文利用谱理论和动态跃迁理论研究了三维Chen系统的动力学行为。利用谱理论，得到了该系统发生两次动态跃迁的临界控制参数c₁和c₂，并给出了相应的PES条件。基于动态跃迁理论，得到了该系统发生两次跃迁的相关结论。首先，该系统在c₁处发生连续型跃迁，原有平衡态失稳，并分歧出两个渐近稳定的奇点。其次，随着控制参数c的持续增大，该系统在c₂处发生两种不同类型的跃迁。在一定条件下，该系统发生跳跃型跃迁，并分歧出不稳定的周期轨道；而在相反的条件下发生连续型跃迁，分歧出渐近稳定的周期轨道。

Abstract: This paper investigates the dynamical behavior of the three-dimensional Chen system using spectral theory and dynamic transition theory. By applying spectral theory, the critical control parameters c₁ and c₂ for the two dynamic transitions of the system are obtained, along with the corresponding PES conditions. Based on dynamic transition theory, conclusions regarding the two transitions are derived. First, the system undergoes a continuous transition at c₁, where the original equilibrium state loses stability and bifurcates into two asymptotically stable singular points. Second, as the control parameter c continues to increase, the system undergoes two different types of transitions at c₂. Under certain conditions, the system undergoes a jump transition and bifurcates into an unstable periodic orbit; under the opposite conditions, a continuous transition occurs, bifurcating into an asymptotically stable periodic orbit.

1. 引言

在自然科学和生物网络等技术领域中，有许多现象无法用线性理论解释，这就促进了非线性科学的出现和蓬勃发展，混沌系统作为非线性科学的核心研究领域，因其深刻揭示了确定性系统中蕴含的内在随机性与复杂动力学行为，一直备受学者关注。自Lorenz系统[1]开启了现代混沌研究的先河以来，一系列三维自治混沌系统，如Rössler系统[2]、Chen系统[3]等被相继发现并被深入研究。1999年，Chen构造了下述三维混沌系统，称为Chen系统[3]，

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a\left(y-x\right),\\ \frac{\text{d}y}{\text{d}t}=\left(c-a\right)x+cy-xz,\\ \frac{\text{d}z}{\text{d}t}=-bz+xy.\end{array}$ (1)

其中$x,y,z$ 是系统变量，常数$a,b,c$ 是系统参数，该系统作为Lorenz系统[1]的对偶系统，在特定参数范围内具有更复杂的拓扑结构。

Ueta和Chen分析了该系统的分岔和动力学性质[4]。Li和Chen借助第一Lyapunov系数研究了Chen系统的Hopf分岔[5]。Chang和Chen利用分岔理论和中心流形定理，给出了Chen系统的Hopf分岔条件[6]。Huang和Liu运用Kosambi-Cartan-Chern几何理论对Chen系统的雅可比稳定性进行了研究，得到原点始终是雅可比不稳定的，另外两个非零平衡点的雅可比稳定性取决于参数的取值，且Chen系统的周期轨道也被证明是雅可比不稳定的[7]。Vijayalakshmi等人利用可微函数的分数阶导数，得到了Chen系统的拉格朗日函数的表达式，并给出Chen系统相关的守恒量的估计[8]。Ignatyev考虑了具有脉冲作用的Chen系统，因为Chen系统的零解在Lyapunov意义上是不稳定的。Ignatyev通过在固定时刻使用脉冲作用或使用在一组特定相空间上发生的脉冲作用使得Chen系统的零解在Lyapunov意义上稳定[9]。

在已有的研究基础上，本文关注三维Chen系统随参数变化而发生的分歧与跃迁现象。传统的分析方法多集中于数值模拟与线性稳定分析，而对于系统在临界参数点附近状态演化的严格数学描述，尤其是对新状态(如平衡点、周期解)的稳定性与跃迁类型的判定，则需要更强大的理论工具。本文采用Ma和Wang提出的动态跃迁理论[10] [11]展开研究。借助谱理论得到了Chen系统发生第一次跃迁和第二次跃迁的临界控制参数；借助动态跃迁理论得到了该系统两次跃迁出的新的稳态解的近似表达式及其稳定性。

2. 预备知识

2.1. 动态跃迁理论

从自然背景角度来说，动态跃迁理论实质上反映了自然界中普遍存在的一种状态变迁现象。自然界中的演化系统总是受到一些因素(在数学中称为参数)的控制，当这些因素相互之间平衡时系统会处于一个动态稳定状态，当旧的平衡被破坏，参变量越过某个临界值时，旧的平衡态不再稳定，系统会跃迁到一个新的状态，这就是跃迁的内容[12]。

值得注意的是，自然界中的非线性耗散系统一般可由下面的非线性演化方程描述，

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=L_{\lambda }u+G\left(u,\lambda \right),\\ u\left(0\right)=\phi \left(x\right).\end{array}$ (2)

其中，$X,X_1$ 是两个Banach空间，$X_1\to X$ 是稠密紧包含。$u:\left[0,+\infty \right)\to X_1$ 是未知函数。$\lambda \in {ℝ}^1$ 为系统参数。$L_{\lambda }:X_1\to X$ 是含参数的线性全连续场，且连续依赖于$\lambda \in {ℝ}^1$ ，并且满足

$L_{\lambda }=-A_{\lambda }+B:X_1\to X$ 是一个扇形算子。 (3)

其中，$A_{\lambda }:X_1\to X$ 是一个线性同胚，$B:X_1\to X$ 是一个线性紧算子。

对某个$0\le \theta <1$ ，假设$G\left(\cdot ,\lambda \right):X_{\theta }\to X$ 是$C^r\left(r\ge 1\right)$ 有界映射，连续依赖于$\lambda \in {ℝ}^1$ ，并满足

$G\left(u,\lambda \right)=o\left({\Vert u\Vert }_{X_{\theta }}\right),0\le \theta <1,\forall \lambda \in {ℝ}^1.$ (4)

其中，$X_{\zeta }$ 是由扇形算子$L_{\lambda }$ 决定的分数次空间。条件(3)和(4)可证明系统(2)的耗散结构。因此，我们假设条件(3)与(4)恒成立。接下来介绍系统(2)发生动态跃迁的数学定义。

定义2.1 [12]如果满足以下条件：

1) 当$\lambda <{\lambda }_0$ 时，$u=0$ 关于方程(2)是局部渐近稳定的；

2) 当$\lambda >{\lambda }_0$ 或$\lambda <{\lambda }_0$ 时，则存在一个独立于${\lambda }_0$ 的${\Sigma }_{\lambda }$ 邻域$U$ ，使得对于任意$y\in U\backslash {\Gamma }_{\lambda }$ ，系统(2)的初始解$u_{\lambda }\left(t,y\right)$ 满足

$\{\begin{array}{l}\overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,y\right)\Vert }_X\ge \delta \left(\lambda \right)>0,\\ \underset{\lambda \to {\lambda }_0+0}{\lim }\delta \left(\lambda \right)=\delta \ge 0.\end{array}$或$\{\begin{array}{l}\overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,y\right)\Vert }_X\ge \delta \left(\lambda \right)>0,\\ \underset{\lambda \to {\lambda }_0-0}{\lim }\delta \left(\lambda \right)>\delta \le 0.\end{array}$

其中，${\Gamma }_{\lambda }$ 是$u=0$ 的稳定流形，对于$\lambda >{\lambda }_0$ 或$\lambda <{\lambda }_0$ ，${\Gamma }_{\lambda }$ 在$U$ 中有余维codim ${\Gamma }_{\lambda }\ge 1$ 。

紧接着介绍一个关于跃迁的基本定理，它为非线性耗散系统(2)提供了跃迁的充分条件和基本分类。令$L_{\lambda }$ 的特征值为$\left\{{\beta }_i\left(\lambda \right)|i=1,2,\cdots \right\}\subseteq ℂ$ ，假设它们满足

$\begin{array}{l}\mathrm{Re}{\beta }_i\left(\lambda \right)\{\begin{array}{l}<0,\lambda <{\lambda }_0,\\ =0,\lambda ={\lambda }_0,\\ >0,\lambda >{\lambda }_0,\end{array}\forall 1\ge i\ge m,\\ \mathrm{Re}{\beta }_i\left({\lambda }_0\right)<0,\forall j\le m+1.\end{array}$ (5)

定理2.1 [12]在条件(5)的假设下，方程(2)从$\left(u,\lambda \right)=\left(0,{\lambda }_0\right)$ 在$\lambda >{\lambda }_0$ 有一个跃迁，并且该跃迁是下面三种类型之一：

1) 连续型跃迁：对$u=0$ 的一个邻域$U\subset X$ 及任何$\phi \in U$ ，方程(2)的解$u_{\lambda }\left(t,\phi \right)$ 满足

$\underset{\lambda \to {\lambda }_0}{\lim }\overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,\phi \right)\Vert }_X=0.$

或者等价的，方程(2)在$\left(0,{\lambda }_0\right)$ 有一个吸引子分歧。

2) 跳跃型跃迁：存在$\epsilon >0$ ，对任意${\lambda }_0<\lambda <{\lambda }_0+\epsilon$ ，存在$U$ 中的一个开稠集$U_{\lambda }\subset U$ ，使得对任意$\phi \in U_{\lambda }$ 有

$\overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,\phi \right)\Vert }_X\ge \delta >0.$

其中，$\delta >0$ 是一个与$\lambda$ 无关的常数。这也称为不连续型跃迁。

3) 混合型跃迁：对任意${\lambda }_0<\lambda <{\lambda }_0+\epsilon$ ，$U$ 能够分解为两个开集$U_1^{\lambda }$ 和$U_2^{\lambda }$ ($U_2^{\lambda }$ 不一定连通)，

$U\backslash \Sigma =U_1^{\lambda }+U_2^{\lambda },U_1^{\lambda }\cap U_1^{\lambda }=\varnothing .$

其中，$\Sigma \subset U$ 是一个余维codim $\Sigma \ge 1$ 的集合，使得(2)的解$u_{\lambda }\left(t,\phi \right)$ 满足

$\{\begin{array}{l}\underset{\phi \in U_1^{\lambda }}{\sup }\overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,\phi \right)\Vert }_X\to 0,\lambda \to {\lambda }_0+0,\\ \overline{\underset{t\to \infty }{\lim }}{\Vert u_{\lambda }\left(t,\phi \right)\Vert }_X\to \delta >0,\forall \phi \in U_2^{\lambda }.\end{array}$

其中，$U_1^{\lambda }$ 称为稳定型区域，$U_2^{\lambda }$ 称为跳跃型区域。

2.2. Routh-Hurwitz判据

定理2.2考虑如下非线性系统，

$\dot{x}=f\left(x,t\right),x\in {ℝ}^n.$ (6)

设在其平衡点上的特征方程的一般式展开式为，

${\lambda }^n+a_1{\lambda }^{n-1}+a_2{\lambda }^{n-2}+\cdots +a_{n-1}\lambda +a_n=0.$ (7)

将上述系数整理成n阶矩阵D，

$D=\left[\begin{array}{cccccccc}{a}_1& 1& 0& 0& 0& 0& \cdots & 0\\ a_3& a_2& a_1& 1& 0& 0& \cdots & 0\\ a_5& a_4& a_3& a_2& a_1& 1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& & ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & a_n\end{array}\right],$ (8)

我们称D的n阶主子式为Routh-Hurwitz行列式，即，

${\Delta }_1=a_1,{\Delta }_2=\left|\begin{array}{cc}{a}_1& 1\\ a_3& a_2\end{array}\right|,{\Delta }_3=\left|\begin{array}{ccc}{a}_1& 1& 0\\ a_3& a_2& a_1\\ a_5& a_4& a_3\end{array}\right|,\cdots$ (9)

方程特征根的实部为负数的充要条件为：n阶Routh-Hurwitz行列式均大于0，即，

${\Delta }_i>0,\left(i=1,2,3,\cdots \right).$ (10)

3. Chen系统的跃迁

在这一部分，本文主要讨论三维Chen系统(1)的动态跃迁。不难计算，$\left(x,y,z\right)=\left(0,0,0\right)$ 为系统(1)的常数稳态解。接下来，本文主要考虑随着控制参数c的变化，系统(1)从$\left(0,0,0\right)$ 处发生的一次动态跃迁的临界条件及其类型；进一步，讨论该系统发生二次跃迁的动力学行为。

3.1. 一次跃迁

3.1.1. PES条件

为了讨论的方便，将系统(1)表示为抽象演化方程的形式。记$\omega ={\left(x,y,z\right)}^{\text{T}}$ ，线性部分表示为

$L\omega =\left(\begin{array}{c}-ax+ay\\ \left(c-a\right)x+cy\\ -bz\end{array}\right),$ (11)

其中，线性矩阵L表示如下，

$L=\left(\begin{array}{ccc}-a& a& 0\\ c-a& c& 0\\ 0& 0& -b\end{array}\right).$

非线性部分为

$G\left(\omega \right)=\left(\begin{array}{c}0\\ -xz\\ xy\end{array}\right).$ (12)

则系统(1)可表示为如下的抽象形式，

$\frac{\text{d}\omega }{\text{d}t}=L\omega +G\left(\omega \right).$

跃迁现象主要由演化系统线性部分的特征值的临界穿越决定。接下来，主要讨论矩阵L的特征值及其符号变化。记矩阵L的特征值为${\beta }_i\left(i=1,2,3\right)$ ，通过计算可以得到，

$\begin{array}{c}\left|\beta I-L\right|=\left|\begin{array}{ccc}\beta +a& -a& 0\\ a-c& \beta -c& 0\\ 0& 0& \beta +b\end{array}\right|\\ =\left(\beta +b\right)\left({\beta }^2+\left(a-c\right)\beta +a^2-2ac\right).\end{array}$ (13)

令$\left|\beta I-L\right|=0$ ，不难得到其中一个特征值${\beta }_3=-b$ 。而另外两个特征值${\beta }_1$ ，${\beta }_2$ 满足下述方程，

${\beta }^2+\left(a-c\right)\beta +a^2-2ac=0.$ (14)

不难计算，当$c>\left(2\sqrt{3}-3\right)a$ 时，判别式$\Delta ={\left(a-c\right)}^2-4\left(a^2-2ac\right)=4c^2-3{\left(a-c\right)}^2>0$ 。此时方程(14)有两个不相等的实根，

$\begin{array}{l}{\beta }_1=\frac{1}{2}\left[c-a+\sqrt{4c^2-3{\left(a-c\right)}^2}\right]=\frac{1}{2}\left[c-a+\sqrt{-3a^2+6ac+c^2}\right],\\ {\beta }_2=\frac{1}{2}\left[c-a-\sqrt{4c^2-3{\left(a-c\right)}^2}\right]=\frac{1}{2}\left[c-a-\sqrt{-3a^2+6ac+c^2}\right].\end{array}$

则矩阵L的特征值可表示如下，

$\begin{array}{l}{\beta }_1=\frac{1}{2}\left[c-a+\sqrt{4c^2-3{\left(a-c\right)}^2}\right]=\frac{1}{2}\left[c-a+\sqrt{-3a^2+6ac+c^2}\right],\\ {\beta }_2=\frac{1}{2}\left[c-a-\sqrt{4c^2-3{\left(a-c\right)}^2}\right]=\frac{1}{2}\left[c-a-\sqrt{-3a^2+6ac+c^2}\right],\\ {\beta }_3=-b.\end{array}$ (15)

接下来计算${\beta }_i\left(i=1,2,3\right)$ 对应的特征向量$e_i\left(i=1,2,3\right)$ ，经计算可得，

$e_1={\left(P-Q,1,0\right)}^{\text{T}},e_2={\left(P+Q,1,0\right)}^{\text{T}},e_3={\left(0,0,1\right)}^{\text{T}}.$

其中，

$P=\frac{a+c}{2\left(a-c\right)},Q=\frac{\sqrt{-3a^2+6ac+c^2}}{2\left(a-c\right)}.$

进一步，通过计算可以得到与特征向量$e_i\left(i=1,2,3\right)$ 对应的共轭向量$e_i^{\ast }\left(i=1,2,3\right)$ ，

$e_1^{*}={\left(-\frac{1}{2Q},\frac{P+Q}{2Q},0\right)}^{\text{T}},e_2^{*}={\left(\frac{1}{2Q},\frac{-P+Q}{2Q},0\right)}^{\text{T}},e_3^{*}={\left(0,0,1\right)}^{\text{T}}.$

并且，$e_i$ 与$e_j^{*}$ 满足下述关系

$\langle e_i,e_j^{*}\rangle =\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array},\left(i,j=1,2,3\right).$ (16)

由(15)中特征值的表达式可知，${\beta }_i<0,\left(i=2,3\right)$ 。随着参数的变化，${\beta }_1$ 的符号会发生正负穿越。令${\beta }_1=0$ 可得

$a-c=\sqrt{4c^2-3\left(a-c^2\right)},$

解得$c=\frac{a}{2}$ ，记$c_1=\frac{a}{2}$ 。通过分析特征值的符号可知，随着控制参数$c$ 的变化，${\beta }_1$ 会发生正负穿越并且有如下的PES条件(即稳定性交换原则)，

$\begin{array}{l}{\beta }_1\{\begin{array}{l}<0,c<c_1,\\ =0,c=c_1,\\ >0,c>c_1,\end{array}\\ {\beta }_i<0,i=2,3.\end{array}$ (17)

3.1.2. 主要结果

由PES条件(17)可知，随着控制参数的变化，系统(1)在原有的平衡态处会发生动态跃迁。在这一部分，本文将给出相应的定理及其详细的证明过程。

定理3.1依据PES条件(17)，关于系统(1)有如下结论成立：

1) 当$0<c<c_1$ 时，系统(1)原有的平衡态${\omega }_0=\left(0,0,0\right)$ 是渐近稳定的；

2) 当$c>c_1$ 时，系统(1)在$\left({\omega }_0,c_1\right)$ 处发生连续型跃迁，并在$c>c_1$ 一侧分歧出两个渐近稳定的奇点；

3) 分歧出的奇点${\omega }^{+},{\omega }^{-}$ 可表示为

$\begin{array}{l}{\omega }^{+}=\left(\sqrt{b\left(2c-a\right)},\sqrt{b\left(2c-a\right)},2c-a\right),\\ {\omega }^{-}=\left(-\sqrt{b\left(2c-a\right)},-\sqrt{b\left(2c-a\right)},2c-a\right).\end{array}$

证明3.1我们分成以下三步来进行证明。

第一步：空间分解

空间${ℝ}^3$ 可被分解为

${ℝ}^3=E_1\oplus E_2,$

其中$E_1=span\left\{e_1\right\}$ ，$E_2=span\left\{e_2,e_3\right\}$ 。

则在$c=c_1$ 附近，方程(1)的解可以表示为

$\omega ={\omega }_1+{\omega }_2,$ (18)

其中，${\omega }_1=x_1{e}_1\in E_1$ ，${\omega }_2=x_2{e}_2+x_3{e}_3\in E_2$ 。

进而，在空间$E_1$ 中，系统(1)可以约化为

$\begin{array}{c}\langle \frac{\text{d}\omega }{\text{d}t},e_1^{\ast }\rangle =\langle L\omega ,e_1^{\ast }\rangle +\langle G\left(\omega \right),e_1^{\ast }\rangle \\ ={\beta }_1{x}_1+\langle G\left(\omega \right),e_1^{\ast }\rangle .\end{array}$ (19)

第二步：方程约化

为了计算约化方程(19)中的最后一项，需要考虑中心流形函数$\Phi \left(x_1\right):E_1\to E_2$ 的影响。也就是在(19)中，取$\omega =x_1{e}_1+\Phi \left(x_1\right)$ 。

在临界参数$c_1$ 附近，基于中心流形的近似计算表达式，通过计算可得

$\Phi \left(x_1\right)=-L^{-1}{P}_{2}G\left(x_1{e}_1\right)+o\left(2\right),$ (20)

其中，$P_2$ 为Leray投影，

$x_1{e}_1=\left(\begin{array}{c}\left(P-Q\right)x_1\\ x_1\\ 0\end{array}\right).$ (21)

由(21)及非线性部分$G\left(\omega \right)$ 的表达式可得，$G\left(x_1{e}_1\right)={\left(0,0,\left(P-Q\right)x_1^2\right)}^{\text{T}}$ 。将其代入(20)可得中心流形函数的近似表达式

$\Phi \left(x_1\right)=-L_{\lambda }^{-1}{P}_{2}G\left(x_1{e}_1\right)+o\left(2\right)=-\frac{P-Q}{{\beta }_3}{x}_1^2{e}_3+o\left(2\right).$

由此，约化方程(19)可进一步表示如下，

$\frac{\text{d}{x}_1}{\text{d}t}={\beta }_1{x}_1+\langle G\left(x_1{e}_1+\Phi \left(x_1\right)\right),e_1^{\ast }\rangle ,$ (22)

其中，$x_1{e}_1+\Phi \left(x_1\right)={\left(\left(P-Q\right)x_1,x_1,-\frac{\left(P-Q\right)x_1^2}{{\beta }_3}\right)}^{\text{T}}+o\left(2\right)$ ，则

$G\left(x_1{e}_1+\Phi \left(x_1\right)\right)={\left(0,\frac{{\left(P-Q\right)}^2{x}_1^3}{{\beta }_3},\left(P-Q\right)x_1^2\right)}^{\text{T}},$

经过计算可得，

$\langle G\left(x_1{e}_1+\Phi \left(x_1\right)\right),e_1^{\ast }\rangle =-\frac{\left(P+Q\right){\left(P-Q\right)}^2}{2Qb}{x}_1^3+o\left(3\right).$ (23)

将(23)代入(22)可得下述约化方程，

$\frac{\text{d}{x}_1}{\text{d}t}={\beta }_1{x}_1+M{x}_1^3+o\left(x_1^3\right),$ (24)

其中，$M=-\frac{\left(P+Q\right){\left(P-Q\right)}^2}{2Qb}=\frac{a\left(\sqrt{-3a^2+6ac+c^2}-\left(a+c\right)\right)}{2\left(a-c\right)b\sqrt{-3a^2+6ac+c^2}}$ 。

设分子为$M_1=a\left(\sqrt{-3a^2+6ac+c^2}-\left(a+c\right)\right)$ ，分母为$M_2=2\left(a-c\right)b\sqrt{-3a^2+6ac+c^2}$ 。当$c$ 接近$\frac{a}{2}$ 时，令$c=\frac{a}{2}+\epsilon$ ，其中$\epsilon$ 是一个很小的数。则分母$M_2\approx a\left(a-c\right)b>0$ 。将$c=\frac{a}{2}+\epsilon$ 代入到分子$M_1$ 中，并利用泰勒近似展开可以得到，

$\begin{array}{c}{M}_1=a\left(\sqrt{a^2/4+7a\epsilon +{\epsilon }^2}-\left(a+a/2+\epsilon \right)\right)\\ \approx a\left(\sqrt{a^2/4+7a\epsilon }-\left(a+a/2+\epsilon \right)\right)\\ \approx a\left(a/2+7\epsilon -3a/2-\epsilon \right)\\ =-a\left(a-6\epsilon \right)<0.\end{array}$

由于$M_1<0$ ，$M_2>0$ 。所以，在$c=c_1=\frac{a}{2}$ 附近$M<0$ 。

第三步：分歧分析

已知方程(1)的分歧及其局部拓扑结构完全由约化方程(24)决定，显然，当$c>c_1$ 时，方程(24)在$c>c_1$ 一侧分歧出两个渐近稳定的奇点，且表示如下

$x_1^{\pm }=\pm \sqrt{\frac{-{\beta }_1}{M}}+o\left(3\right).$ (25)

则系统(1)在$c_1$ 处发生连续型跃迁，当$c>c_1$ 时，方程(1)在$c>c_1$ 一侧分歧出两个渐近稳定的奇点，且表示如下

$\begin{array}{l}{\omega }^{+}=\left(\sqrt{b\left(2c-a\right)},\sqrt{b\left(2c-a\right)},2c-a\right),\\ {\omega }^{-}=\left(-\sqrt{b\left(2c-a\right)},-\sqrt{b\left(2c-a\right)},2c-a\right).\end{array}$ (26)

并且${\omega }^{\pm }$ 的稳定性与$x_1^{\pm }$ 的稳定性相同。

3.2. 二次跃迁

接下来，本文继续讨论随着控制参数$c$ 的增大，系统(1)发生二次跃迁的动力学行为。因为系统(1)从${\omega }^{+}$ 和${\omega }^{-}$ 处的动态跃迁分析过程相似。因此，本文只分析系统(1)从${\omega }^{+}$ 处发生的动态跃迁。为了方便，令${\omega }^{+}=\left(B,B,C\right)$ ，并对方程(1)做如下平移。令

$x^{\prime }=x-B,y^{\prime }=y-B,z^{\prime }=z-C,$ (27)

则$x=x^{\prime }+B,y=y^{\prime }+B,z=z^{\prime }+C$ 。将其代入方程(1)，可得到平移过后的方程为，

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a\left(y-x\right),\\ \frac{\text{d}y}{\text{d}t}=-ax+cy-Bz-xz,\\ \frac{\text{d}z}{\text{d}t}=Bx+By-bz+xy.\end{array}$ (28)

则系统(1)在${\omega }^{+}$ 处发生的动态跃迁就转化为方程(28)在$\left(0,0,0\right)$ 处的动态跃迁。

为了方便，记方程(28)的线性部分为

$ℒ\omega =\left(\begin{array}{c}-ax+ay\\ -ax+cy-Bz\\ Bx+By-bz\end{array}\right),$ (29)

其中，$\omega ={\left(x,y,z\right)}^{\text{T}}$ ，矩阵

$ℒ=\left(\begin{array}{ccc}-a& -a& 0\\ -a& c& -B\\ B& B& -b\end{array}\right),$

非线性部分为

$G\left(\omega \right)=\left(\begin{array}{c}0\\ -xz\\ xy\end{array}\right).$ (30)

记$ℒ$ 的特征值为$\beta$ ，则有

$\begin{array}{c}\left|\beta I-ℒ\right|=\left|\begin{array}{ccc}\beta +a& -a& 0\\ a-c& \beta -c& 0\\ 0& 0& \beta +b\end{array}\right|\\ ={\beta }^3+\left(a+b-c\right){\beta }^2+\left(a^2+bc-ac\right)\beta +3abc-2a^{2}b\\ \triangleq {\beta }^3+m_1{\beta }^2+m_2\beta +m_3.\end{array}$ (31)

即系统(28)的特征方程如下所示

${\beta }^3+m_1{\beta }^2+m_2\beta +m_3=0,$ (32)

其中，$m_1=a+b-c,m_2=a^2+bc-ac,m_3=ab\left(3c-2a\right)$ 。

不难看出$m_i>0,i=1,2,3$ 当且仅当以下条件成立

$\{\begin{array}{l}c<a+b,\\ ac<a^2+bc,\\ 2a/3<c.\end{array}$

设${\beta }_k\left(c\right)\left(k=1,2,3\right)$ 为方程(32)的根，接下来，我们根据Routh-Hurwitz定理讨论${\beta }_k\left(c\right)\left(k=1,2,3\right)$ 的符号。已知，${\beta }_k\left(c\right)\left(k=1,2,3\right)$ 有负实部当且仅当以下条件成立

${\Delta }_1=m_1>0,{\Delta }_2=\left|\begin{array}{cc}{m}_1& 1\\ m_3& m_2\end{array}\right|>0,m_3>0.$ (33)

注意，$m_3=-{\beta }_1{\beta }_2{\beta }_3$ ，$m_3>0$ 意味着该系统不会出现简单实特征值发生正负穿越的情况。同时${\Delta }_1>0$ ，${\Delta }_2={\displaystyle {\prod }_{1=i<j}^3\left({\beta }_i+{\beta }_j\right)}$ 。因此，如果${\Delta }_2=0$ ，则系统(28)出现一对简单复特征根的实部发生正负穿越。给定$a=36$ ，$b=3$ ，数值结果表明，当$c_2=27$ 时，${\Delta }_2=0$ ，如图1所示。可以看到在$c=c_2=27$ 附近，${\Delta }_2$ 的符号由正变为负，这意味系统(28)在临界值$c_2$ 处有一对简单复特征根其实部发生正负穿越。并有下述的PES条件成立，

$\begin{array}{l}\mathrm{Re}{\beta }_i\left(c\right)\{\begin{array}{l}<0,c<c_2,\\ =0,c=c_2,\\ >0,c>c_2,\end{array}i=1,2,\\ \mathrm{Im}{\beta }_i\left(c_2\right)\ne 0,i=1,2\\ \mathrm{Re}{\beta }_3\left(c_2\right)<0.\end{array}$ (34)

Figure 1. Critical value of the second transition

图1. 二次跃迁的临界值

下面，我们讨论系统(1)发生的二次跃迁的类型。假设${\beta }_1=\gamma -i\eta ={\overline{\beta }}_2$ 。令$e_1$ 和$e_2$ 为矩阵$ℒ$ 对应于${\beta }_1$ 的特征向量的实部和虚部，$e_3$ 为对应于${\beta }_3$ 的特征向量。则空间${ℝ}^3$ 可分解为以下形式

${ℝ}^3={\tilde{E}}_1\oplus {\tilde{E}}_2,$

其中，${\tilde{E}}_1=span\left\{e_1,e_2\right\},{\tilde{E}}_2=span\left\{e_3\right\}$。

在临界参数$c_2$ 附近，方程(28)的解$\omega$ 可以被分解为

$\omega ={\omega }_1+{\omega }_2,$

其中，${\omega }_1=x_1{e}_1+x_2{e}_2\in {\tilde{E}}_1,{\omega }_2=x_3{e}_3\in {\tilde{E}}_2$。则在${\tilde{E}}_1$ 中，方程(28)可被约化为如下形式，

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=\gamma x_1-\eta x_2+\left(G\left({\omega }_1+{\omega }_2\right),e_1^{\ast }\right),\\ \frac{\text{d}{x}_2}{\text{d}t}=\gamma x_2+\eta x_1+\left(G\left({\omega }_1+{\omega }_2\right),e_2^{\ast }\right).\end{array}$ (35)

其中，$e_i^{\ast }\left(i=1,2\right)$ 是$e_i\left(i=1,2\right)$ 的共轭特征向量。

为了得到(35)中最后一项的具体表达式，我们需要计算中心流形函数$\Phi \left({\omega }_1\right):{\tilde{E}}_1\to {\tilde{E}}_2$ 。此时，系统(28)的解可近似表示为$\omega ={\omega }_1+\Phi \left({\omega }_1\right)+o\left({\omega }_1^3\right)$ 。在$c_2$ 附近，利用中心流形函数近似表达式可以得到以下简化方程，

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=\gamma x_1-\eta x_2+{\displaystyle \sum _{2\le p+q\le 3}{a}_{pq}^1{x}_1^p{x}_2^q}+o\left({\left|\left(x_1,x_2\right)\right|}^3\right),\\ \frac{\text{d}{x}_2}{\text{d}t}=\gamma x_2+\eta x_1+{\displaystyle \sum _{2\le p+q\le 3}{a}_{pq}^2{x}_1^p{x}_2^q}+o\left({\left|\left(x_1,x_2\right)\right|}^3\right).\end{array}$ (36)

其中，$\eta \left(c_2\right)\ne 0$ 并且${\displaystyle {\sum }_{2\le p+q\le 3}^{}{a}_{pq}^i{x}_1^p{x}_2^q}=\left(G\left({\omega }_1+\Phi \left({\omega }_1\right)\right),e_i^{\ast }\right)\left(i=1,2\right)$ 。

利用极坐标变换，令$\left(x_1,x_2\right)=\left(\rho \cos \theta ,\rho \sin \theta \right)$ ，并将其代入(36)可得到下述方程

$\frac{\text{d}\rho }{\text{d}\theta }=\frac{\gamma \rho +{\displaystyle {\sum }_{k=2}^3{\rho }^k{u}_k\left(\sin \theta ,\cos \theta \right)+o\left({\rho }^3\right)}}{\eta -{\displaystyle {\sum }_{k=2}^3{\rho }^{k-1}{v}_k\left(\sin \theta ,\cos \theta \right)+o\left({\rho }^2\right)}},$ (37)

其中，

$\begin{array}{l}{u}_k\left(\sin \theta ,\cos \theta \right)={\displaystyle \sum _{p+q=k}{a}_{pq}^1{\cos }^{p+1}\theta {\sin }^q\theta }+a_{pq}^2{\cos }^p\theta {\sin }^{q+1}\theta ,\\ v_k\left(\sin \theta ,\cos \theta \right)={\displaystyle \sum _{p+q=k}{a}_{pq}^1{\cos }^p\theta {\sin }^{q+1}\theta }-a_{pq}^2{\cos }^{p+1}\theta {\sin }^q\theta .\end{array}$

在$\rho =0$ 附近，方程(37)可以转化为

$\frac{1}{{\rho }^2}\frac{\text{d}\rho }{\text{d}\theta }=\frac{1}{\eta }\left[\frac{\gamma }{\rho }+{\displaystyle \sum _{k=2}^3{\rho }^{k-2}{g}_k+o\left(\rho \right)}\right],$ (38)

其中，

$\begin{array}{l}{g}_2=u_2+{\eta }^{-1}\gamma v_2,\\ g_3=u_3+{\eta }^{-1}\gamma v_3+{\eta }^{-1}{u}_2{v}_2+{\eta }^{-2}\gamma v_2^2.\end{array}$ (39)

令$\rho \left(\theta ,d\right)$ 是方程(38)以初始值$\rho \left(0,d\right)=d$ 的解。方程两边同时对$\theta$ 积分，通过计算可以得到，

$\frac{\rho \left(\theta ,d\right)-d}{d\rho \left(\theta ,d\right)}=\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }\left[\frac{\gamma }{\rho }+{\displaystyle \sum _{k=2}^3{\rho }^{k-2}{g}_k+o\left(\rho \right)}\right]\text{d}s,}$

这意味着$\rho \left(\theta ,d\right)$ 对于$d$ 总有如下泰勒展开式，

$\rho \left(\theta ,d\right)=d+{\kappa }_2\left(\theta ,d\right)d^2+{\kappa }_3\left(\theta ,d\right)d^3+o\left(d^3\right).$ (40)

联立(38)和(40)，可以得到，

$\frac{\text{d}\rho }{\text{d}\theta }=\frac{\gamma }{\eta }d+\left(\frac{2\gamma }{\eta }{\kappa }_2+\frac{g_2}{\eta }\right)d^2+\left[\frac{\gamma }{\eta }\left({\kappa }_2^2+2{\kappa }_3\right)+2{\kappa }_2\frac{g_2}{\eta }+\frac{g_3}{\eta }\right]d^3+o\left(d^3\right).$ (41)

将(41)两端关于$\theta$ 积分，则有，

$\begin{array}{c}\rho \left(\theta ,d\right)=d+\frac{\gamma }{\eta }d\theta +\frac{d^2}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\text{d}s}+\frac{d^3}{\eta }{\displaystyle {\int }_0^{\theta }\left(2{\kappa }_2{g}_2+g_3\right)\text{d}s}\\ +d^2{\displaystyle {\int }_0^{\theta }\frac{2\gamma }{\eta }{\kappa }_2\text{d}s}+d^3{\displaystyle {\int }_0^{\theta }\frac{\gamma }{\eta }\left({\kappa }_2^2+2{\kappa }_3\right)\text{d}s}+o\left(d^3\right).\end{array}$ (42)

对比(40)和(42)，由对应系数相等可以得到，

${\kappa }_2=\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\text{d}s},{\kappa }_3=\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }\left(2{\kappa }_2{g}_2+g_3\right)\text{d}s}.$ (43)

将(38)式两端关于$\theta$ 从0到$2\pi$ 积分，并联立(43)可以得到，

$\begin{array}{c}\frac{\rho \left(2\pi ,d\right)-\rho \left(0,d\right)}{\rho \left(2\pi ,d\right)}=\frac{\gamma \tilde{\rho }}{\eta }+\frac{d}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_2\text{d}\theta }+\frac{d^2}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_3\text{d}\theta }+\frac{d^3}{\eta }{\displaystyle {\int }_0^{2\pi }\left(g_3\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\text{d}s}\right)\text{d}\theta }\\ +\frac{2d^4}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_3\left(\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\left(\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\text{d}s}\right)\text{d}s}\right)\text{d}\theta }\\ +\frac{d^4}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_3\left(\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_3\text{d}s}\right)\text{d}\theta }+o\left(d^4\right),\end{array}$ (44)

其中，$\tilde{\rho }=2\pi +o\left(d\right)$ 。进一步，通过直接计算可得，

$\begin{array}{l}{\displaystyle {\int }_0^{2\pi }{g}_2\text{d}\theta }=0,\\ {\displaystyle {\int }_0^{2\pi }{g}_3\left({\displaystyle {\int }_0^{\theta }{g}_2\text{d}s}\right)\text{d}\theta }=\frac{2}{3}{a}_{02}^2{\displaystyle {\int }_0^{2\pi }{g}_3\text{d}\theta }+o\left(d\right),\\ {\displaystyle {\int }_0^{2\pi }{g}_3\left({\displaystyle {\int }_0^{\theta }{g}_3\text{d}s}\right)\text{d}\theta }=0.\end{array}$ (45)

因此，(44)可以等价表示为，

$\frac{\rho \left(2\pi ,d\right)-\rho \left(0,d\right)}{\rho \left(2\pi ,d\right)}=\frac{\gamma \tilde{\rho }}{\eta }+{\xi }_2{d}^2+{\xi }_3{d}^3+{\xi }_4{d}^4+o\left(d^4\right),$ (46)

其中，

$\begin{array}{l}{\xi }_2=\frac{1}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_3\text{d}\theta },{\xi }_3=\frac{2}{3}{a}_{02}^2{\xi }_2,\\ {\xi }_4=\frac{2}{\eta }{\displaystyle {\int }_0^{2\pi }{g}_3\left(\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\left(\frac{1}{\eta }{\displaystyle {\int }_0^{\theta }{g}_2\text{d}s}\right)\text{d}s}\right)\text{d}\theta }.\end{array}$ (47)

由于约化方程(36)分歧出的周期解与下述方程的正根$d_0$ 密切相关，

$\frac{\gamma \tilde{\rho }}{\eta }+{\xi }_2{d}^2+{\xi }_3{d}^3+{\xi }_4{d}^4+o\left(d^4\right)=0.$ (48)

下面，我们对(48)展开分析。当$d$ 充分接近$d_0$ 时，如果

$\rho \left(2\pi ,d\right)-\rho \left(0,d\right)\{\begin{array}{l}<0,d>d_0,\\ >0,d<d_0.\end{array}$ (49)

则从$d_0$ 处分歧出的周期轨道是稳定的。如果

$\rho \left(2\pi ,d\right)-\rho \left(0,d\right)\{\begin{array}{l}>0,d>d_0,\\ <0,d<d_0.\end{array}$ (50)

则从$d_0$ 处分歧出的周期轨道是不稳定的。

为了确定分歧解的稳定性和跃迁的类型，令

$N=\frac{\gamma \tilde{\rho }}{\eta }+{\xi }_2{d}^2+{\xi }_3{d}^3+{\xi }_4{d}^4+o\left(d^4\right).$ (51)

已知方程(28)的跃迁及其局部拓扑结构完全由(36)确定，并且对于任意很小的正数$d$ ，${\omega }^{+}$ 在临界值$c_2$ 的稳定性完全由N的符号决定，称N为跃迁数。当$N<0$ 时，系统的动态跃迁是连续型的；而当$N>0$ 时，系统的动态跃迁是跳跃型的。

在临界值$c_2$ 处，$\gamma =0$ ，此时N的符号完全由${\xi }_2\left(c_2\right)$ 决定。当${\xi }_2\left(c_2\right)<0$ 时，系统(28)发生连续型跃迁。当${\xi }_2\left(c_2\right)>0$ 时，系统(28)发生跳跃型跃迁。如果${\xi }_2\left(c_2\right)=0$ ，则${\xi }_3=0$ ，我们用${\xi }_4\left(c_2\right)$ 来判断跃迁类型。

对于系统(28)的二次跃迁，我们有下述结论。

定理3.2当${\xi }_2\left(c_2\right)>0$ 或${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)>0$ ，关于系统(28)有如下结论成立：

1) 系统(28)在${\omega }^{+}$ 处发生跳跃型跃迁，并在$c<c_2$ 一侧分歧出一个不稳定周期轨道${\Gamma }_1$ ；

2) 存在一个亚临界跃迁数$c^{\ast }\left(c_1<c^{\ast }\le c_2\right)$ ，在该跃迁数处存在一个周期轨道奇异分离。该系统从${\omega }^{+}$ 处分歧出非零吸引子${\Gamma }_2$ 。而当$c>c_2$ 时，系统(28)不存在从${\omega }^{+}$ 分歧的周期解。

证明3.2在此假设下，很容易得到对于任意很小的数$d$ ，(51)中定义的N在$c=c_2$ 处为正，这意味着${\omega }^{+}$ 是不稳定的，并且在$c=c_2$ 处发生跳跃型跃迁。且当${\xi }_2\left(c_2\right)>0$ 或${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)>0$ ，(48)在$c<c_2$ 一侧具有下述正实根，

${\Gamma }_1={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_2}\right)}^{\frac{1}{2}},$ 当${\xi }_2\left(c_2\right)>0.$

或者

${\Gamma }_1={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_4}\right)}^{\frac{1}{4}},$ 当${\xi }_2\left(c_2\right)=0,{\xi }_4\left(c_2\right)>0.$

结合(50)可知分歧出的${\Gamma }_1$ 是不稳定的。

基于跃迁理论可以分析得到，存在一个亚临界跃迁数$c^{\ast }\left(c_1<c^{\ast }\le c_2\right)$ ，当$c=c^{\ast }$ 时，系统从${\omega }^{+}$ 处分歧出一个非零吸引子。当$c>c_2$ 时，方程(51)无根，即系统在$c>c_2$ 一侧不存在从${\omega }^{+}$ 分歧的周期解。

定理3.3如果${\xi }_2\left(c_2\right)<0$ 或${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)<0$ ，系统(28)在$c=c_2$ 处发生连续型跃迁。并在$c>c_2$ 一侧从稳态解${\omega }^{+}$ 处分歧出稳定的周期轨道$\Gamma$ 。此外，周期轨道可以近似表达为，

1) 当${\xi }_2\left(c_2\right)<0$ 时，

$\omega \left(t\right)={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_2}\right)}^{\frac{1}{2}}\left[\cos \left(\eta t\right)e_1+\sin \left(\eta t\right)e_2\right]+o\left(\left|\gamma \left(c\right)\right|\right).$

2) 当${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)<0$ 时，

$\omega \left(t\right)={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_4}\right)}^{\frac{1}{4}}\left[\cos \left(\eta t\right)e_1+\sin \left(\eta t\right)e_2\right]+o\left(\left|\gamma \left(c\right)\right|\right).$

证明3.3当${\xi }_2\left(c_2\right)<0$ 或者${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)<0$ ，很容易得到对于任意很小的数$d$ ，(51)中定义的$N$ 在$c=c_2$ 处为负，这意味着${\omega }^{+}$ 是稳定的，并且从$\left({\omega }^{+},c_2\right)$ 处发生的动态跃迁是连续型的。此外，不难得到，当$c>c_2$ 时，(48)只有一个正根。

当${\xi }_2\left(c_2\right)<0$ ，$\Gamma ={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_2}\right)}^{\frac{1}{2}}$ ，或者当${\xi }_2\left(c_2\right)=0$ ，${\xi }_4\left(c_2\right)<0$ ，$\Gamma ={\left(\frac{-\tilde{\rho }\gamma }{\eta {\xi }_4}\right)}^{\frac{1}{4}}$ 。

结合(49)可以得到分歧出的周期轨道$\Gamma$ 是稳定的。

4. 总结

本文基于谱理论和动态跃迁理论对三维Chen系统的动力学行为进行了系统分析。理论研究表明，随着参数变化，系统出现两次动态跃迁。

本文借助谱理论得到了Chen系统发生一次跃迁的临界控制参数$c_1$ 和二次跃迁的临界控制参数$c_2$ ，以及相应的PES条件。借助动态跃迁理论和中心流形约化，得到了该系统发生动态跃迁的相关结论。首先，当参数$c$ 穿越第一个临界值$c_1$ 时，系统发生连续型跃迁，从原有的常数稳态解处分歧出两个新的渐近稳定的奇点。其次，随着参数$c$ 的增大，当参数$c$ 穿越第二个临界值$c_2$ 时，Chen系统出现二次跃迁行为。其跃迁类型由具体的参数条件决定：在特定条件下发生跳跃型跃迁，系统分歧出不稳定的周期轨道；而在相反的条件下则发生连续型跃迁，分歧出渐近稳定的周期轨道。

值得注意的是，本文的“连续型”与“跳跃型”跃迁，并非纯粹的数学抽象。在物理世界中，它们通常对应着截然不同的失稳机制和现象。

一次跃迁(在$c_1$ 处，连续型)，系统从唯一平衡态平滑地分歧出两个稳定的新平衡点。这在物理上可类比于一个对称破缺的软失稳过程。例如，在磁体模型中，当温度降低通过居里点时，系统从无磁状态连续地过渡到两个等价的、具有相反磁化的稳定态之一。系统演化路径没有突变，初值的小扰动仅决定最终落入哪一个分支。

二次跃迁有两种情形，当发生跳跃型跃迁时：系统从平衡态分歧出不稳定的周期轨道。这强烈暗示了系统存在一个亚临界的Hopf分岔。在物理上，这对应于一种硬激发现象。系统不会自发地小幅振荡并平滑增长至大幅周期运动；相反，要激发出一个稳定的大幅振荡(可能存在于更远的吸引子上)，需要一个有限幅度的强扰动。这种行为常见于非线性电子电路、光学双稳态系统以及流体对流中，并伴随滞后效应，即参数正向和反向变化时，系统状态变化的路径不同。当发生连续型跃迁时：系统分歧出稳定的周期轨道。这对应一个超临界的Hopf分岔。在物理上，这意味着当参数越过临界值，平衡态失稳后，系统会自发地产生小幅振荡，并逐渐增长到一个稳定的极限环振幅。这是许多自持振荡(如化学振荡反应、神经元电活动、机械颤振)的经典产生机制。

基金项目

本项工作由重庆市自然科学基金面上项目(No. CSTB2022NSCQ-MSX1609)和重庆交通大学科研启动项目(No. 2020023034)资助完成。

NOTES

^*通讯作者。