DOI:10.12677/aam.2024.136265,PDF,HTML,XML,被引量下载: 79浏览: 165科研立项经费支持
作者:张彦超,王鑫乐,赵安旭,赵楠楠^*：长安大学理学院，陕西 西安
关键词:混沌R?ssler系统；临界耦合强度；最大Lyapunov指数；同步；Chaotic R?ssler System；Critical Coupling Strength；Maximum Lyapunov Exponent；Synchronization

文章引用：张彦超, 王鑫乐, 赵安旭, 赵楠楠. 网络拓扑对混沌系统同步能力的影响[J]. 应用数学进展, 2024, 13(6): 2761-2770. https://doi.org/10.12677/aam.2024.136265

1. 引言

自20世纪60年代以来，非线性科学领域蓬勃发展。在非线性动力学崭露头角的同时，同步理论也逐渐成形。1673年，Huygens观察到了两个悬挂在同一横梁上的弱耦合摆钟同步运动的现象，这标志着同步概念的最早提出。随后的研究主要集中在耦合周期系统上。直到1990年，Pecora和Carroll证明了两个耦合的混沌系统可以实现同步，引发了混沌同步研究的热潮[1][2]。此后，研究者发现了各种类型的耦合混沌系统同步状态，包括完全同步(complete synchronization)、相位同步(phase synchronization)、滞后同步(lag synchronization)、期望同步(anticipating synchronization)、射影同步(projective synchronization)、广义同步(generalized synchronization)、阵发性滞后同步(intermittent lag synchronization)、弱相位同步(imperfect phase synchronization)和几乎完全同步(almost complete synchronization)等。因此，同步研究已经从周期振子的锁相同步扩展到混沌系统的同步，并从两个或三个耦合系统的同步发展到了复杂网络系统的同步[3]-[6]。在复杂网络中，节点的同步指的是网络中各个节点的动态行为呈现出一致性[7]。同步现象在自然界中广泛存在，包括物理、生物、化学以及社会系统，在这些系统中都可以观察到，并且具有丰富的应用场景。例如，Josephson-Junction阵列[8]、半导体激光器阵列[9]、萤火虫的闪烁[10]、心脏起搏细胞[11]、帕金森疾病[12]等都涉及到同步现象，并且还有其他一些应用[13][14]。

1998年，Pecora等人提出了著名的主稳定函数法来解决线性耦合振子网络的同步稳定性问题[15]。该方法适用于任何线性耦合网络，并至今仍然是研究网络同步的重要方法之一。2004年，Timme等人提出了网络同步速度的数学表达式，并指出网络的连通性限制了同步速度[16]。2009年，Yan等人的研究指出时滞网络同步的临界值与最大拉普拉斯特征值成反比，而非时滞网络的收敛速度与代数连通度(第二最小拉普拉斯特征值)成正比。这一研究表明了网络同步对网络拓扑结构的依赖性[17]。2011年，舒永乐等人分析了一个五维混沌系统的最终界，为混沌同步研究提供了理论基础[18]。2014年，行鸿彦等人得出了统一变形混沌系统的全局指数吸引集，并将结果应用到了统一混沌系统的异构同步中[19]。2017年，刘莹等人采用线性反馈控制方法实现了强迫布鲁塞尔振子的全局指数同步[20]。2019年，舒睿等人讨论了星型网络中耦合Kuramoto振子的同步优化问题，分别就单个同步集团和多个同步集团两种情形，得到了最小和最大同步临界耦合强度[21]。2024年，陈松等人对一个复杂混沌系统进行了理论分析，并设计了合适的控制器实现了，两个相同混沌系统的全局指数同步[22]。

本文旨在探讨网络拓扑对混沌系统同步能力的影响问题。我们将以Rössler系统为例，通过数值模拟方法观察单个系统的混沌行为，然后构建规范网络模型，研究多个节点相互耦合的Rössler系统在不同网络拓扑结构下的同步现象。通过主稳定函数法和同步误差函数的应用，我们将计算出系统同步发生的最低耦合强度阈值，并进一步比较不同网络拓扑结构下的临界耦合强度，以评估网络拓扑对混沌系统同步能力的影响。本文的研究不仅有助于深入理解混沌系统同步现象的形成机制，还为设计和优化复杂网络系统提供了重要参考。

2. 单个Rössler模型

单个罗斯勒(Rössler)系统模型

