不同频率涨落驱动下全局耦合谐振子的集体动力学行为*

姜磊 赖莉 蔚涛† 罗懋康2)

1) (四川大学数学学院, 成都 610064)

2) (四川大学空天科学与工程学院, 成都 610064)

对多粒子耦合系统而言, 环境涨落对各粒子的作用在实际情况中往往是互异的, 为此, 本文研究不同频率涨落驱动下全局耦合过阻尼谐振子系统中的集体动力学行为, 包括稳定、同步和随机共振.通过随机平均法推导得出粒子行为的统计同步性, 进而得到了系统平均场与单粒子行为在统计意义下的等价性.并且, 利用该同步性进一步求解得到了输出幅值增益和系统稳定的充要条件.前者为分析系统随机共振行为奠定了理论基础, 后者给出了本文所得结论的适应范围.仿真表明, 耦合强度 ε 的增加或系统规模N的增大会带来两方面的影响: 首先, 稳定区域逐渐增大, 同步时间逐渐缩短; 其次, 系统的有序性增强, 需要更大的噪声强度提供更强的随机性来与之实现最优匹配, 从而关于噪声强度 σ 的随机共振峰逐渐右移, 反之亦然.

1 引 言

噪声在自然和工程系统中无处不在, 而受噪声影响的复杂系统[1]的集体动力学行为一直是非线性科学的研究重点.稳定性作为系统的一个基本结构特性, 一直是随机系统研究所关心的基本问题[2,3].例如, 在天气预报中, 预报的质量就取决于系统模型的稳定性, 即在一段时间后的天气状况对初始观测数据扰动的依赖[4].在飞行器控制系统的设计中, 系统的稳定性则对飞行器的安全性与可操控程度有至关重要的影响[5].

同步作为复杂系统集体行为的一种最基本的表现形式广泛存在于如钟摆、乐器、生物系统、神经等各种自然科学与工程系统中[6,7].作为一种集体行为, 同步会对系统产生重要的影响[8].例如, 鱼群鸟群的同步行为可减少个体被捕食的概率, 电网中电机的同步可以大大提高输电效率; 而另一方面, 过桥时人群的同步踩踏会使桥面产生共振断裂, 大脑神经元突发性异常同步发电会导致癫痫发作.因此, 对复杂系统的同步进行理论分析以掌握其机制往往具有重要意义.

随机共振由于表现出噪声与非线性系统中微弱周期信号的协作现象而受到广泛深入的研究[9−11].这一现象表明, 一定强度的噪声可以放大系统对微弱周期信号的响应幅值[12,13].广义上讲,随机共振指的是某些系统响应函数(如矩、信噪比、幅值增益等)随系统的某些特征参数(如周期信号的频率、噪声强度、相关率等)非单调变化的现象[14].随机共振产生的条件为: 系统的非线性性,微弱周期信号与噪声.因此, 早期关于随机共振的研究主要集中在受加性白噪声驱动的非线性系统中[15−23].由于大部分非线性模型难以解析求解, 对其共振机理的分析往往难以深入, 从而使得随机共振在物理、生物、工程等领域的应用具有较大局限.近期研究成果表明, 在乘性噪声驱动下的线性系统也会出现随机共振现象[14,24−34].并且, 这类线性系统通常可以理论求解, 这便为共振机理研究的深入与随机共振应用的推广提供了有利条件.

频率涨落通常由系统外噪声引起.关于频率涨落谐振子的动力学行为的研究, 从Chandrasekhar[35]的研究开始, 一直受到广泛的关注.Berdichevsky和Gitterman[34]最先给出如下受双态噪声驱动的过阻尼系统的解析解并分析了该系统的随机共振现象:

由于方程(1)在波动障碍穿越[36]、酶动力学[37]及核磁共振[38]等化学、生物、物理多个领域存在应用, 所以关于这类模型的研究具有较高的理论和应用价值.文献[25, 39]将方程(1)中的双态噪声改为非对称双态噪声, 对新的模型进行解析求解, 分析了模型的随机共振现象.并且以一阶线性电路系统为例给出了该模型在新的物理场景中的应用.

以上关于随机共振的研究主要集中在非耦合系统中.最近, 随着复杂网络的兴起, 耦合系统中随机共振的研究正受到越来越多的关注.研究表明, 耦合系统中不仅会产生随机共振现象而且耦合关系还会丰富系统的共振行为[19−23,28,40,41].Pikovsky等[22]发现不同类型的耦合随机系统中都存在系统规模共振.Li[28]研究发现, 在合适的参数下, 调整局部耦合系统中的耦合强度与粒子个数, 既能加强也可以削弱随机共振的强度.Yu等[41]在对双耦合谐振子共振行为的研究中发现, 调整耦合强度不仅可以增强随机共振, 还能改变系统的共振形式.特别地, Yang等[3]将模型(1)推广为N粒子全局耦合模型, 通过理论推导及仿真分析了N个粒子的集体动力学行为.但在该推广模型中, 系统中各粒子的频率涨落被建模为全同噪声.而在实际环境中, 环境涨落在同一时刻对不同粒子的影响应该是有所差异的.因此, 将各粒子的频率涨落建模为独立同分布的乘性噪声更能反映随机耦合系统的真实动力学行为.

