不同頻率漲落驅(qū)動下全局耦合諧振子的集體動力學(xué)行為*

姜磊 賴?yán)?蔚濤? 羅懋康2)

1) (四川大學(xué)數(shù)學(xué)學(xué)院, 成都 610064)

2) (四川大學(xué)空天科學(xué)與工程學(xué)院, 成都 610064)

對多粒子耦合系統(tǒng)而言, 環(huán)境漲落對各粒子的作用在實際情況中往往是互異的, 為此, 本文研究不同頻率漲落驅(qū)動下全局耦合過阻尼諧振子系統(tǒng)中的集體動力學(xué)行為, 包括穩(wěn)定、同步和隨機共振.通過隨機平均法推導(dǎo)得出粒子行為的統(tǒng)計同步性, 進而得到了系統(tǒng)平均場與單粒子行為在統(tǒng)計意義下的等價性.并且, 利用該同步性進一步求解得到了輸出幅值增益和系統(tǒng)穩(wěn)定的充要條件.前者為分析系統(tǒng)隨機共振行為奠定了理論基礎(chǔ), 后者給出了本文所得結(jié)論的適應(yīng)范圍.仿真表明, 耦合強度 ε 的增加或系統(tǒng)規(guī)模N的增大會帶來兩方面的影響: 首先, 穩(wěn)定區(qū)域逐漸增大, 同步時間逐漸縮短; 其次, 系統(tǒng)的有序性增強, 需要更大的噪聲強度提供更強的隨機性來與之實現(xiàn)最優(yōu)匹配, 從而關(guān)于噪聲強度 σ 的隨機共振峰逐漸右移, 反之亦然.

1 引 言

噪聲在自然和工程系統(tǒng)中無處不在, 而受噪聲影響的復(fù)雜系統(tǒng)[1]的集體動力學(xué)行為一直是非線性科學(xué)的研究重點.穩(wěn)定性作為系統(tǒng)的一個基本結(jié)構(gòu)特性, 一直是隨機系統(tǒng)研究所關(guān)心的基本問題[2,3].例如, 在天氣預(yù)報中, 預(yù)報的質(zhì)量就取決于系統(tǒng)模型的穩(wěn)定性, 即在一段時間后的天氣狀況對初始觀測數(shù)據(jù)擾動的依賴[4].在飛行器控制系統(tǒng)的設(shè)計中, 系統(tǒng)的穩(wěn)定性則對飛行器的安全性與可操控程度有至關(guān)重要的影響[5].

同步作為復(fù)雜系統(tǒng)集體行為的一種最基本的表現(xiàn)形式廣泛存在于如鐘擺、樂器、生物系統(tǒng)、神經(jīng)等各種自然科學(xué)與工程系統(tǒng)中[6,7].作為一種集體行為, 同步會對系統(tǒng)產(chǎn)生重要的影響[8].例如, 魚群鳥群的同步行為可減少個體被捕食的概率, 電網(wǎng)中電機的同步可以大大提高輸電效率; 而另一方面, 過橋時人群的同步踩踏會使橋面產(chǎn)生共振斷裂, 大腦神經(jīng)元突發(fā)性異常同步發(fā)電會導(dǎo)致癲癇發(fā)作.因此, 對復(fù)雜系統(tǒng)的同步進行理論分析以掌握其機制往往具有重要意義.

隨機共振由于表現(xiàn)出噪聲與非線性系統(tǒng)中微弱周期信號的協(xié)作現(xiàn)象而受到廣泛深入的研究[9?11].這一現(xiàn)象表明, 一定強度的噪聲可以放大系統(tǒng)對微弱周期信號的響應(yīng)幅值[12,13].廣義上講,隨機共振指的是某些系統(tǒng)響應(yīng)函數(shù)(如矩、信噪比、幅值增益等)隨系統(tǒng)的某些特征參數(shù)(如周期信號的頻率、噪聲強度、相關(guān)率等)非單調(diào)變化的現(xiàn)象[14].隨機共振產(chǎn)生的條件為: 系統(tǒng)的非線性性,微弱周期信號與噪聲.因此, 早期關(guān)于隨機共振的研究主要集中在受加性白噪聲驅(qū)動的非線性系統(tǒng)中[15?23].由于大部分非線性模型難以解析求解, 對其共振機理的分析往往難以深入, 從而使得隨機共振在物理、生物、工程等領(lǐng)域的應(yīng)用具有較大局限.近期研究成果表明, 在乘性噪聲驅(qū)動下的線性系統(tǒng)也會出現(xiàn)隨機共振現(xiàn)象[14,24?34].并且, 這類線性系統(tǒng)通?？梢岳碚撉蠼? 這便為共振機理研究的深入與隨機共振應(yīng)用的推廣提供了有利條件.

頻率漲落通常由系統(tǒng)外噪聲引起.關(guān)于頻率漲落諧振子的動力學(xué)行為的研究, 從Chandrasekhar[35]的研究開始, 一直受到廣泛的關(guān)注.Berdichevsky和Gitterman[34]最先給出如下受雙態(tài)噪聲驅(qū)動的過阻尼系統(tǒng)的解析解并分析了該系統(tǒng)的隨機共振現(xiàn)象:

由于方程(1)在波動障礙穿越[36]、酶動力學(xué)[37]及核磁共振[38]等化學(xué)、生物、物理多個領(lǐng)域存在應(yīng)用, 所以關(guān)于這類模型的研究具有較高的理論和應(yīng)用價值.文獻[25, 39]將方程(1)中的雙態(tài)噪聲改為非對稱雙態(tài)噪聲, 對新的模型進行解析求解, 分析了模型的隨機共振現(xiàn)象.并且以一階線性電路系統(tǒng)為例給出了該模型在新的物理場景中的應(yīng)用.

以上關(guān)于隨機共振的研究主要集中在非耦合系統(tǒng)中.最近, 隨著復(fù)雜網(wǎng)絡(luò)的興起, 耦合系統(tǒng)中隨機共振的研究正受到越來越多的關(guān)注.研究表明, 耦合系統(tǒng)中不僅會產(chǎn)生隨機共振現(xiàn)象而且耦合關(guān)系還會豐富系統(tǒng)的共振行為[19?23,28,40,41].Pikovsky等[22]發(fā)現(xiàn)不同類型的耦合隨機系統(tǒng)中都存在系統(tǒng)規(guī)模共振.Li[28]研究發(fā)現(xiàn), 在合適的參數(shù)下, 調(diào)整局部耦合系統(tǒng)中的耦合強度與粒子個數(shù), 既能加強也可以削弱隨機共振的強度.Yu等[41]在對雙耦合諧振子共振行為的研究中發(fā)現(xiàn), 調(diào)整耦合強度不僅可以增強隨機共振, 還能改變系統(tǒng)的共振形式.特別地, Yang等[3]將模型(1)推廣為N粒子全局耦合模型, 通過理論推導(dǎo)及仿真分析了N個粒子的集體動力學(xué)行為.但在該推廣模型中, 系統(tǒng)中各粒子的頻率漲落被建模為全同噪聲.而在實際環(huán)境中, 環(huán)境漲落在同一時刻對不同粒子的影響應(yīng)該是有所差異的.因此, 將各粒子的頻率漲落建模為獨立同分布的乘性噪聲更能反映隨機耦合系統(tǒng)的真實動力學(xué)行為.

鑒于上述分析, 本文研究受不同頻率漲落驅(qū)動的全局耦合系統(tǒng)的集體動力學(xué)行為.其中, 頻率漲落建模為獨立同分布的乘性對稱雙態(tài)噪聲.這樣建模主要是基于雙態(tài)噪聲在理論研究和工程應(yīng)用中的普適性考慮.一方面, 雙態(tài)噪聲有界且形式簡單易于實現(xiàn), 常被應(yīng)用于各種物理和生物系統(tǒng)中[42?44].另一方面, 在極限條件下, 雙態(tài)噪聲可轉(zhuǎn)化為高斯白噪聲和白散粒噪聲.這使得雙態(tài)噪聲具有一定的理論研究價值[45].更重要的是, 受乘性雙態(tài)噪聲影響的線性模型通常可以精確求解, 便于進行分析研究.在此基礎(chǔ)上, 本文通過隨機平均法推導(dǎo)出系統(tǒng)各粒子的同步性.并根據(jù)同步性求出系統(tǒng)輸出幅值增益G的表達式及穩(wěn)定的充要條件.最后通過數(shù)值仿真分析了系統(tǒng)穩(wěn)定、同步及隨機共振等集體行為.

2 模 型

本文研究受乘性雙態(tài)噪聲和余弦信號共同作用的N粒子全局耦合過阻尼系統(tǒng), 該系統(tǒng)所滿足的微分方程如下:

其中xi(t) 為第i個粒子在t時刻的位移,N為系統(tǒng)規(guī)模,為全局耦合項,ε≥0為耦合強度.系統(tǒng)的全局耦合性體現(xiàn)在系統(tǒng)中的每個粒子都與其他粒子直接相連.A0cos(?t) 為外部周期驅(qū)動力,ω≥0 表示系統(tǒng)的固有頻率.該頻率的漲落被建模為對稱雙態(tài)噪聲ξi(t).該噪聲為一個雙態(tài)Markov過程, 在±σ中隨機取值, 其均值和相關(guān)函數(shù)分別為

其中σ表示噪聲幅值,λ表示噪聲相關(guān)率,τc=1/λ為平均等待時間.

由于噪聲ξi(t) 隨機波動, 各粒子所處的勢場隨機地在

之間切換.這類隨機勢場在很多類型的系統(tǒng)中均有發(fā)現(xiàn).例如生物體內(nèi)的ATP-ADP循環(huán)勢場, ATP將能量儲存在化學(xué)鍵中, 通過水解反應(yīng), 化學(xué)鍵中的一些電子進入到更低的能量態(tài).這一過程產(chǎn)生了ADP及生物所需的能量.反過來, ADP又可以通過消耗一定的能量生成ATP.這一循環(huán)過程產(chǎn)生了生物系統(tǒng)中的ATP-ADP漲落勢場[46?48].

另一方面, 在液體環(huán)境中, 粒子往往表現(xiàn)為過阻尼運動, 即其慣性項可被忽略[49].因此, 與文獻[3]類似, 模型(2)可被看作一個受到周期驅(qū)動的耦合蛋白質(zhì)分子馬達在液體環(huán)境下ATP-ADP勢場中的運動模型.文獻[28]在局部耦合情況下對這類模型的動力學(xué)行為進行了分析, 并列舉了josephson結(jié)構(gòu)陣列[50]、神經(jīng)網(wǎng)絡(luò)[51]及耦合蛋白子馬達[52,53]作為其潛在應(yīng)用場景.另外需要指出, 若令系統(tǒng)規(guī)模N=1 , 則模型(2)退化為文獻[34]中的單粒子運動模型;N=2 時, 模型(2)即變?yōu)槲覀冎拔恼轮醒芯康膭恿W(xué)模型[54]; 若令ξ1(t)=···=ξN(t), 則本文的模型又可退化為之前Yang等[3]研究的全局耦合模型.因此, 從理論和應(yīng)用兩方面考慮, 模型(2)都具有較高的研究價值.