$\{\begin{array}{l}\stackrel{\dot{}}{x}=-\left(y+z\right)\\ \stackrel{\dot{}}{y}=x+ay\\ \stackrel{\dot{}}{z}=b+z\left(x-c\right)\end{array}$(2.1)

其中当 $a=b=0.2$， $c=9$，Rössler系统是混沌的。吸引子相图描绘的是系统的解曲线在某一时间段的运动轨迹，通过观察不同维度上混沌吸引子的相图，可以初步判断系统所处的状态。选取初值 $\left(x_0,y_0,z_0\right)=\left(1,0,0\right)$，利用龙格–库塔法在MATLAB内，可求解系统混沌吸引子的二维和三维相图，如图1和图2所示。

Figure 1.Two-dimensional plot of system chaotic attractor

图1.系统混沌吸引子的二维图

Figure 2.Three-dimensional plot of system chaotic attractor

图2.系统混沌吸引子的三维图

3. 网络化Rössler模型

多个振子相互耦合的Rössler模型

${\stackrel{\dot{}}{x}}_i=f\left(x_i\right)-\sigma {\displaystyle \underset{j=1}{\overset{N}{\sum }}{g}_{ij}H{x}_j}$(3.1)

其中 $i=1,2,\cdots ,N$。 $x_i\in {ℝ}^n$是第i个振子的振荡状态， $g_{ij}$是邻接矩阵的元素，且 $g_{ij}=g_{ji}=-1$当第i个振子和第j个相连接，其余的都为零，H是内部耦合矩阵。

全局同步的条件： $x_1\left(t\right)=x_2\left(t\right)=\cdots =x_N\left(t\right)=s\left(t\right)$， $\stackrel{\dot{}}{s}\left(t\right)=F\left(s\left(t\right)\right)$。将(3.1)改写为：

$\stackrel{\dot{}}{x}=F\left(x\right)-\sigma G\otimes I_n\cdot H\left(x\right)$(3.2)

其中 $x={\left[x_1^{\text{T}},x_2^{\text{T}},\cdots ,x_N^{\text{T}}\right]}^{\text{T}}$， $I_n$表示一个 $n\times n$的单位阵(n是状态维度)。当 $H\left(x\right)$为线性的时候，退化为Hx(H是一个常数矩阵)。

$F\left(x\right)={\left[\begin{array}{c}f\left(x_1\right)\\ f\left(x_2\right)\\ ⋮\\ f\left(x_N\right)\end{array}\right]}_{\left(N\cdot n\right)\times 1},H\left(x\right)={\left[\begin{array}{c}h\left(x_1\right)\\ h\left(x_2\right)\\ ⋮\\ h\left(x_N\right)\end{array}\right]}_{\left(N\cdot n\right)\times 1}$(3.3)

令 $\eta ={\left[\begin{array}{cccc}{\eta }_1^{\text{T}}& {\eta }_2^{\text{T}}& \cdots & {\eta }_N^{\text{T}}\end{array}\right]}^{\text{T}}$， ${\eta }_i\left(t\right)=x_i\left(t\right)-s\left(t\right)$，则有：

$\stackrel{\dot{}}{\eta }=\left(I_N\otimes Df\right)\eta -\sigma \left(G\otimes Dh\right)\eta$(3.4)

其中 $Df={\frac{\partial f\left(x_i\right)}{\partial x_i}|}_{x_i=s\left(t\right)}$， $Dh={\frac{\partial h\left(x_i\right)}{\partial x_i}|}_{x_i=s\left(t\right)}$，做变换 $\xi =\left(T^{-1}\otimes I_n\right)\eta$，T满足 $T^{-1}GT=\Gamma$， $\Gamma$是一个对角阵。即 ${\stackrel{\dot{}}{\xi }}_i=\left[Df-\sigma {\gamma }_{i}Dh\right]{\xi }_i$，其中 ${\gamma }_i$是G的特征值。得到主稳定方程：

$\stackrel{\dot{}}{\xi }=\left[Df-\alpha Dh\right]\xi$(3.5)

全局同步的分析步骤：

1) 独立于具体网络的拓扑结构，通过研究方程(3.5)的最大Lyapunov指数 ${\lambda }_{\max }$作为 $\alpha$的函数，从而得到主稳定函数MSF，即 ${\lambda }_{\max }={\lambda }_{\max }\left(\alpha \right)$。

2) 在点 $a=\sigma {\gamma }_i,i=2,3,\cdots ,N$( ${\gamma }_i$是由G导出的特征值)的符号，如果所有点都稳定，那么在给定的耦合强度下稳定。