鉴于上述分析, 本文研究受不同频率涨落驱动的全局耦合系统的集体动力学行为.其中, 频率涨落建模为独立同分布的乘性对称双态噪声.这样建模主要是基于双态噪声在理论研究和工程应用中的普适性考虑.一方面, 双态噪声有界且形式简单易于实现, 常被应用于各种物理和生物系统中[42−44].另一方面, 在极限条件下, 双态噪声可转化为高斯白噪声和白散粒噪声.这使得双态噪声具有一定的理论研究价值[45].更重要的是, 受乘性双态噪声影响的线性模型通常可以精确求解, 便于进行分析研究.在此基础上, 本文通过随机平均法推导出系统各粒子的同步性.并根据同步性求出系统输出幅值增益G的表达式及稳定的充要条件.最后通过数值仿真分析了系统稳定、同步及随机共振等集体行为.

2 模 型

本文研究受乘性双态噪声和余弦信号共同作用的N粒子全局耦合过阻尼系统, 该系统所满足的微分方程如下:

其中xi(t) 为第i个粒子在t时刻的位移,N为系统规模,为全局耦合项,ε≥0为耦合强度.系统的全局耦合性体现在系统中的每个粒子都与其他粒子直接相连.A0cos(Ωt) 为外部周期驱动力,ω≥0 表示系统的固有频率.该频率的涨落被建模为对称双态噪声ξi(t).该噪声为一个双态Markov过程, 在±σ中随机取值, 其均值和相关函数分别为

其中σ表示噪声幅值,λ表示噪声相关率,τc=1/λ为平均等待时间.

由于噪声ξi(t) 随机波动, 各粒子所处的势场随机地在

之间切换.这类随机势场在很多类型的系统中均有发现.例如生物体内的ATP-ADP循环势场, ATP将能量储存在化学键中, 通过水解反应, 化学键中的一些电子进入到更低的能量态.这一过程产生了ADP及生物所需的能量.反过来, ADP又可以通过消耗一定的能量生成ATP.这一循环过程产生了生物系统中的ATP-ADP涨落势场[46−48].

另一方面, 在液体环境中, 粒子往往表现为过阻尼运动, 即其惯性项可被忽略[49].因此, 与文献[3]类似, 模型(2)可被看作一个受到周期驱动的耦合蛋白质分子马达在液体环境下ATP-ADP势场中的运动模型.文献[28]在局部耦合情况下对这类模型的动力学行为进行了分析, 并列举了josephson结构阵列[50]、神经网络[51]及耦合蛋白子马达[52,53]作为其潜在应用场景.另外需要指出, 若令系统规模N=1 , 则模型(2)退化为文献[34]中的单粒子运动模型;N=2 时, 模型(2)即变为我们之前文章中研究的动力学模型[54]; 若令ξ1(t)=···=ξN(t), 则本文的模型又可退化为之前Yang等[3]研究的全局耦合模型.因此, 从理论和应用两方面考虑, 模型(2)都具有较高的研究价值.

在本文的推导过程中将用到指数关联噪声所满足的Shapiro-Loginov(S-L)公式[55]的推广版本,为此, 先对S-L公式做简要介绍, 再用该公式推导出本文所需的公式版本.

定理1(Shapiro-Loginov公式)η(t) 是一个噪声.假设如下条件成立:

(1)〈η(t)〉=0,〈η(t)ξ(s)〉=σ2e−λ|t−s|,

(2)ψ(t) 是指数相关噪声η(t) 的泛函,

则有如下公式成立:

为了对模型(2)进行求解, 以下将(5)式推广到多个相互独立指数关联噪声的情形.

定理2(Shapiro-Loginov公式)ηi(t),(i=1···k)为k个相互独立的噪声.假设如下条件成立:

(1)〈ηi(t)〉=0,〈ηi(t)ηi(s)〉=σ2e−λ|t−s|,i=1···k,

(2)ψ(t) 是指数相关噪声ηi(t) 的泛函,

则有如下公式成立:

证明显然,η1(t)η2(t)···ηk(t) 仍然是一个噪声(即随机过程).以下证明该噪声满足定理1的两个条件.由ηi(t),(i=1···k) 的相互独立性及条件(1)可得:

故η1(t)η2(t)···ηk(t) 仍然是一个指数关联的随机噪声.由条件(2)可知,ψ(t) 是噪声η1(t)η2(t)···ηk(t)的泛函.

故由定理1得(6)式成立.

由(6)式直接可得:

其中

为集合{ξi1,ξi2,···,ξik}所包含元素的个数.1≤i1,i2,···,ik≤N,i1,i2,···,ik互不相等.