在本文的推導(dǎo)過程中將用到指數(shù)關(guān)聯(lián)噪聲所滿足的Shapiro-Loginov(S-L)公式[55]的推廣版本,為此, 先對S-L公式做簡要介紹, 再用該公式推導(dǎo)出本文所需的公式版本.

定理1(Shapiro-Loginov公式)η(t) 是一個噪聲.假設(shè)如下條件成立:

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

(2)ψ(t) 是指數(shù)相關(guān)噪聲η(t) 的泛函,

則有如下公式成立:

為了對模型(2)進行求解, 以下將(5)式推廣到多個相互獨立指數(shù)關(guān)聯(lián)噪聲的情形.

定理2(Shapiro-Loginov公式)ηi(t),(i=1···k)為k個相互獨立的噪聲.假設(shè)如下條件成立:

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

(2)ψ(t) 是指數(shù)相關(guān)噪聲ηi(t) 的泛函,

則有如下公式成立:

證明顯然,η1(t)η2(t)···ηk(t) 仍然是一個噪聲(即隨機過程).以下證明該噪聲滿足定理1的兩個條件.由ηi(t),(i=1···k) 的相互獨立性及條件(1)可得:

故η1(t)η2(t)···ηk(t) 仍然是一個指數(shù)關(guān)聯(lián)的隨機噪聲.由條件(2)可知,ψ(t) 是噪聲η1(t)η2(t)···ηk(t)的泛函.

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

由(6)式直接可得:

其中

為集合{ξi1,ξi2,···,ξik}所包含元素的個數(shù).1≤i1,i2,···,ik≤N,i1,i2,···,ik互不相等.

3 理論結(jié)果

3.1 同步性

本節(jié)主要研究上述耦合系統(tǒng)中各粒子行為是否具有同步性, 由于各粒子所受驅(qū)動噪聲相互獨立, 因而, 這里所說的同步性只能是統(tǒng)計意義下的同步.為此, 設(shè)法構(gòu)造包含變量〈xi〉,i=1,2,···,N的封閉微分方程組.用 1,ξi1ξi2···ξik分別乘以(2)式兩端, 因為形如ξi1ξi2···ξik這種k個互不相同的噪聲相乘的形式的因子共有個.當(dāng)k取遍1到N, 再算上常數(shù)項因子 1 , 共有 2N個線性無關(guān)的乘因子.因而可以得到 2NN個隨機微分方程, 利用S-L公式取平均, 即得到包含 2NN個變量的2NN個線性無關(guān)的非齊次常微分方程組.未知量如表1所列.

3.1.1 變量分組

接下來對表1中的變量進行分組, 將相等的變量分到一組, 以約簡變量個數(shù).各組變量間的等價關(guān)系會在下一步予以證明.

表1 變量表Table 1.Variable table.

1) 由于模型(2)中各粒子位置對稱, 表1中第一行中的各變量應(yīng)該是兩兩相等的.同理, 最后一行中的各變量也應(yīng)該是兩兩相等的.分別用U0,UN表示兩行元素構(gòu)成的集合.即:

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中的變量應(yīng)該是不相等(同步)的.

3.1.2 組內(nèi)變量等價性證明

通過上述步驟, 將表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,這些集合的總元素個數(shù)為 2NN?2N個.顯然這2NN?2N個元素線性無關(guān).下面構(gòu)造一個包含這2NN?2N個變量的 2NN?2N個線性無關(guān)的齊次常微分方程組.

對于新得到的 2NN?2N個變量中任意一個變量〈ξm1ξm2···ξmkxi〉?〈ξn1ξn2···ξnkxj〉.要得到包含該變量的微分方程, 只需用ξm1ξm2···ξmk,ξn1ξn2···ξnk分別乘以

并兩端相減, 可得

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

對(15)式取均值并應(yīng)用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個線性無關(guān)的齊次常微分方程組.由于本文主要關(guān)注粒子在長時間范圍內(nèi)的解的狀況, 而方程組的初始值對粒子的運動狀態(tài)的影響會隨著t的增大逐漸消失.所以這里不妨假設(shè)該方程組初始條件為0.利用Laplace變換法可得上述方程組((16)式和(17)式)只有零解.由此推導(dǎo)出各集合U0,···,Uk1,Uk2,···,UN中變量的相等關(guān)系.

3.1.3 同步性結(jié)論與意義

通過上述分析, 得到了各集合U0,···,Uk1,Uk2,···,UN中變量的相等關(guān)系.也即同步性結(jié)論如下:

上述同步性結(jié)論具有以下兩條重要意義.

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

又因為系統(tǒng)的平均場

從而, 由(20)式可知

對復(fù)雜耦合系統(tǒng)集體動力學(xué)行為的研究往往圍繞其平均場行為展開, (22)式表明: 在本模型中,系統(tǒng)的平均場行為與單粒子行為具有統(tǒng)計一致性,從而, 本文可通過研究單粒子的動力學(xué)行為來開展對系統(tǒng)集體動力學(xué)行為的研究.