${\lambda }_{\max }$对a的函数依赖性导致了三种不同的情况，如图3左图所示。第一种情况(类型I)使得 ${\lambda }_{\max }$是恒正，因此，网络节点永远不可能同步。在第二种情况下(类型II)， ${\lambda }_{\max }$在足够强的耦合作用下总会同步，即 ${\alpha }_1<\sigma {\gamma }_i$， $i=2,3,\cdots ,N$。在第三种情况(类型III)中， ${\lambda }_{\max }$仅在a的封闭区间 $\left[{\alpha }_1,{\alpha }_2\right]$内为负，即 ${\alpha }_1<\sigma {\gamma }_i<{\alpha }_2$， $i=2,3,\cdots ,N$时，网络节点才能同步。

Rössler系统的最大Lyapunov指数的计算，对方程于(3.5)，其中

$Df=\left[\begin{array}{ccc}0& -1& -1\\ 1& b& 0\\ z& 0& x-c\end{array}\right],Dh=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],\xi =\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$(3.6)

利用G. Benettin计算方法，可以得到II型曲线( ${\alpha }_1=0.2$)，如图3右图所示

Figure 3.MSF curve plot and maximum Lyapunov exponent plot of the Rössler system

图3.MSF曲线图和Rössler系统的最大Lyapunov指数图

4. 数值求解

下面我们以具有N(N= 10)个节点相互线性耦合的Rössler系统为例，网络结构如图4规范网络(example)所示。

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_i=-\left(y_i+z_i\right)\\ {\stackrel{\dot{}}{y}}_i=x_i+\alpha y_i-\sigma {\displaystyle \underset{j=1}{\overset{N}{\sum }}{g}_{ij}{y}_j}\\ {\stackrel{\dot{}}{z}}_i=b+z_i\left(x_i-c\right))\end{array}$(4.1)

拉普拉斯矩阵(Laplace) $L=D-A$。其中A是邻接矩阵，当i与j相互连接时， $a_{ij}=a_{ji}=1$，否则 $a_{ij}=0$。

其中D是度矩阵，元素 $d_{ii}={\displaystyle \underset{j}{\overset{N}{\sum }}{a}_{ij}}$，其余元素为零。因为拉普拉斯矩阵是零行和矩阵，那么 ${\lambda }_{\min }=0$。对其所有特征值按照递增排序，则 $0={\lambda }_1<{\lambda }_2<\cdots <{\lambda }_N$。

记 ${\bar{\sigma }}_0={\alpha }_1/{\lambda }_2$为系统发生同步的理论临界耦合强度值，当 $\sigma >{\bar{\sigma }}_0$时，系统始终发生全局同步现象。我们可以知道example网络的理论临界耦合强度 ${\bar{\sigma }}_0=1.0379$。

Figure 4.Standard network diagram

图4.规范网络图

接着我们定义同步误差函数 $e\left(t\right)$，对应(3.1)式，我们以 $x_i$表示第i个振子t时刻的振动状态，且 $x_i={\left[x_1,x_2,\cdots ,x_N\right]}^{\text{T}}$，则有：

$e\left(t\right)={\left(\frac{1}{N\left(N-1\right)}{\displaystyle \underset{ij}{\sum }{\Vert x_i\left(t\right)-x_j\left(t\right)\Vert }^2}\right)}^{\frac{1}{2}}$(4.2)

当 $e\left(t\right)\to 0$时，系统处于完全同步态。通过修改(4.1)式中的参数 $\sigma$，观察最后时刻 $e\left(t\right)$的值，可以找到数值上临界耦合强度的值，参照图5中的右图。

Figure 5.Left—synchronization error time plot, right—synchronization error coupling strength plot

图5.左图为同步误差时间图，右图为同步误差耦合强度图

如图5所示，我们可以观察到Rössler系统发生同步，所需要的临界耦合强度为1.06与理论值绝对误差为0.0221。

5. 不同网络同步能力比较

下面我们计算，Rössler系统在其余四种网络结构下的临界耦合强度，包含数值解和理论解。其详细过程雷同上述规范网络过程，不再赘述。

ER随机网络，构造方法如下， $p=0.03$， $N=10$， $\sigma =1.2$，p是每对节点的连接概率，一个节点具有d条边的概率为

$p_d=\left(\begin{array}{c}N-1\\ d\end{array}\right)p^d{\left(1-p\right)}^{N-1-d}$(5.1)

平均每个节点的边为 $\langle d\rangle =p\left(N-1\right)$条。

ER网络下，如图6左图所示，可以发现网络化系统实现同步。数值上临界耦合强度为0.32，理论上临界耦合强度为0.329，绝对误差为0.029。