3 理论结果

3.1 同步性

本节主要研究上述耦合系统中各粒子行为是否具有同步性, 由于各粒子所受驱动噪声相互独立, 因而, 这里所说的同步性只能是统计意义下的同步.为此, 设法构造包含变量〈xi〉,i=1,2,···,N的封闭微分方程组.用 1,ξi1ξi2···ξik分别乘以(2)式两端, 因为形如ξi1ξi2···ξik这种k个互不相同的噪声相乘的形式的因子共有个.当k取遍1到N, 再算上常数项因子 1 , 共有 2N个线性无关的乘因子.因而可以得到 2NN个随机微分方程, 利用S-L公式取平均, 即得到包含 2NN个变量的2NN个线性无关的非齐次常微分方程组.未知量如表1所列.

3.1.1 变量分组

接下来对表1中的变量进行分组, 将相等的变量分到一组, 以约简变量个数.各组变量间的等价关系会在下一步予以证明.

表1 变量表Table 1.Variable table.

1) 由于模型(2)中各粒子位置对称, 表1中第一行中的各变量应该是两两相等的.同理, 最后一行中的各变量也应该是两两相等的.分别用U0,UN表示两行元素构成的集合.即:

2) 若k̸=l, 则于是, 我们记

3) 对于Uk(这里假定 1 ≤k≤N−1 )中的元素〈ξi1ξi2···ξikxi〉 , 若i∈{i1,i2,···,ik}, 则表示组成该乘因子的噪声中既有直接作用于粒子xi的, 也有先作用于其他粒子再通过耦合作用间接作用于粒子xi的.若i∈/{i1,i2,···,ik}, 则表示组成该乘因子的噪声都是间接作用于粒子xi的.为此将集合Uk分解为两个集合Uk1,Uk2的不交并:

那么Uk1中的变量和Uk2中的变量应该是不相等(同步)的.

3.1.2 组内变量等价性证明

通过上述步骤, 将表1中的变量划分到了2N个集合中, 并猜测各集合中的变量相等.下面来证明这一猜测.

1) 对于集合Uk1, 从中任意选定一个元素v1,将Uk1中除v1以外的元素分别与v1作差, 得到一个新的集合:且v1∈Uk1,v∈Uk1}.显然Uk1中任意两个元素的差属于Vk1或者可以表示为Vk1中两个元素的差, 且有|Vk1|+1=|Uk1|.

2) 对剩余的 2N−1 个集合进行与1)完全相同的操作.

得到 2N个集合V0,···,Vk1,Vk2,···,VN,这些集合的总元素个数为 2NN−2N个.显然这2NN−2N个元素线性无关.下面构造一个包含这2NN−2N个变量的 2NN−2N个线性无关的齐次常微分方程组.

对于新得到的 2NN−2N个变量中任意一个变量〈ξm1ξm2···ξmkxi〉−〈ξn1ξn2···ξnkxj〉.要得到包含该变量的微分方程, 只需用ξm1ξm2···ξmk,ξn1ξn2···ξnk分别乘以

并两端相减, 可得

1) 若i∈/{m1,m2,···,mk},j∈/{n1,n2,···,nk}, 则有

对(15)式取均值并应用S-L公式(9)式得

显然(16)式中各均值号中的变量均可由Vk1等 2N个集合中的 2NN−2N个元素线性表示.

2) 若i∈{m1,m2,···,mk},j∈{n1,n2,···,nk}, 按照1)中同样的步骤可得:

其中

显然(17)式中均值号中的变量也都可以由Vk1等2N个集合中的 2NN−2N个元素线性表示.

(16)式或(17)式即为包含变量〈ξm1ξm2···ξmkxi〉−〈ξn1ξn2···ξnkxj〉的微分方程.由此, 对待求的 2NN−2N个变量, 可以得到包含这些变量的2NN−2N个线性无关的齐次常微分方程组.由于本文主要关注粒子在长时间范围内的解的状况, 而方程组的初始值对粒子的运动状态的影响会随着t的增大逐渐消失.所以这里不妨假设该方程组初始条件为0.利用Laplace变换法可得上述方程组((16)式和(17)式)只有零解.由此推导出各集合U0,···,Uk1,Uk2,···,UN中变量的相等关系.

3.1.3 同步性结论与意义

通过上述分析, 得到了各集合U0,···,Uk1,Uk2,···,UN中变量的相等关系.也即同步性结论如下:

上述同步性结论具有以下两条重要意义.

1) 从同步性(19)式可以看出

又因为系统的平均场

从而, 由(20)式可知

对复杂耦合系统集体动力学行为的研究往往围绕其平均场行为展开, (22)式表明: 在本模型中,系统的平均场行为与单粒子行为具有统计一致性,从而, 本文可通过研究单粒子的动力学行为来开展对系统集体动力学行为的研究.