2) 為研究單粒子的動力學(xué)行為, 需對其一階矩〈xi〉 進行求解.同步性(19)式的另一個重要意義在于: 可極大地簡化模型的求解過程, 使對〈xi〉 的解析求解成為可能.在下一節(jié)的的分析中, 將利用該同步性((19)式)完成對〈xi〉 和相應(yīng)輸出幅值增益G的具體解析求解.

3.2 輸出幅值增益

本節(jié)利用3.1節(jié)所得結(jié)果推導(dǎo)一階矩〈xi〉 及系統(tǒng)的輸出幅值增益G的表達式.首先將可能出現(xiàn)的變量用新的符號表示如下:

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

對(2)式中第i個式子兩邊取平均, 并利用同步性得到含有變量y0的方程如下:

(24)式出現(xiàn)了新的耦合項y1,1, 為了建立包含該變量的微分方程, 用ξi1ξi2···ξik乘以第i個式子兩邊.這里i∈{i1,i2,···,ik}, 則有

對(25)式兩邊取平均并利用S-L公式(9)式得

對(26)進行整理得

由于(27)式中又出現(xiàn)了新的變量yk,2.同理, 當(dāng)i∈/{i1,i2,···,ik}時, 按照同樣的步驟, 可得如下包 含變量yk,2的微分方程式:

最后, 用ξ1ξ2···ξN乘以(2)式中第i個式子兩邊,計算可得包含變量yN的微分方程式:

所得的 2N個式子((24)式, (27)式, (28)式,(29)式)可改寫成如下矩陣形式:

其中,

其中T表示矩陣的轉(zhuǎn)置.

A為三對角矩陣, 即除了三條對角線上的元素外其余元素都為零.

對(30)式作Laplace變換, 有

不妨記通過求解線性方程組(35),可以求得向量.特別地, 有

其中

其中h0(t) 為H0(s) 的Laplace逆變換.

另一方面, 系統(tǒng)((24)式, (27)式, (28)式,(29)式)可看作一個線性時不變系統(tǒng).根據(jù)線性時不變系統(tǒng)的響應(yīng)理論, 該系統(tǒng)的輸出信號是一個與輸入信號僅在幅值和相位上有差異的余弦信號.即

其中A和?分別為〈xi〉 的振幅和相位.將(40)式作Laplace變換并對比(37)式可得

最后求得輸出振幅增益為

根據(jù)(42)式, 通過數(shù)值仿真分析不同的系統(tǒng)參數(shù)變化對G的影響, 可以判斷系統(tǒng)有無隨機共振現(xiàn)象產(chǎn)生.在此之前, 首先對系統(tǒng)(2)的穩(wěn)定性做一些簡要分析.

3.3 系統(tǒng)穩(wěn)定性

系統(tǒng)(2)穩(wěn)定等價于特征方程

的根都分布在左半復(fù)平面, 即矩陣A的特征值都分布在右半復(fù)平面.在此基礎(chǔ)上, 本節(jié)建立了系統(tǒng)(2)穩(wěn)定的一個充分條件和一個充要條件.

定理3(充分條件) 當(dāng)σ2<ω(ω+λ) 時, 系統(tǒng)(2)是穩(wěn)定的.

證明1) 首先當(dāng)σ=0 時, 由同步性(20)式可知, 系統(tǒng)(2)為含一個變量的常微分方程, 又由于ω>0, 系統(tǒng)顯然是穩(wěn)定的.

2) 當(dāng)σ>0 時, 證明矩陣A的所有特征值都為正實數(shù).以此推出系統(tǒng)(2)的穩(wěn)定性.

(i) 先對ε>0 的情況予以證明.此時, 有bi+1ci>0,1≤i≤N?1.故A相似于如下實對稱三對角矩陣[56]:

A與S有相同的特征值且都為實數(shù).下面證明當(dāng)定理條件滿足時,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) 當(dāng)ε=0 時, 矩陣A退化為分塊對角矩陣.用與(i)中同樣的步驟證明即可.

下面, 為給出系統(tǒng)(2)穩(wěn)定的充要條件, 將特征方程(43)式看成是關(guān)于σ2和s的函數(shù)來展開分析, 從而將其改記為f(σ2,s).特別地, 若將σ2看成參數(shù), 則f(σ2,s) 是關(guān)于變量s的 2N次多項式.

定理4(充要條件) 系統(tǒng)(2)穩(wěn)定當(dāng)且僅當(dāng)

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

證明因為矩陣A相似于實對稱矩陣(定理3), 故對任意固定的σ2,f(σ2,s)=0 有 2N個實根.將這 2N個實根從大到小排列為

則si(σ2)(1≤i≤2N) 為σ2的連續(xù)函數(shù)[57].顯然, 系統(tǒng)(2)穩(wěn)定等價于s1(σ2)<0.故由定理3 知, 當(dāng)0<σ2<ω(ω+λ) 時 ,s1(σ2)<0.另一方面, 隨著σ的增加, 系統(tǒng)勢函數(shù)的不穩(wěn)定性逐漸增強; 從而, 當(dāng)σ充分大時, 系統(tǒng)必然是不穩(wěn)定的, 此時s1(σ2)>0.所以根據(jù)s1(σ2) 的連續(xù)性可知, 存在σ>0 , 使得s1(σ2)=0.記