Figure 6.Left—synchronization error time plot, right—synchronization error coupling strength plot

图6.左图为同步误差时间图，右图为同步误差耦合强度图

小世界(WS)网络，构造方法如下， $N=10$， $p=3$， $\rho =0.1$，p为邻近数， $\rho$为重新布线概率，首先生成p邻环，每个节点由2p个邻居，再以 $\rho$的概率重新链接.

在小世界网络中，如图7所示，可以发现网络化系统实现同步。数值上临界耦合强度为0.05，理论上临界耦合强度为0.0485，绝对误差为0.0015。

BA网络，构造方法如下， $N=10$， $m=2$， $n_0=3$， $\sigma =1.2$，以 $n_0$个相互链接的节点开始，每新增一个节点，与m个旧节点链接( $m<n_0$)概率为

${\Pi }_i=\frac{d_i}{{\displaystyle \underset{j}{\sum }{d}_j}}$(5.2)

其中 $j=1,2,\cdots ,n$，n为旧节点数。

在BA网络中，如图8所示，可以发现网络化系统实现同步. 数值上临界耦合强度为0.16，理论上临界耦合强度为0.1694，绝对误差为0.0094。

Figure 7.Left—synchronization error time plot, right—synchronization error coupling strength plot

图7.左图为同步误差时间图，右图为同步误差耦合强度图

Figure 8.Left-synchronization error time plot, right-synchronization error coupling strength plot

图8.左图为同步误差时间图，右图为同步误差耦合强度图

Figure 9.Left—synchronization error time plot, Right—synchronization error coupling strength plot

图9.左图为同步误差时间图，右图为同步误差耦合强度图

环形网络下，如图9所示，可以发现网络化系统实现同步. 数值上临界耦合强度为0.54，理论上临界耦合强度为0.5236，绝对误差为0.0164。

最后我们对以上五种网络的同步能力进行比较

根据图10我们可以得到小世界(WS)网络实现同步所需的临界耦合强度最小，即同步能力最强，而规则网络(Example)所需的临界耦合强度最大，即同步能力最弱。猜测这种能力与网络的平均度有关，平均度定义为：

$\langle d\rangle =\frac{{\displaystyle \sum diag\left(L\right)}}{N}$(5.3)

其中L为网络图的拉普拉斯矩阵， $diag\left(L\right)$表示取L对角线上的元素，N为节点数。

Figure 10.Comparison of theoretical critical coupling strengths and numerical solutions across different networks

图10.不同网络下，临界耦合强度，数值解理论界比较图

通过表1可以观察到，随着平均度 $\langle d\rangle$的增加，网络的同步能力基本也在增大，同步能力的增强与平均度呈正相关性。

Table 1.Comparison of theoretical critical coupling strengths and average degrees across different networks

表1.不同网络下，理论临界耦合强度，平均度比较图

Network	WS	BA	ER	Circle	Example
平均度 $\langle d\rangle$	6	3.4	3	2	2.2
理论临界耦合强度 ${\bar{\sigma }}_0$	0.0485	0.1694	0.329	0.5236	1.0379

6. 结论

本文讨论了网络拓扑对混沌系统同步能力的影响问题。首先通过龙格–库塔法求解了单个Rössler系统模型，观察到了混沌现象。其次以规范网络为例，多个节点相互耦合的Rössler系统为模型，观察到同步现象的发生，同时通过主稳定函数法和同步误差函数分别求解出系统发生同步现象，所需要的最低耦合强度阈值。最后通过比较不同网络拓扑结构下，临界耦合强度的大小，判断系统同步能力的强弱。并猜测网络拓扑对混沌系统的影响与网络的平均度密切相关，平均度越大，网络同步能力越强，所需要的临界耦合强度越小。这与直觉也密切符合，平均每个节点所影响的节点数越多，混沌系统更容易发生同步现象。

基金项目

由大学生创新创业训练计划项目资金支持(项目编号：S202310710314)。

NOTES

^*通讯作者。

参考文献