2) 为研究单粒子的动力学行为, 需对其一阶矩〈xi〉 进行求解.同步性(19)式的另一个重要意义在于: 可极大地简化模型的求解过程, 使对〈xi〉 的解析求解成为可能.在下一节的的分析中, 将利用该同步性((19)式)完成对〈xi〉 和相应输出幅值增益G的具体解析求解.

3.2 输出幅值增益

本节利用3.1节所得结果推导一阶矩〈xi〉 及系统的输出幅值增益G的表达式.首先将可能出现的变量用新的符号表示如下:

其中 1 ≤i≤N, 1 ≤k≤N−1.

对(2)式中第i个式子两边取平均, 并利用同步性得到含有变量y0的方程如下:

(24)式出现了新的耦合项y1,1, 为了建立包含该变量的微分方程, 用ξi1ξi2···ξik乘以第i个式子两边.这里i∈{i1,i2,···,ik}, 则有

对(25)式两边取平均并利用S-L公式(9)式得

对(26)进行整理得

由于(27)式中又出现了新的变量yk,2.同理, 当i∈/{i1,i2,···,ik}时, 按照同样的步骤, 可得如下包 含变量yk,2的微分方程式:

最后, 用ξ1ξ2···ξN乘以(2)式中第i个式子两边,计算可得包含变量yN的微分方程式:

所得的 2N个式子((24)式, (27)式, (28)式,(29)式)可改写成如下矩阵形式:

其中,

其中T表示矩阵的转置.

A为三对角矩阵, 即除了三条对角线上的元素外其余元素都为零.

对(30)式作Laplace变换, 有

不妨记通过求解线性方程组(35),可以求得向量.特别地, 有

其中

其中h0(t) 为H0(s) 的Laplace逆变换.

另一方面, 系统((24)式, (27)式, (28)式,(29)式)可看作一个线性时不变系统.根据线性时不变系统的响应理论, 该系统的输出信号是一个与输入信号仅在幅值和相位上有差异的余弦信号.即

其中A和ϕ分别为〈xi〉 的振幅和相位.将(40)式作Laplace变换并对比(37)式可得

最后求得输出振幅增益为

根据(42)式, 通过数值仿真分析不同的系统参数变化对G的影响, 可以判断系统有无随机共振现象产生.在此之前, 首先对系统(2)的稳定性做一些简要分析.

3.3 系统稳定性

系统(2)稳定等价于特征方程

的根都分布在左半复平面, 即矩阵A的特征值都分布在右半复平面.在此基础上, 本节建立了系统(2)稳定的一个充分条件和一个充要条件.

定理3(充分条件) 当σ2<ω(ω+λ) 时, 系统(2)是稳定的.

证明1) 首先当σ=0 时, 由同步性(20)式可知, 系统(2)为含一个变量的常微分方程, 又由于ω>0, 系统显然是稳定的.

2) 当σ>0 时, 证明矩阵A的所有特征值都为正实数.以此推出系统(2)的稳定性.

(i) 先对ε>0 的情况予以证明.此时, 有bi+1ci>0,1≤i≤N−1.故A相似于如下实对称三对角矩阵[56]:

A与S有相同的特征值且都为实数.下面证明当定理条件满足时,S为正定矩阵.

令S=S1+S2.其 中S1为ε=0 时S退 化 得到的矩阵.S2=S−S1.显然S1,S2都是分块对角矩阵.其中S1的每个分块是如下形式的二阶方阵:

其中 0 ≤k≤N−1.S2的每个分块为零矩阵或如下形式的二阶方阵:

其 中 1 ≤k≤N−1.显然,S1,S2都为实对称矩阵.易证S2半正定.以下证明定理条件满足时,S1正定.因为分块矩阵各个分块的特征值的并集即为该分块矩阵的全部特征值.故要让S1正定只需S1每个分块的特征值为正即可.易得, 这只需要

(ii) 当ε=0 时, 矩阵A退化为分块对角矩阵.用与(i)中同样的步骤证明即可.

下面, 为给出系统(2)稳定的充要条件, 将特征方程(43)式看成是关于σ2和s的函数来展开分析, 从而将其改记为f(σ2,s).特别地, 若将σ2看成参数, 则f(σ2,s) 是关于变量s的 2N次多项式.

定理4(充要条件) 系统(2)稳定当且仅当

其中为f(σ2,0)=0 的最小正根.

证明因为矩阵A相似于实对称矩阵(定理3), 故对任意固定的σ2,f(σ2,s)=0 有 2N个实根.将这 2N个实根从大到小排列为

则si(σ2)(1≤i≤2N) 为σ2的连续函数[57].显然, 系统(2)稳定等价于s1(σ2)<0.故由定理3 知, 当0<σ2<ω(ω+λ) 时 ,s1(σ2)<0.另一方面, 随着σ的增加, 系统势函数的不稳定性逐渐增强; 从而, 当σ充分大时, 系统必然是不稳定的, 此时s1(σ2)>0.所以根据s1(σ2) 的连续性可知, 存在σ>0 , 使得s1(σ2)=0.记

