0 引言

【研究意义】近年来，复杂网络的发展速度十分迅猛，复杂网络也逐渐成为物理学、计算机科学、数学等多领域交叉的热点研究问题，对复杂网络中同步现象的研究对社会实践有着重要的影响与指导意义^[1-2]。复杂网络上的耦合相位振子的同步化对研究网络中振子出现自发的集体行为现象有重大意义。【前人研究进展】自Kuramoto模型提出以来，该模型已有广泛的应用，如研究一阶与二阶相互作用项对同步的影响^[3]，或者讨论基于该模型在极限环振子中加入时间延迟因素后的影响^[4]，又或者讨论噪声^[5]或错位^[6]对爆炸式同步的影响等等。在众多的工作中，Gómez-Gardeóes等^[7]从微观层面研究在复杂网络中爆炸式同步是如何产生的，并在星形网络中给出分析研究，重现无标度网络中的结果；Sakaguchi等^[8]在规则网络上分别讨论相位振子之间的错位和相互拖拽对系统同步的影响。【本研究切入点】在现实生活中，仅有Kuramoto模型是不能反应更复杂的系统的。在无标度网络中加入相互拖拽因子，可以针对某些更复杂的问题给出理论依据；而对在星形网络加入相互拖拽的研究，则可更深刻地理解无标度网络中相互拖拽的影响。【拟解决的关键问题】在原有Kuramoto模型基础上，加入相互拖拽因子，并进行数值仿真，分析在无标度网络中相互拖拽对同步的影响。由于星形网络有无标度网络的主要性质^[7]，故在考虑减少分析困难的基础上，对星形网络中加入相互拖拽因子对同步的影响做理论分析，以求更深入理解无标度网络下的情形。

1 模型的建立

1970年Kuramot模型被提出^[9-11]。在一个无权重、无向的网络中，每一个振子的相位定义为θ_i(t) (i=1, 2, …, N)，并且每个振子θ_i(t)随时间演变的方程满足：

$ {\dot \theta _i} = {\omega _i} + \lambda \sum\limits_{j = 1}^N {{A_{ij}}{\rm{sin}}\left( {{\theta _j}-{\theta _i}} \right)}, \;\;i = 1, 2, \cdots, N, $

(1)

其中，ω_i为本征频率，A_ij为邻接矩阵(当振子的第i节点与第j节点连接时，A_ij=1;当振子的第i节点与第j节点无连接时，A_ij=0)，λ代表振子i与j间的耦合强度大小。

原始的Kuramot模型假定所有的振子都是连接在一起的，即对于所有的i与j，都有A_ij=1，并发现当振子的本征频率ω_i等于振子的度k_i时，会发生爆炸式同步^[7]。同时使用序参量r来量化同步的程度，而序参量r被定义为

$ r\left( t \right){e^{i\psi \left( t \right)}} = \frac{1}{N}\sum\limits_{j = 1}^N {{e^{i{\theta _j}}}}, $

(2)

序参量r满足0≤r≤1，当r=1时，网络中的振子达到完全同步；当r=0时，则网络中振子全体展现出不一致的状态^[10-11]。

在无标度网络上，考虑Sakaguchi模型^[8]，加入拖拽因子sin β，讨论sin β对同步的影响。该模型为

$ \begin{array}{l} \;\;\;\;\;\;\;\;{{\dot \theta }_i} = {\omega _i} + \lambda \sum\limits_{j = 1}^N {{A_{ij}}{\rm{sin}}\left( {\sin \left( {{\theta _j}-{\theta _i}} \right) + \sin \;\beta } \right)}, \\ i = 1, 2, \cdots, N。\end{array} $

(3)

根据三角函数的性质，只需讨论$\beta \in \left( {-\frac{\pi }{2}, \frac{\pi }{2}} \right)$的范围即可。设定振子ω_i等于节点的度k_i，度分布满足幂律关系p(k)∝k^-ν, v=3。

为探究sin β对同步的影响，可以采用有100个振子的无标度网络，同时为从微观层次对同步进程进行研究，定义平均有效频率为

$ \omega _i^{eff} = \frac{1}{T}\int\limits_\tau ^{\tau + T} {{{\dot \theta }_i}\left( \tau \right)d\tau 。} $

(4)

将时间设为T=10 000，在无标度网络上，采用四阶龙格库塔方法，设置时间步长为dt=0.01，求平均有效频率＜ω_i＞和序参量r随耦合强度λ的变化关系。

2 模型的数值模拟

对方程(3) 进行数值模拟，可以得到平均有效频率＜ω_i＞和耦合强度λ的变化关系(图 1)以及序参量r随耦合强度λ的变化关系(图 2)。从图 1与图 2可以看出，对比β为π的情况，当sin β的值为正值，sin β对同步起到提前的作用，并且sin β值越大，临界耦合强度λ_c越提前；但是当sin β为负值时，对同步则会产生滞后的影响，并且sin β值越小，临界耦合强度λ_c越滞后。