[1]	Pecora, L.M. and Carroll, T.L. (1991) Driving Systems with Chaotic Signals.Physical Review A, 44, 2374-2383. https://doi.org/10.1103/physreva.44.2374
[2]	Pecora, L.M. and Carroll, T.L. (1990) Synchronization in Chaotic Systems.Physical Review Letters, 64, 821-824. https://doi.org/10.1103/physrevlett.64.821
[3]	Barahona, M. and Pecora, L.M. (2002) Synchronization in Small-World Systems.Physical Review Letters, 89, Article 054101. https://doi.org/10.1103/physrevlett.89.054101
[4]	Masuda, N. and Aihara, K. (2004) Global and Local Synchrony of Coupled Neurons in Small-World Networks.Biological Cybernetics, 90, 302-309. https://doi.org/10.1007/s00422-004-0471-9
[5]	Li, X. and Chen, G. (2003) Synchronization and Desynchronization of Complex Dynamical Networks: An Engineering Viewpoint.IEEE Transactions on Circuits and Systems I:Fundamental Theory and Applications, 50, 1381-1390. https://doi.org/10.1109/tcsi.2003.818611
[6]	Belykh, I., de Lange, E. and Hasler, M. (2005) Synchronization of Bursting Neurons: What Matters in the Network Topology.Physical Review Letters, 94, Article 188101. https://doi.org/10.1103/physrevlett.94.188101
[7]	Ditto, W. (2002) Keeping in Sync.Nature, 415, 736-737. https://doi.org/10.1038/415736b
[8]	Wiesenfeld, K., Colet, P. and Strogatz, S. (1998) Frequency Locking in Josephson Arrays: Connection with the Kuramoto Model.Physical Review E, 57, 1563-1569. https://doi.org/10.1103/physreve.57.1563
[9]	Argyris, A., Syvridis, D., Larger, L., Annovazzi-Lodi, V., Colet, P., Fischer, I.,et al. (2005) Chaos-Based Communications at High Bit Rates Using Commercial Fibre-Optic Links.Nature, 438, 343-346. https://doi.org/10.1038/nature04275
[10]	Ermentrout, B. (1991) An Adaptive Model for Synchrony in the FireflyPteroptyxmalaccae.Journal of Mathematical Biology, 29, 571-585. https://doi.org/10.1007/bf00164052
[11]	Vinogradova, T.M., Lyashkov, A.E., Zhu, W., Ruknudin, A.M., Sirenko, S., Yang, D.,et al. (2006) High Basal Protein Kinase A—Dependent Phosphorylation Drives Rhythmic Internal Ca2+Store Oscillations and Spontaneous Beating of Cardiac Pacemaker Cells.Circulation Research, 98, 505-514. https://doi.org/10.1161/01.res.0000204575.94040.d1
[12]	Stam, C.J. (2005) Nonlinear Dynamical Analysis of EEG and MEG: Review of an Emerging Field.Clinical Neurophysiology, 116, 2266-2301. https://doi.org/10.1016/j.clinph.2005.06.011
[13]	MacLeod, K. and Laurent, G. (1996) Distinct Mechanisms for Synchronization and Temporal Patterning of Odor-Encoding Neural Assemblies.Science, 274, 976-979. https://doi.org/10.1126/science.274.5289.976
[14]	Prakash, M. and Gershenfeld, N. (2007) Microfluidic Bubble Logic.Science, 315, 832-835. https://doi.org/10.1126/science.1136907
[15]	Pecora, L.M. and Carroll, T.L. (1998) Master Stability Functions for Synchronized Coupled Systems.Physical Review Letters, 80, 2109-2112. https://doi.org/10.1103/physrevlett.80.2109
[16]	Timme, M., Wolf, F. and Geisel, T. (2004) Topological Speed Limits to Network Synchronization.Physical Review Letters, 92, Article 074101. https://doi.org/10.1103/physrevlett.92.074101
[17]	Yan, G., Chen, G., Lü, J. and Fu, Z. (2009) Synchronization Performance of Complex Oscillator Networks.Physical Review E, 80, Article 056116. https://doi.org/10.1103/physreve.80.056116
[18]	舒永录, 张付臣, 杨洪亮. 一个新的多维超混沌系统及其性质研究[J]. 四川大学学报(自然科学版), 2011, 48(4): 857-864.
[19]	行鸿彦, 冒海微, 徐伟. 基于全局指数吸引集的统一变形混沌系统同步[J]. 信息与控制, 2014, 43(4): 385-391.
[20]	刘莹, 王贺元, 陈荟颖. 强迫布鲁塞尔振子动力学行为和全局指数同步的数值仿真[J]. 动力学与控制学报, 2017, 15(5): 423-429.
[21]	舒睿, 陈伟, 肖井华. 多个耦合星型网络的同步优化[J]. 物理学报, 2019, 68(18): 65-74.
[22]	陈松, 张付臣, 肖敏. 一个复杂混沌系统的分析与同步控制[J]. 系统科学与数学, 2024, 44(5): 1311-1323.