由于本文只关心系统稳定时的各种动力学行为, 所以后续对同步行为与随机共振行为的分析都是在稳定性条件(48)式成立的前提下进行的.

4 数值仿真算法

本节利用随机Taylor展开算法对模型(2)进行数值模拟, 以此验证算法的有效性.对于充分小的时间步长 ∆t, 模型(2)在离散时间下的近似表达式为

其中

其中, ∆t=tn−tn−1为时间步长, ∆Xi为单个时间步长内双态噪声的增量, 可通过文献[58]中的方法仿真求得.

为了验证仿真算法的有效性, 需要将仿真结果与理论结果进行对比.首先, 进行K次Monte Carlo实验并记第i个粒子第k次Monte Carlo实验的轨迹为则可得到第i个粒子位移一阶矩〈xi(t)〉的数值近似为

之后, 根据(40)式可知, 输出信号〈xi(t)〉 为与输入信号A0cos(Ωt) 有相同频率的余弦信号.所以对〈xi(t)〉num作Fourier变换即可得到〈xi(t)〉 的仿真幅值Anum.最后利用(42)式即得输出幅值增益的数值仿真结果Gnum.图1给出了仿真算法的一个具体的例子.其中, 图1的第一排子图为10个粒子进行 1 04次Monte Carlo仿真所得平均位移, 各粒子的初始位移服从标准正态分布.对该平均位移作Fourier变换得到第二排子图对应的频域值.由此得到在该组参数下输出幅值增益的数值仿真结果Gnum=1.11708 与理论值G=1.11677 的相对误差仅为 2.78×10−4.

图1 粒子平均位移及对应频域值的数值仿真结果, 其中ε=1 , A 0=1 , σ =1.3,N=10 , ω =0.5 , λ =1 , Ω=π/4 , ∆ t=10−3 , K=104,T=120Fig.1.The average realization and the corresponding frequency domain representation with ε =1 , A 0=1 , σ=1.3,N=10, ω =0.5 , λ =1 , Ω =π/4 , ∆ t=10−3 , K=104,T=120.

以下分别讨论时间步长 ∆t、Monte Carlo仿真次数K和总仿真时长T对数值算法收敛性的影响.图2(a)—图2(c)分别描绘了Gnum作为上述三个仿真参数的函数的曲线图.可以看出, 随着仿真∆t的减小及K和T的增加, 仿真结果逐渐趋近于理论结果.且当 ∆t≤0.001,K≥10000 ,T≥120时, 仿真结果与理论结果的误差小于 6×10−4, 对应的相对误差仅为 5.372×10−4.因此, 若不特别指出, 本文数值仿真结果都为固定参数∆t=0.001,K=10000,T=120 下仿真求得.

图2 仿真结果与理论结果对比图, 其中Ω=π/4,ω=0.5,λ=1,ε=1,A0=1,σ=2 (a) K =10000,T=120 ;(b) ∆ t=0.001,T=120 ; ( c) K=10000,∆t=0.001Fig.2.Comparison of theoretical and simulation results with Ω =π/4,ω=0.5,λ=1,ε=1,A0=1,σ=2 :(a) K =10000,T=120 ; ( b) ∆ t=0.001,T=120 ;(c) K =10000,∆t=0.001.

5 集体动力学行为分析

5.1 系统稳定性

在3.3节中, 推导得出了粒子位移一阶矩稳定的充要条件.这里进一步通过仿真分析稳定性区域随参数的变化及其对粒子运动行为的影响.图3(a)给出了不同的N所对应的系统稳定性区域图像.其中黑色区域为N=2 所对应的稳定性区域; 黑色区域与深灰色区域的并为N=5 所对应的系统稳定性区域; 黑色区域, 深灰色区域与浅灰色区域的并集则为N=10 时的稳定性区域.从中可以看出,系统的稳定性区域随耦合强度ε及系统规模N的增加而增大.为了更直观地反映系统稳定性区域随N的变化, 图3(b)给出了在一组固定的参数下,系统稳定性区域边界随N的变化曲线.可以看出,稳定性区域随着N的增大单调递增.关于图3中现象的一个合理的物理解释为: 粒子规模N及耦合强度ε决定粒子间的耦合力Fi, 而耦合力为系统输出提供有序性.噪声为粒子位移引入无序性, 且噪声幅值σ越大, 粒子位移的无序性越强.有序性和无序性相互竞争决定系统输出的稳定性.当σ∈[0,σmin)时, 耦合力提供的有序性占主导地位.当σ=σmin时, 有序与无序的竞争达到一个临界点,使得初值的影响不会随时间增加而消失.系统输出开始变得不稳定.当σ>σmin时, 噪声引入的无序性占主导地位, 系统输出无界.由于增加N及ε可以增加耦合力Fi, 使得系统输出的有序性增强.此时就需要更大的σ才能使得系统输出达到稳定性的临界点.所以, 系统的稳定性区域会随着粒子规模N及耦合强度ε的增加而增大.