由于本文只關(guān)心系統(tǒng)穩(wěn)定時的各種動力學(xué)行為, 所以后續(xù)對同步行為與隨機共振行為的分析都是在穩(wěn)定性條件(48)式成立的前提下進行的.

4 數(shù)值仿真算法

本節(jié)利用隨機Taylor展開算法對模型(2)進行數(shù)值模擬, 以此驗證算法的有效性.對于充分小的時間步長 ?t, 模型(2)在離散時間下的近似表達式為

其中

其中, ?t=tn?tn?1為時間步長, ?Xi為單個時間步長內(nèi)雙態(tài)噪聲的增量, 可通過文獻[58]中的方法仿真求得.

為了驗證仿真算法的有效性, 需要將仿真結(jié)果與理論結(jié)果進行對比.首先, 進行K次Monte Carlo實驗并記第i個粒子第k次Monte Carlo實驗的軌跡為則可得到第i個粒子位移一階矩〈xi(t)〉的數(shù)值近似為

之后, 根據(jù)(40)式可知, 輸出信號〈xi(t)〉 為與輸入信號A0cos(?t) 有相同頻率的余弦信號.所以對〈xi(t)〉num作Fourier變換即可得到〈xi(t)〉 的仿真幅值A(chǔ)num.最后利用(42)式即得輸出幅值增益的數(shù)值仿真結(jié)果Gnum.圖1給出了仿真算法的一個具體的例子.其中, 圖1的第一排子圖為10個粒子進行 1 04次Monte Carlo仿真所得平均位移, 各粒子的初始位移服從標(biāo)準(zhǔn)正態(tài)分布.對該平均位移作Fourier變換得到第二排子圖對應(yīng)的頻域值.由此得到在該組參數(shù)下輸出幅值增益的數(shù)值仿真結(jié)果Gnum=1.11708 與理論值G=1.11677 的相對誤差僅為 2.78×10?4.

圖1 粒子平均位移及對應(yīng)頻域值的數(shù)值仿真結(jié)果, 其中ε=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仿真次數(shù)K和總仿真時長T對數(shù)值算法收斂性的影響.圖2(a)—圖2(c)分別描繪了Gnum作為上述三個仿真參數(shù)的函數(shù)的曲線圖.可以看出, 隨著仿真?t的減小及K和T的增加, 仿真結(jié)果逐漸趨近于理論結(jié)果.且當(dāng) ?t≤0.001,K≥10000 ,T≥120時, 仿真結(jié)果與理論結(jié)果的誤差小于 6×10?4, 對應(yīng)的相對誤差僅為 5.372×10?4.因此, 若不特別指出, 本文數(shù)值仿真結(jié)果都為固定參數(shù)?t=0.001,K=10000,T=120 下仿真求得.

圖2 仿真結(jié)果與理論結(jié)果對比圖, 其中?=π/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 集體動力學(xué)行為分析

5.1 系統(tǒng)穩(wěn)定性

在3.3節(jié)中, 推導(dǎo)得出了粒子位移一階矩穩(wěn)定的充要條件.這里進一步通過仿真分析穩(wěn)定性區(qū)域隨參數(shù)的變化及其對粒子運動行為的影響.圖3(a)給出了不同的N所對應(yīng)的系統(tǒng)穩(wěn)定性區(qū)域圖像.其中黑色區(qū)域為N=2 所對應(yīng)的穩(wěn)定性區(qū)域; 黑色區(qū)域與深灰色區(qū)域的并為N=5 所對應(yīng)的系統(tǒng)穩(wěn)定性區(qū)域; 黑色區(qū)域, 深灰色區(qū)域與淺灰色區(qū)域的并集則為N=10 時的穩(wěn)定性區(qū)域.從中可以看出,系統(tǒng)的穩(wěn)定性區(qū)域隨耦合強度ε及系統(tǒng)規(guī)模N的增加而增大.為了更直觀地反映系統(tǒng)穩(wěn)定性區(qū)域隨N的變化, 圖3(b)給出了在一組固定的參數(shù)下,系統(tǒng)穩(wěn)定性區(qū)域邊界隨N的變化曲線.可以看出,穩(wěn)定性區(qū)域隨著N的增大單調(diào)遞增.關(guān)于圖3中現(xiàn)象的一個合理的物理解釋為: 粒子規(guī)模N及耦合強度ε決定粒子間的耦合力Fi, 而耦合力為系統(tǒng)輸出提供有序性.噪聲為粒子位移引入無序性, 且噪聲幅值σ越大, 粒子位移的無序性越強.有序性和無序性相互競爭決定系統(tǒng)輸出的穩(wěn)定性.當(dāng)σ∈[0,σmin)時, 耦合力提供的有序性占主導(dǎo)地位.當(dāng)σ=σmin時, 有序與無序的競爭達到一個臨界點,使得初值的影響不會隨時間增加而消失.系統(tǒng)輸出開始變得不穩(wěn)定.當(dāng)σ>σmin時, 噪聲引入的無序性占主導(dǎo)地位, 系統(tǒng)輸出無界.由于增加N及ε可以增加耦合力Fi, 使得系統(tǒng)輸出的有序性增強.此時就需要更大的σ才能使得系統(tǒng)輸出達到穩(wěn)定性的臨界點.所以, 系統(tǒng)的穩(wěn)定性區(qū)域會隨著粒子規(guī)模N及耦合強度ε的增加而增大.