3 理论分析

为简单起见，可以在星形网络中做理论分析，因为星形网络的结构特点抓住了无标度网络的主要特征^[7]。在星形网络中，依旧设定振子ω_i等于节点的度k_i。

在星形网络中，建立一个旋转坐标系，该系统的平均相位可以表达为ψ(t)=ψ(0)+Ωt，其中Ω为星形网络振子的平均频率^[7]，可以设置ψ(0)=0。如此可以定义ϕ_h=θ_h-ψ=θ_h-Ωt，ϕ_j=θ_j-ψ=θ_j-Ωt，其中θ_h表示中心节点振子的相位，θ_j为边缘节点振子的相位。

设星形网络节点数为(k+1)，中心节点振子与边缘节点振子的运动方程可以写成：

$ \begin{array}{l} \;\;\;\;\;\;\;\dot \phi = \left( {{\omega _h}-\mathit{\Omega }} \right) + \lambda \sum\limits_{i = 1}^k {\left( {\sin \left( {{\phi _j}-{\phi _h}} \right) + \sin \;\beta } \right)}, \\ \;\;\;\;\;\;\;\dot \phi = \left( {\omega-\mathit{\Omega }} \right) + \lambda \left( {\sin \left( {{\phi _j} - {\phi _h}} \right) + \sin \;\beta } \right), j = 1, \\ 2, \cdots, k, \end{array} $

(5)

且

$ {\omega _h} = k, \omega = 1, \mathit{\Omega = }\frac{{2k}}{{k + 1}}, \beta \in \left( {-\frac{\pi }{2}, \frac{\pi }{2}} \right)。$

现在计算锁定相位${\dot \phi _h} = 0$时的情况。利用(2) 式，可以计算得到

$ r{e^{i\psi \left( t \right)}} = \frac{{{e^{i{\theta _h}}} + \sum\limits_{j = 1}^k {\left( {{e^{i{\theta _j}}}} \right)} }}{{1 + k}}, $

两边同乘e^-iθ_h，可得

$ r{e^{-i{{\dot \phi }_h}}} = \frac{{1 + \sum\limits_{j = 1}^k {{e^{i\left( {{\theta _j}-{\theta _h}} \right)}}} }}{{1 + k}}。$

取虚数部分可得

$ -r \cdot \sin \;{\phi _h} = \frac{{\sum\limits_{j = 1}^k {\sin \left( {{\theta _j}-{\theta _h}} \right)} }}{{1{\rm{ + }}k}}, $

并结合式(5) 可得

$ \begin{array}{l} \;\;\;\;\;\;\;\;{{\dot \phi }_h} =-\lambda \cdot \left( {1 + k} \right) \cdot r\sin {\phi _h} + \left( {{\omega _h}-\Omega } \right) + \\ \lambda \cdot k\sin \beta 。\end{array} $

当${{\dot \phi }_h} = 0$时，有

$ \sin \;{\phi _h} = \frac{{\left( {{\omega _h}-\mathit{\Omega }} \right) + \lambda \cdot k\sin \beta }}{{\lambda \cdot \left( {1 + k} \right) \cdot r}}, $

而锁定相位${{\dot \phi }_j} = 0$时，有

$ \sin \left( {{\phi _h}{\rm{-}}{\phi _j}} \right) = \frac{{\mathit{\Omega-}\omega-\lambda \sin \beta }}{\lambda }, $

所以，可以计算得到

$ \begin{array}{l} \;\;\;\;\cos \;{\phi _j} = \frac{{\mathit{\Omega- }\omega \mathit{- }\lambda \sin \;\beta }}{\lambda }\sin {\phi _h} \pm \\ \frac{{\sqrt {\left[{{\lambda ^2}-{{\left( {\mathit{\Omega-}\omega \mathit{-}\lambda \sin \beta } \right)}^2} \cdot \left( {1 - \sin \;{\phi _h}} \right)} \right]} }}{\lambda }。\end{array} $

由此，可得$\lambda \ge \frac{{\mathit{\Omega }-\omega }}{{1 + \sin \;\beta }}$，所以星形网络在锁相情况下的临界耦合强度为${\lambda _c} = \frac{{\mathit{\Omega }-\omega }}{{1 + \sin \;\beta }}$。当sin β＜0时，λ_c会滞后，即同步滞后；当sin β＞0时，λ_c会提前，即同步提前。

4 结论

本研究在无标度网络及星形网络中讨论加入拖拽因子后的Kuramoto模型，通过数值模拟，发现在无标度网络中，会发生爆炸式同步，并且拖拽因子对临界耦合强度产生影响，当拖拽因子为正值时，临界耦合强度减小，同步提前；反之，当拖拽因子为负值时，临界耦合强度增大，同步滞后。理论上分析星形网络中也存在爆炸式同步，拖拽因子对同步起到了与无标度网络中相同的作用。