图3 ( a) 不同的N所对应的系统稳定性区域, 其中ω=0.5,λ=1 ; ( b) 系统稳定性区域边界随N的变化曲线, 其中 ε =1 , 其他参 数同图(a)Fig.3.( a) System stability region corresponding to different N with ω =0.5,λ=1 ; ( b) curve of the boundary of the system stability region with N with ε =1 , other parameters are the same as panel ( a).

为了直观反应稳定性区域对粒子平均位移的影响.不失一般性, 这里固定N=10 并从图3(a)中选择A(1.5,1),B(2,1) 两点进行仿真实现.结果如图4所示.从图4(a)可以看出, 粒子平均位移有界, 此时系统是稳定的.图4(b)中的系统输出随着时间的增大向无穷远发散, 此时系统是不稳定的.这与定理4的结论相符.

图4 粒子平均位移实现, 其中ω=0.5,λ=1,Ω=π/4,ε=1,A0=1,N=10 (a)σ=1.5(b)σ=2Fig.4.The average displacements of particles with ω=0.5,λ=1,Ω=π/4,ε=1,A0=1,N=10:(a)σ=1.5 ; ( b) σ =2.

由于本文只关心系统稳定时的各种动力学行为, 所以以下两节的仿真都是在稳定性条件成立的前提下进行的.

5.2 粒子的集体同步行为

本节通过数值仿真对3.1节中的结论(20)式予以验证.即随着系统时间的演化,N个粒子的平均位移会达到同步状态.

图5给出了系统(2)在不同参数条件下不同粒子平均位移的仿真结果和对应的方差, 其中各粒子的初始位置服从标准正态分布.从图5(a)—图5(c)的第一行子图中可以看出, 具有不同初始值的粒子随着时间的增加最终都会达到同步状态.这与3.1节中的结论相符, 即方程组(16)式, (17)式的初值条件不改变长时间后系统中粒子的同步状态.为直观反映粒子的同步速度, 引入粒子平均位移〈xi(t)〉num的方差∆(t) 如下:

其中

图5(a)—图5(c)的第二行子图给出了方差随时间变化的曲线图.从中可直观地看出, 图5(a)中曲线的下降速度慢于图5(b)和图5(c).且当t=1 时,三幅图对应的方差∆分别是4.21×10−3,9.2×10−7,1.2×10−7.下面通过定义同步时间t0, 对不同参数下系统(2)的同步速度进行定量描述.由于系统(2)中各噪声ξi独立同分布, 且实验仿真次数有限, 所以∆随着时间增加不会降为零.因此, 这里假设∆(t)≤10−5时, 即认为系统达到同步状态, 并记相应的同步时间为:

将图5(a)—图5(c)第二行子图中各方差曲线图横纵坐标以对数形式表示得到第三行子图.三幅子图中系统达到同步的时间分别是 3.53,0.72,0.55.对比三幅图中的同步时间可知, 增加耦合强度ε和系统规模N可显著加快粒子的同步速度.这是由于一方面增加耦合强度ε使得粒子间的耦合力

图5 系统(2)在不同参数条件下的仿真实现和所对应的方差, 其中ω=0.5,λ=1,σ=0.7,A0=1,Ω=π/4 (a)ε=1,N=2 ; ( b) ε =4,N=2 ; ( c) ε =1,N=10.三幅子图的顶部图中不同颜色的实线代表不同粒子的平均位移Fig.5.The realization of the system (2) under different parameter conditions and the corresponding variance with ω=0.5,λ=1,σ=0.7,A0=1,Ω=π/4 : ( a) ε =1,N=2 ; ( b) ε =4,N=2 ; ( c) ε =1,N=10.The solid lines in different colors in the first panel of the three subfigure represent the average displacement of different particles.

加大, 而耦合力的作用是使得各粒子相互吸引, 这使得系统可以更快同步.另一方面由于系统中的耦合关系为全局耦合, 增加N也会加大粒子间的吸引力, 从而让各粒子更快达到同步状态.另外, 从图6也可以定量地看出, 随着N的增加, 系统的同步时间t0逐渐变小.

图6 同步时间 t0 随群体规模N的变化曲线, 所选系统参数与图5相同Fig.6.The change curve of synchronization time t0 with system size N.The selected system parameters are the same as the Fig.5.

最后需要指出, 耦合强度只改变系统的同步时间而不改变系统达到稳定时的同步状态.这是由于系统中各粒子的位移主要受周期外力、粒子所处势场、各粒子的初始位置以及粒子间的耦合力四个因素影响.这四个因素中, 决定粒子同步的本质因素为周期外力和粒子所处势场.一方面, 各粒子受相同的周期外力作用; 另一方面, 各粒子有相同的本征频率, 频率噪声有相同的统计性质, 这使得各粒子所处的势场从统计意义看是相同的.所以一段时间后, 系统初值对各粒子位移的影响消失.耦合强度的存在只是加速了粒子的同步速度, 当ε=0 时各粒子所满足的运动学方程相同, 此时各粒子仍然会同步.故耦合强度不改变系统的同步性.