圖3 ( a) 不同的N所對應(yīng)的系統(tǒng)穩(wěn)定性區(qū)域, 其中ω=0.5,λ=1 ; ( b) 系統(tǒng)穩(wěn)定性區(qū)域邊界隨N的變化曲線, 其中 ε =1 , 其他參 數(shù)同圖(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).

為了直觀反應(yīng)穩(wěn)定性區(qū)域?qū)αＷ悠骄灰频挠绊?不失一般性, 這里固定N=10 并從圖3(a)中選擇A(1.5,1),B(2,1) 兩點進行仿真實現(xiàn).結(jié)果如圖4所示.從圖4(a)可以看出, 粒子平均位移有界, 此時系統(tǒng)是穩(wěn)定的.圖4(b)中的系統(tǒng)輸出隨著時間的增大向無窮遠發(fā)散, 此時系統(tǒng)是不穩(wěn)定的.這與定理4的結(jié)論相符.

圖4 粒子平均位移實現(xiàn), 其中ω=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.

由于本文只關(guān)心系統(tǒng)穩(wěn)定時的各種動力學(xué)行為, 所以以下兩節(jié)的仿真都是在穩(wěn)定性條件成立的前提下進行的.

5.2 粒子的集體同步行為

本節(jié)通過數(shù)值仿真對3.1節(jié)中的結(jié)論(20)式予以驗證.即隨著系統(tǒng)時間的演化,N個粒子的平均位移會達到同步狀態(tài).

圖5給出了系統(tǒng)(2)在不同參數(shù)條件下不同粒子平均位移的仿真結(jié)果和對應(yīng)的方差, 其中各粒子的初始位置服從標(biāo)準(zhǔn)正態(tài)分布.從圖5(a)—圖5(c)的第一行子圖中可以看出, 具有不同初始值的粒子隨著時間的增加最終都會達到同步狀態(tài).這與3.1節(jié)中的結(jié)論相符, 即方程組(16)式, (17)式的初值條件不改變長時間后系統(tǒng)中粒子的同步狀態(tài).為直觀反映粒子的同步速度, 引入粒子平均位移〈xi(t)〉num的方差?(t) 如下:

其中

圖5(a)—圖5(c)的第二行子圖給出了方差隨時間變化的曲線圖.從中可直觀地看出, 圖5(a)中曲線的下降速度慢于圖5(b)和圖5(c).且當(dāng)t=1 時,三幅圖對應(yīng)的方差?分別是4.21×10?3,9.2×10?7,1.2×10?7.下面通過定義同步時間t0, 對不同參數(shù)下系統(tǒng)(2)的同步速度進行定量描述.由于系統(tǒng)(2)中各噪聲ξi獨立同分布, 且實驗仿真次數(shù)有限, 所以?隨著時間增加不會降為零.因此, 這里假設(shè)?(t)≤10?5時, 即認(rèn)為系統(tǒng)達到同步狀態(tài), 并記相應(yīng)的同步時間為:

將圖5(a)—圖5(c)第二行子圖中各方差曲線圖橫縱坐標(biāo)以對數(shù)形式表示得到第三行子圖.三幅子圖中系統(tǒng)達到同步的時間分別是 3.53,0.72,0.55.對比三幅圖中的同步時間可知, 增加耦合強度ε和系統(tǒng)規(guī)模N可顯著加快粒子的同步速度.這是由于一方面增加耦合強度ε使得粒子間的耦合力

圖5 系統(tǒng)(2)在不同參數(shù)條件下的仿真實現(xiàn)和所對應(yīng)的方差, 其中ω=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.

加大, 而耦合力的作用是使得各粒子相互吸引, 這使得系統(tǒng)可以更快同步.另一方面由于系統(tǒng)中的耦合關(guān)系為全局耦合, 增加N也會加大粒子間的吸引力, 從而讓各粒子更快達到同步狀態(tài).另外, 從圖6也可以定量地看出, 隨著N的增加, 系統(tǒng)的同步時間t0逐漸變小.

圖6 同步時間 t0 隨群體規(guī)模N的變化曲線, 所選系統(tǒng)參數(shù)與圖5相同F(xiàn)ig.6.The change curve of synchronization time t0 with system size N.The selected system parameters are the same as the Fig.5.

最后需要指出, 耦合強度只改變系統(tǒng)的同步時間而不改變系統(tǒng)達到穩(wěn)定時的同步狀態(tài).這是由于系統(tǒng)中各粒子的位移主要受周期外力、粒子所處勢場、各粒子的初始位置以及粒子間的耦合力四個因素影響.這四個因素中, 決定粒子同步的本質(zhì)因素為周期外力和粒子所處勢場.一方面, 各粒子受相同的周期外力作用; 另一方面, 各粒子有相同的本征頻率, 頻率噪聲有相同的統(tǒng)計性質(zhì), 這使得各粒子所處的勢場從統(tǒng)計意義看是相同的.所以一段時間后, 系統(tǒng)初值對各粒子位移的影響消失.耦合強度的存在只是加速了粒子的同步速度, 當(dāng)ε=0 時各粒子所滿足的運動學(xué)方程相同, 此時各粒子仍然會同步.故耦合強度不改變系統(tǒng)的同步性.

5.3 隨機共振