5.3 随机共振

随机共振现象是噪声的随机性与驱动力的有序性在非线性作用下相互协作的结果.在本文的模型中, 随机性由噪声提供, 有序性由正弦驱动力和耦合作用力两部分组成.由于正弦驱动力的作用在以往的文献中已展开过大量充分的讨论, 因而本文主要分析噪声与耦合力的协作现象.它们相互协作产生随机共振的机理如下: 噪声ξi(i=1,2,···,N)为系统提供随机性, 噪声幅值σ越大, 则随机性越强; 耦合力Fi(i=1,2,···,N) 为系统提供有序性,耦合强度ε或系统规模N越大, 则有序性越强.在“乘性”这种非线性作用方式下, 有序与无序相互竞争, 当两者间达到某种最优匹配时, 系统产生随机共振, 从而使得系统输出达到最大.

在上述共振机理的指导下, 接下来通过数值分析的方式, 逐一对输出幅值增益G关于耦合强度ε、系统规模N以及噪声幅值σ的共振行为进行具体分析.

5.3.1 关于耦合强度ε的参数随机共振

图7(a)和图7(b)分别给出了不同的σ,N下,G随ε的变化曲线.可以看出, 除了σ=0 外, 其他曲线都有不同强度的共振峰.且共振峰主要出现在低耦合区域.随着耦合强度ε的增大, 每条曲线都趋于一个定值.这是因为当参数σ,N固定时, 系统输出有序性主要由ε确定.且系统输出有序性随ε增加而增强.当ε过小时, 过弱的耦合力不能有效协调粒子的运动.当ε过大时, 过强的耦合力会使粒子相互吸引组成一个整体而完全抑制噪声带来的随机性的影响.两种情况都限制了噪声与输入信号的协作效应.因此, 使得输出幅值增益G最大的ε一定是介于两者之间的.

图7 不同的 σ , N下 G (ε) 的变化曲线, 其中ω=0.5,λ=1,A0=1,Ω=π/4 (a)N=10 ;(b)σ=0.7Fig.7.The curve of G (ε) under different σ and N with ω=0.5,λ=1,A0=1,Ω=π/4 : ( a) N =10 ; ( b) σ =0.7.

图7 (a)中, 当σ=0 时,G取常值.这是由于当σ=0时, 噪声消失.此时模型(2)退化为确定性模型, 统计意义下的同步〈x1〉=〈x2〉=···=〈xN〉 变为x1=x2=···=xN.耦 合 力Fi(i=1,2,···,N)使各粒子相互吸引组成一个整体, 之后变为零.因此,ε对G不产生影响.另一方面, 将σ=0 代入(42)式直接计算可得:

当σ≥0 时, 共振峰的高度随着σ的增大而增大.这说明增加噪声幅值σ可加大随机共振的强度.

图7(b)中, 随着系统规模N的增加, 共振峰变得越来越尖锐, 同时共振峰的位置左移, 共振峰的高度不断增加.产生这种现象的原因如下: 由于系统中的耦合关系为全局耦合, 增加N会加大粒子间的吸引力, 所以当N较大时, 随着ε的增加,粒子间的吸引力增大得更快.因此G随ε的变化曲线也更陡.同样地, 更强的耦合力进一步抑制了噪声带来的随机性, 使得G不能达到最优.所以当噪声幅值σ固定时, 增加N会使得共振峰左移.另外注意到, 图7(b)中所有曲线都有共同的起点.这是由于当粒子间的耦合作用消失时(ε=0 ), 模型(2)变为N个彼此独立且相同的系统.此时, 改变系统规模N并不改变系统中各粒子的受力状况.因此N的变化不影响G.同样地, 将ε=0 代入(42)式可得:

所以, 从(58)式也可得出上述结论.

5.3.2 关于系统规模N的参数随机共振

图8(a)和图8(b)分别给出了不同的σ,ε下,G随N的变化曲线.可以看出, 除了σ=0,ε=0外, 其他曲线都有不同强度的单峰共振.这表明在合适的参数下, 存在最优的系统规模N使得输出幅值增益G达到最优.这里系统中有序与无序竞争关系的变化主要由系统规模N的变化引起.在给定参数ε,σ的条件下, 系统规模N的增加使粒子间的耦合力Fi(i=1,2,···,N) 变大, 进而使得系统输出有序逐渐增大.当系统规模N较小时, 过弱的耦合力不能有效协调粒子的运动, 即噪声引入的随机性过大.当系统规模N较大时, 过强的耦合力会完全抑制噪声带来的随机性的影响, 使系统(2)近似于确定性系统.这两种极端情况都会限制粒子的平均运动.因此, 使得输出幅值增益G最大的N一定是介于两者之间的.