隨機共振現(xiàn)象是噪聲的隨機性與驅(qū)動力的有序性在非線性作用下相互協(xié)作的結(jié)果.在本文的模型中, 隨機性由噪聲提供, 有序性由正弦驅(qū)動力和耦合作用力兩部分組成.由于正弦驅(qū)動力的作用在以往的文獻中已展開過大量充分的討論, 因而本文主要分析噪聲與耦合力的協(xié)作現(xiàn)象.它們相互協(xié)作產(chǎn)生隨機共振的機理如下: 噪聲ξi(i=1,2,···,N)為系統(tǒng)提供隨機性, 噪聲幅值σ越大, 則隨機性越強; 耦合力Fi(i=1,2,···,N) 為系統(tǒng)提供有序性,耦合強度ε或系統(tǒng)規(guī)模N越大, 則有序性越強.在“乘性”這種非線性作用方式下, 有序與無序相互競爭, 當(dāng)兩者間達到某種最優(yōu)匹配時, 系統(tǒng)產(chǎn)生隨機共振, 從而使得系統(tǒng)輸出達到最大.

在上述共振機理的指導(dǎo)下, 接下來通過數(shù)值分析的方式, 逐一對輸出幅值增益G關(guān)于耦合強度ε、系統(tǒng)規(guī)模N以及噪聲幅值σ的共振行為進行具體分析.

5.3.1 關(guān)于耦合強度ε的參數(shù)隨機共振

圖7(a)和圖7(b)分別給出了不同的σ,N下,G隨ε的變化曲線.可以看出, 除了σ=0 外, 其他曲線都有不同強度的共振峰.且共振峰主要出現(xiàn)在低耦合區(qū)域.隨著耦合強度ε的增大, 每條曲線都趨于一個定值.這是因為當(dāng)參數(shù)σ,N固定時, 系統(tǒng)輸出有序性主要由ε確定.且系統(tǒng)輸出有序性隨ε增加而增強.當(dāng)ε過小時, 過弱的耦合力不能有效協(xié)調(diào)粒子的運動.當(dāng)ε過大時, 過強的耦合力會使粒子相互吸引組成一個整體而完全抑制噪聲帶來的隨機性的影響.兩種情況都限制了噪聲與輸入信號的協(xié)作效應(yīng).因此, 使得輸出幅值增益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)中, 當(dāng)σ=0 時,G取常值.這是由于當(dāng)σ=0時, 噪聲消失.此時模型(2)退化為確定性模型, 統(tǒng)計意義下的同步〈x1〉=〈x2〉=···=〈xN〉 變?yōu)閤1=x2=···=xN.耦 合 力Fi(i=1,2,···,N)使各粒子相互吸引組成一個整體, 之后變?yōu)榱?因此,ε對G不產(chǎn)生影響.另一方面, 將σ=0 代入(42)式直接計算可得:

當(dāng)σ≥0 時, 共振峰的高度隨著σ的增大而增大.這說明增加噪聲幅值σ可加大隨機共振的強度.

圖7(b)中, 隨著系統(tǒng)規(guī)模N的增加, 共振峰變得越來越尖銳, 同時共振峰的位置左移, 共振峰的高度不斷增加.產(chǎn)生這種現(xiàn)象的原因如下: 由于系統(tǒng)中的耦合關(guān)系為全局耦合, 增加N會加大粒子間的吸引力, 所以當(dāng)N較大時, 隨著ε的增加,粒子間的吸引力增大得更快.因此G隨ε的變化曲線也更陡.同樣地, 更強的耦合力進一步抑制了噪聲帶來的隨機性, 使得G不能達到最優(yōu).所以當(dāng)噪聲幅值σ固定時, 增加N會使得共振峰左移.另外注意到, 圖7(b)中所有曲線都有共同的起點.這是由于當(dāng)粒子間的耦合作用消失時(ε=0 ), 模型(2)變?yōu)镹個彼此獨立且相同的系統(tǒng).此時, 改變系統(tǒng)規(guī)模N并不改變系統(tǒng)中各粒子的受力狀況.因此N的變化不影響G.同樣地, 將ε=0 代入(42)式可得:

所以, 從(58)式也可得出上述結(jié)論.

5.3.2 關(guān)于系統(tǒng)規(guī)模N的參數(shù)隨機共振

圖8(a)和圖8(b)分別給出了不同的σ,ε下,G隨N的變化曲線.可以看出, 除了σ=0,ε=0外, 其他曲線都有不同強度的單峰共振.這表明在合適的參數(shù)下, 存在最優(yōu)的系統(tǒng)規(guī)模N使得輸出幅值增益G達到最優(yōu).這里系統(tǒng)中有序與無序競爭關(guān)系的變化主要由系統(tǒng)規(guī)模N的變化引起.在給定參數(shù)ε,σ的條件下, 系統(tǒng)規(guī)模N的增加使粒子間的耦合力Fi(i=1,2,···,N) 變大, 進而使得系統(tǒng)輸出有序逐漸增大.當(dāng)系統(tǒng)規(guī)模N較小時, 過弱的耦合力不能有效協(xié)調(diào)粒子的運動, 即噪聲引入的隨機性過大.當(dāng)系統(tǒng)規(guī)模N較大時, 過強的耦合力會完全抑制噪聲帶來的隨機性的影響, 使系統(tǒng)(2)近似于確定性系統(tǒng).這兩種極端情況都會限制粒子的平均運動.因此, 使得輸出幅值增益G最大的N一定是介于兩者之間的.

圖8(a)中, 當(dāng)ε=0 時,G取常值.當(dāng)ε>0 時,隨著耦合強度ε的增加, 使共振曲線趨于穩(wěn)定值的N減小.這與圖7(b)中的結(jié)論一致.另外還可以看出, 隨著ε的增加, 共振曲線下移.這表明增加耦合強度ε會減弱系統(tǒng)關(guān)于N的參數(shù)隨機共振強度.