图8(a)中, 当ε=0 时,G取常值.当ε>0 时,随着耦合强度ε的增加, 使共振曲线趋于稳定值的N减小.这与图7(b)中的结论一致.另外还可以看出, 随着ε的增加, 共振曲线下移.这表明增加耦合强度ε会减弱系统关于N的参数随机共振强度.

图8(b)中, 当σ=0 时,G取常值.这是由于此时模型(2)退化为确定性模型, 而耦合力Fi(i=1,2,···,N)使各粒子相互吸引组成一个整体, 之后变为零.此时, 不论粒子个数多少, 系统中的各粒子都只受周期外力A0cos(Ωt) 的作用.所以改变系统规模N不影响输出幅值增益.该结论也与(57)式一致.当σ>0 时, 随着σ的增加, 共振峰右移, 共振强度加大.这是由于更大的σ会使粒子的运动轨迹产生更强的随机性, 过强的随机性又会抑制系统的输出.这种情况下, 要使系统达到最优输出就需要加强粒子间的吸引力使粒子运动轨迹更加有序.所以当σ增加时, 使系统输出达到最大的N增大.另外还可以看出, 增加噪声幅值σ可以加强系统的参数随机共振强度.

图8 不 同 的 σ , ε 下 G (N) 的 变 化 曲 线, 其 中ω=0.5,λ=1,A0=1,Ω=π/4 (a)σ=0.7 ;(b)ε=1Fig.8.The curve of G (N) under different σ and ε with ω=0.5,λ=1,A0=1,Ω=π/4: ( a) σ =0.7 ; ( b) ε =1.

5.3.3 关于噪声幅值σ的随机共振

图9(a)和图9(b)分别给出了不同的ε,N下,G随σ的变化曲线.可以看出, 所有曲线都表现出不同强度的共振现象.这表明, 在合适的参数下,系统关于噪声幅值σ会产生随机共振.从图9(a)和图9(b)还可以看出, 当σ=0 时, 所有曲线都有共同的起点.这一现象与图7(a), 图8(b)及(57)式得到的结论一致.

图9(a)中, 随着耦合强度ε从0增加到1, 共振强度变大且共振峰位置右移.当ε从1继续增加到10时, 共振强度逐渐减弱, 并且也逐渐趋向于一个极限值, 而共振峰的位置相对变化较小.这说明增加耦合强度既能增强也可以削弱随机共振现象.

图9(b)中, 随着系统规模N的增加, 共振强度变大且共振峰位置右移.出现这一现象的原因与图8(b)类似.更大的系统规模N使系统各粒子间的相互吸引力变强, 而过强的吸引力进一步抑制了噪声带来的随机性, 使系统(2)产生与确定性系统类似的有序输出.因此, 随着N的增加, 使系统输出G达到最大的σ变大, 共振峰右移.图9(b)表明, 增加系统规模N可增强系统的随机共振强度.

图9 不同的 ε , N下 G (σ) 的变化曲线, 其中ω=0.5,λ=1,A0=1,Ω=π/4 (a)N=10;(b)ε=1Fig.9.The curve of G (σ) under different ε and N with ω=0.5,λ=1,A0=1,Ω=π/4: ( a) N =10 ; ( b) ε =1.

6 结 论

本文研究了不同频率噪声驱动下过阻尼全局耦合系统的集体动力学行为.理论方面, 主要通过随机平均法来展开研究.首先, 得到了各粒子行为的统计同步性, 该同步性使得系统平均场行为与单粒子行为具有统计一致性, 从而可通过对单粒子行为的研究来完成对平均场行为(也即集体行为)的研究.其次, 利用该同步性, 完成了对系统输出幅值增益的解析求解, 为进一步分析系统的随机共振行为奠定了理论基础.最后, 推导给出了系统稳定的充分条件和充要条件, 为本文所得结论的适用范围给出了参考.数值仿真方面, 主要通过随机Taylor展开算法进行研究.首先, 分析了粒子规模N及耦合强度ε对稳定性区域与同步时间的影响, 结果表明, 随着耦合强度ε的增加或系统规模N的增大,粒子间的耦合力增大, 系统有序性增强, 从而使得稳定区域逐渐增大, 同步时间逐渐缩短.其次, 对系统的随机共振行为展开了研究.噪声为系统提供随机性, 耦合力为系统提供有序性, 两者相互竞争,从而使得系统输出关于噪声强度σ、耦合强度ε与系统规模N均表现出随机共振行为.随着耦合强度的增加或系统规模的增大, 系统的有序性增强,需要更大的噪声强度提供更强的随机性来与之实现最优匹配, 从而关于噪声强度σ的共振峰逐渐右移.反之, 随着噪声强度σ的增大, 关于耦合强度ε与系统规模N的共振峰也会右移.本文的研究工作虽然围绕特定耦合模型展开, 但相关结论具有很强的普适性, 对其他复杂耦合系统的研究具有很好的理论指导和参考意义.后续, 将在本文工作的基础上, 进一步结合实际需求展开深入研究.