圖8(b)中, 當(dāng)σ=0 時,G取常值.這是由于此時模型(2)退化為確定性模型, 而耦合力Fi(i=1,2,···,N)使各粒子相互吸引組成一個整體, 之后變?yōu)榱?此時, 不論粒子個數(shù)多少, 系統(tǒng)中的各粒子都只受周期外力A0cos(?t) 的作用.所以改變系統(tǒng)規(guī)模N不影響輸出幅值增益.該結(jié)論也與(57)式一致.當(dāng)σ>0 時, 隨著σ的增加, 共振峰右移, 共振強度加大.這是由于更大的σ會使粒子的運動軌跡產(chǎn)生更強的隨機性, 過強的隨機性又會抑制系統(tǒng)的輸出.這種情況下, 要使系統(tǒng)達到最優(yōu)輸出就需要加強粒子間的吸引力使粒子運動軌跡更加有序.所以當(dāng)σ增加時, 使系統(tǒng)輸出達到最大的N增大.另外還可以看出, 增加噪聲幅值σ可以加強系統(tǒng)的參數(shù)隨機共振強度.

圖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 關(guān)于噪聲幅值σ的隨機共振

圖9(a)和圖9(b)分別給出了不同的ε,N下,G隨σ的變化曲線.可以看出, 所有曲線都表現(xiàn)出不同強度的共振現(xiàn)象.這表明, 在合適的參數(shù)下,系統(tǒng)關(guān)于噪聲幅值σ會產(chǎn)生隨機共振.從圖9(a)和圖9(b)還可以看出, 當(dāng)σ=0 時, 所有曲線都有共同的起點.這一現(xiàn)象與圖7(a), 圖8(b)及(57)式得到的結(jié)論一致.

圖9(a)中, 隨著耦合強度ε從0增加到1, 共振強度變大且共振峰位置右移.當(dāng)ε從1繼續(xù)增加到10時, 共振強度逐漸減弱, 并且也逐漸趨向于一個極限值, 而共振峰的位置相對變化較小.這說明增加耦合強度既能增強也可以削弱隨機共振現(xiàn)象.

圖9(b)中, 隨著系統(tǒng)規(guī)模N的增加, 共振強度變大且共振峰位置右移.出現(xiàn)這一現(xiàn)象的原因與圖8(b)類似.更大的系統(tǒng)規(guī)模N使系統(tǒng)各粒子間的相互吸引力變強, 而過強的吸引力進一步抑制了噪聲帶來的隨機性, 使系統(tǒng)(2)產(chǎn)生與確定性系統(tǒng)類似的有序輸出.因此, 隨著N的增加, 使系統(tǒng)輸出G達到最大的σ變大, 共振峰右移.圖9(b)表明, 增加系統(tǒng)規(guī)模N可增強系統(tǒng)的隨機共振強度.

圖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 結(jié) 論

本文研究了不同頻率噪聲驅(qū)動下過阻尼全局耦合系統(tǒng)的集體動力學(xué)行為.理論方面, 主要通過隨機平均法來展開研究.首先, 得到了各粒子行為的統(tǒng)計同步性, 該同步性使得系統(tǒng)平均場行為與單粒子行為具有統(tǒng)計一致性, 從而可通過對單粒子行為的研究來完成對平均場行為(也即集體行為)的研究.其次, 利用該同步性, 完成了對系統(tǒng)輸出幅值增益的解析求解, 為進一步分析系統(tǒng)的隨機共振行為奠定了理論基礎(chǔ).最后, 推導(dǎo)給出了系統(tǒng)穩(wěn)定的充分條件和充要條件, 為本文所得結(jié)論的適用范圍給出了參考.數(shù)值仿真方面, 主要通過隨機Taylor展開算法進行研究.首先, 分析了粒子規(guī)模N及耦合強度ε對穩(wěn)定性區(qū)域與同步時間的影響, 結(jié)果表明, 隨著耦合強度ε的增加或系統(tǒng)規(guī)模N的增大,粒子間的耦合力增大, 系統(tǒng)有序性增強, 從而使得穩(wěn)定區(qū)域逐漸增大, 同步時間逐漸縮短.其次, 對系統(tǒng)的隨機共振行為展開了研究.噪聲為系統(tǒng)提供隨機性, 耦合力為系統(tǒng)提供有序性, 兩者相互競爭,從而使得系統(tǒng)輸出關(guān)于噪聲強度σ、耦合強度ε與系統(tǒng)規(guī)模N均表現(xiàn)出隨機共振行為.隨著耦合強度的增加或系統(tǒng)規(guī)模的增大, 系統(tǒng)的有序性增強,需要更大的噪聲強度提供更強的隨機性來與之實現(xiàn)最優(yōu)匹配, 從而關(guān)于噪聲強度σ的共振峰逐漸右移.反之, 隨著噪聲強度σ的增大, 關(guān)于耦合強度ε與系統(tǒng)規(guī)模N的共振峰也會右移.本文的研究工作雖然圍繞特定耦合模型展開, 但相關(guān)結(jié)論具有很強的普適性, 對其他復(fù)雜耦合系統(tǒng)的研究具有很好的理論指導(dǎo)和參考意義.后續(xù), 將在本文工作的基礎(chǔ)上, 進一步結(jié)合實際需求展開深入研究.