王艷麗, 李向紅,,3, 王 敏, 申永軍,3
(1.石家莊鐵道大學(xué) 機(jī)械工程學(xué)院,石家莊 050043;2. 石家莊鐵道大學(xué) 數(shù)理系,石家莊 050043;3.石家莊鐵道大學(xué) 交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國(guó)家重點(diǎn)實(shí)驗(yàn)室,石家莊 050043)
多尺度耦合現(xiàn)象已經(jīng)廣泛存在于機(jī)械、化工、生物等各個(gè)領(lǐng)域,成為目前非線性動(dòng)力學(xué)領(lǐng)域的研究熱點(diǎn)和前沿課題之一,其中簇發(fā)振動(dòng)[1]是多尺度耦合系統(tǒng)的典型動(dòng)力學(xué)行為,表現(xiàn)為在一個(gè)周期內(nèi)系統(tǒng)交替出現(xiàn)沉寂態(tài)和激發(fā)態(tài)。文獻(xiàn)[2-4]基于多尺度效應(yīng)和非線性動(dòng)力學(xué)穩(wěn)定性與分岔理論,深入研究了神經(jīng)系統(tǒng)中神經(jīng)元放電活動(dòng)的分岔模式。文獻(xiàn)[5-6]在非光滑系統(tǒng)、周期激勵(lì)系統(tǒng)、化學(xué)反應(yīng)振蕩系統(tǒng)的多尺度耦合效應(yīng)方面做了大量的工作,比如利用微分包含理論深入研究了非光滑系統(tǒng)尺度效應(yīng)的分岔機(jī)制,提出了轉(zhuǎn)換相圖和包絡(luò)分析等方法,揭示了慢變激勵(lì)下的簇發(fā)振蕩模式及其演化過(guò)程。此外,國(guó)外很多學(xué)者也對(duì)簇發(fā)振動(dòng)行為做了一些研究。Kouayep等[7]研究了一種光電振子,在電路中引入Colpitts振子會(huì)產(chǎn)生高頻電信號(hào),并表現(xiàn)出周期性和簇發(fā)性等動(dòng)力學(xué)行為。Herve等[8]研究了單擺臂和雙擺臂通過(guò)磁場(chǎng)耦合到非線性電路振動(dòng)臺(tái)上的響應(yīng),結(jié)果表明由于電子信號(hào)的傳遞,擺臂的位移會(huì)發(fā)生簇發(fā)振動(dòng),并且振動(dòng)的形狀周期和振幅取決于各種控制參數(shù)。然而,目前的研究主要集中在整數(shù)階系統(tǒng)上,相比較下,分?jǐn)?shù)階系統(tǒng)的快慢效應(yīng)研究成果較少。
近年來(lái),分?jǐn)?shù)階微積分系統(tǒng)被廣泛的應(yīng)用到高能物理、反常擴(kuò)散、復(fù)雜黏彈性材料力學(xué)本構(gòu)關(guān)系、系統(tǒng)控制等若干領(lǐng)域[9],其理論研究和應(yīng)用研究受到了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注。比如,王軍等[10-12]研究了一些含有分?jǐn)?shù)階微分的非線性與線性振子,采用平均法得到了系統(tǒng)的一階近似解,證明了分?jǐn)?shù)階微分項(xiàng)不僅起著阻尼的作用同時(shí)還起著剛度的作用。Xu等[13]提出了具有多個(gè)時(shí)變時(shí)滯環(huán)節(jié)的分?jǐn)?shù)階競(jìng)爭(zhēng)神經(jīng)網(wǎng)絡(luò),并研究了這類神經(jīng)網(wǎng)絡(luò)的全局漸近穩(wěn)定性,為揭示穩(wěn)定性判據(jù)與網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)之間的密切關(guān)系提供了一個(gè)新視角。Kandasamy等[14]研究了分?jǐn)?shù)階脈沖四元數(shù)值神經(jīng)網(wǎng)絡(luò)的多重Mittag-Leffler穩(wěn)定性問(wèn)題并仿真驗(yàn)證了所得到的主要結(jié)果的有效性。Sekerci[15]研究了不同分?jǐn)?shù)階定義下的食餌-捕食者系統(tǒng)的區(qū)別與聯(lián)系。Sllva等[16]研究了兩個(gè)基于分?jǐn)?shù)階微積分提出的常微分方程系統(tǒng),解決了腫瘤細(xì)胞和免疫系統(tǒng)之間的動(dòng)力學(xué)問(wèn)題。
Brusselator模型[17]是一類描述自催化反應(yīng)的典型振子,早在1968年由Prigogine等提出。由于催化劑的存在,使反應(yīng)過(guò)程中不同反應(yīng)物的反應(yīng)速率存在較大的差別,即涉及快慢不同的時(shí)間尺度。在之后的幾十年里,越來(lái)越多的學(xué)者投入到對(duì)Brusselator振子分岔、混沌以及近似解的研究中[18-24]。Rech[25]研究了周期強(qiáng)迫的四參數(shù)Brusselator振子,利用四維參數(shù)空間的一個(gè)給定截面,證明了該系統(tǒng)存在多穩(wěn)態(tài)。Capone等[26]研究了全局Brusselator振子解的有界性,分析了離散化對(duì)線性穩(wěn)定性的影響,發(fā)現(xiàn)在常微分方程系統(tǒng)中的穩(wěn)定性與有擾動(dòng)下的穩(wěn)定性截然不同,并進(jìn)行了非線性穩(wěn)定性分析。Wang等[27]在研究分?jǐn)?shù)階微分系統(tǒng)的同步現(xiàn)象時(shí)首次提出了分?jǐn)?shù)階Brusselator系統(tǒng),發(fā)現(xiàn)當(dāng)分?jǐn)?shù)階階次總和小于0.97時(shí),分?jǐn)?shù)階Brusselator振子存在一個(gè)極限環(huán)。Jena等[28]采用分?jǐn)?shù)階約化差分變換法求解分?jǐn)?shù)階Brusselator振子的近似解,并對(duì)其進(jìn)行了誤差分析,發(fā)現(xiàn)該方法能夠以級(jí)數(shù)形式快速收斂。目前圍繞Brusselator振子整數(shù)階成果較多,相比之下,其關(guān)于多尺度耦合的分?jǐn)?shù)階系統(tǒng)方面的研究尚在逐步發(fā)展階段。
本文主要考慮外部周期擾動(dòng)因素影響反應(yīng)過(guò)程中催化劑濃度的情形下,分?jǐn)?shù)階Brusselator振子的簇發(fā)行為及其分岔機(jī)制,研究不同參數(shù)對(duì)系統(tǒng)的影響規(guī)律,為實(shí)際工程領(lǐng)域提供理論指導(dǎo)。
分?jǐn)?shù)階Brusselator振子
(1)
(2)
式中:n-1 實(shí)際上,催化過(guò)程易受到外界因素干擾及影響,因此要考慮外在周期因素。本文從外部擾動(dòng)影響反應(yīng)過(guò)程中催化劑濃度的情形出發(fā),將這種擾動(dòng)因素近似為外部周期擾動(dòng),從而進(jìn)行動(dòng)力學(xué)分析,其模型為 (3) 式中:γ和ω分別為外部周期擾動(dòng)的幅值和頻率,令θ=ωt,則式(3)可轉(zhuǎn)化為廣義自治系統(tǒng) (4) 當(dāng)ω遠(yuǎn)小于原系統(tǒng)的固有頻率,式(3)會(huì)存在量級(jí)上差異的兩個(gè)時(shí)間尺度,使其非線性行為更為復(fù)雜。此時(shí)由于周期性擾動(dòng),該反應(yīng)會(huì)變?yōu)榭炻齼蓚€(gè)時(shí)間尺度耦合的系統(tǒng),在式(4)中變量x和y為快子系統(tǒng),θ為慢子系統(tǒng)?;诳炻治龇ǎ瑢⒙兞恳暈榭熳酉到y(tǒng)的分岔參數(shù)。令F=γcosθ為快子系統(tǒng)的慢變參數(shù),則快子系統(tǒng)表示為 (5) 由 (6) 得到快子系統(tǒng)式(5)的平衡點(diǎn) 下面將按照A+F是否等于0分兩種情況進(jìn)行討論。 平衡點(diǎn)E0處的特征方程為 λ2+[(A+F)2+1-B]λ+(A+F)2=0 (7) 式(5)在平衡點(diǎn)E0處的Jacobian矩陣為 (8) 進(jìn)一步分析分?jǐn)?shù)階系統(tǒng)發(fā)生Hopf分岔的條件,選取系統(tǒng)參數(shù)F為分岔參數(shù),定義函數(shù) 根據(jù)整數(shù)階系統(tǒng)發(fā)生Hopf分岔的條件,進(jìn)而推廣到分?jǐn)?shù)階系統(tǒng)發(fā)生Hopf分岔[29]的條件,見(jiàn)命題1。 命題1當(dāng)B≥1,且 若滿足下列條件: 條件1式(5)對(duì)應(yīng)的特征方程有一對(duì)復(fù)共軛根λ1,2(F)=θ(F)±iω(F); 下面證明此命題: λ1,2=θ(F)±iω(F) (9) 條件1得證。 (2) 在命題所給條件下,特征方程僅有一對(duì)共軛復(fù)根,那么有 因此, (11) 基于分?jǐn)?shù)階微分方程穩(wěn)定性理論[30],可得式(5)失穩(wěn)的邊界條件 (12) 具體表達(dá)式為 (13) 條件2得證。 (3) 由于 (15) 兩邊關(guān)于F求導(dǎo),可得 其中, θ(F)=[B-1-(A+F)2]/2, θ′(F)=-(A+F), (17) 因?yàn)锳+F≠0,所以 θ′(F*)ω(F*)-ω′(F*)θ(F*)= (18) 條件3得證。綜上,命題1成立。 選取參數(shù)A=1,B=3.8,可得到快子系統(tǒng)的階次q和慢變參數(shù)F之間的雙參數(shù)分岔圖。如圖1(a)所示,區(qū)域(I)、區(qū)域(II)和區(qū)域(III)內(nèi)平衡點(diǎn)是穩(wěn)定的,當(dāng)跨越臨界線進(jìn)入?yún)^(qū)域(IV)時(shí),平衡點(diǎn)失穩(wěn)同時(shí)出現(xiàn)穩(wěn)定的極限環(huán),圖1(b)和圖1(c)分別給出了區(qū)域(II)和區(qū)域(IV)中的系統(tǒng)時(shí)間歷程圖。圖1(a)中A和B表示q=0.75時(shí)對(duì)應(yīng)的兩個(gè)臨界點(diǎn)。可見(jiàn),式(5)有兩個(gè)Hopf分岔點(diǎn),并且隨著分?jǐn)?shù)階階次的減小,兩個(gè)Hopf分岔點(diǎn)之間的距離越來(lái)越小。 當(dāng)q=1時(shí),圖2給出了快子系統(tǒng)式(5)關(guān)于慢變參數(shù)F的分岔圖,其中平衡線F1-F2有五種類型的平衡點(diǎn),分別是位于F1-H1上的穩(wěn)定焦點(diǎn),H1-S1上的不穩(wěn)定焦點(diǎn),S1-F0,F(xiàn)0-S2上為不穩(wěn)定結(jié)點(diǎn)(這里不包括F0點(diǎn),即F=-1的情況),S2-H2上的不穩(wěn)定焦點(diǎn),H2-F2上的穩(wěn)定焦點(diǎn)。其中H1和H2為Hopf分岔臨界點(diǎn),參數(shù)值為FH1=-2.67,F(xiàn)H2=0.67,當(dāng)F在H1-F0和F0-H2(這里不包括F0點(diǎn)的情況)快子系統(tǒng)存在穩(wěn)定的極限環(huán),圖3給出了F為-1.5和-0.5時(shí)的相圖。 圖1 系統(tǒng)的分岔圖和時(shí)間歷程圖Fig.1 Bifurcation diagram and time history diagram of the system 圖2 快子系統(tǒng)的平衡點(diǎn)關(guān)于慢變參數(shù)F的分岔圖Fig.2 The bifurcation diagram with respect to parameter F 圖3 極限環(huán)的相圖Fig.3 The phase diagram of limit cycle 對(duì)于整數(shù)階的情形,若A+F=0,此時(shí)F=-1。對(duì)于式(5),直線x=0上任何點(diǎn)均為系統(tǒng)的奇點(diǎn),即F=-1時(shí),縱坐標(biāo)y軸是奇線。顯然,式(5)的平衡點(diǎn)為(0,m),m為任意實(shí)數(shù)。在(0,m)處對(duì)系統(tǒng)線性化可得矩陣 (19) 特征值是0和-4.8,此時(shí)不能直接判別平衡點(diǎn)的穩(wěn)定性,需要借助中心流形定理進(jìn)行判斷。構(gòu)造一個(gè)T矩陣,其列為M的特征向量,即 (20) 此時(shí),顯然不是標(biāo)準(zhǔn)形,因此需進(jìn)一步標(biāo)準(zhǔn)化來(lái)判斷其平衡點(diǎn)的穩(wěn)定性,故有 (21) 應(yīng)用變量代換 (22) 可得 (23) 在所給參數(shù)條件下整數(shù)階Brusselator式(5)可變?yōu)?/p> (24) 將式(23)代入式(24)中可以進(jìn)一步推出 (25) 此時(shí), Ec={(x1,x2)∈R2x2=0},x軸 (26a) Es={(x1,x2)∈R2x1=0},y軸 (26b) 從而可知式(25)是標(biāo)準(zhǔn)形式,因此 (27) (28) 中心流形方程為 h′(x1){-18.24[h(x1)]3-4.8x1[h(x1)]2}+ 18.24[h(x1)]3+4.8x1[h(x1)]2+4.8[h(x1)]=0 (29) 圖4 奇線Fig.4 Singular line 對(duì)于F=-1時(shí),分?jǐn)?shù)階的穩(wěn)定性情況,可以利用數(shù)值模擬進(jìn)行驗(yàn)證。圖5是對(duì)應(yīng)分?jǐn)?shù)階情況下,F(xiàn)=-1時(shí)的時(shí)間歷程圖??梢园l(fā)現(xiàn),對(duì)于分?jǐn)?shù)階系統(tǒng),與整數(shù)階系統(tǒng)有相同的結(jié)論,即在F=-1處,分?jǐn)?shù)階系統(tǒng)也存在穩(wěn)定的奇線x=0。 圖5 當(dāng)F=-1時(shí)的時(shí)間歷程圖Fig.5 The time history for F=-1 本部分研究系統(tǒng)簇發(fā)振動(dòng)行為與機(jī)理,以及不同參數(shù)對(duì)簇發(fā)振動(dòng)的影響。 圖6是參數(shù)γ=0.8,ω=0.03,q=1時(shí)的時(shí)間歷程圖??梢杂?jì)算出,整個(gè)周期振動(dòng)頻率恰好是外部周期擾動(dòng)頻率0.03。為了進(jìn)一步揭示其簇發(fā)振動(dòng)的產(chǎn)生機(jī)理,給出了關(guān)于變量x與慢過(guò)程F=0.8cos(0.03t)的轉(zhuǎn)換相圖,并將其轉(zhuǎn)換相圖和分岔圖疊加得到圖7。 圖7進(jìn)一步表明了簇發(fā)振動(dòng)的產(chǎn)生機(jī)理。假設(shè)在一個(gè)周期的演變過(guò)程中,相軌跡從A點(diǎn)出發(fā),在A點(diǎn)時(shí)平衡點(diǎn)是不穩(wěn)定的,并受到快子系統(tǒng)穩(wěn)定極限環(huán)的吸引,導(dǎo)致系統(tǒng)呈現(xiàn)出大幅度振動(dòng)而處于激發(fā)態(tài)。極限環(huán)會(huì)在經(jīng)過(guò)Hopf分叉點(diǎn)H2后消失,只有穩(wěn)定的平衡點(diǎn)吸引子,系統(tǒng)的大幅振動(dòng)將逐漸收斂至穩(wěn)定平衡線,此時(shí)系統(tǒng)完成了激發(fā)態(tài)演變到沉寂態(tài)的過(guò)程。最后當(dāng)軌線到達(dá)B時(shí),慢變過(guò)程F=0.8cos(0.03t)達(dá)到最大值,軌線將按照F減小的方向運(yùn)動(dòng)。此時(shí)系統(tǒng)仍處于沉寂態(tài),并且在經(jīng)過(guò)分岔點(diǎn)H2,系統(tǒng)并不會(huì)立刻進(jìn)入激發(fā)態(tài),而是由于分岔延遲在點(diǎn)P,開(kāi)始大幅振動(dòng),系統(tǒng)再次處于激發(fā)態(tài),直至完成了一個(gè)周期的運(yùn)動(dòng),該簇發(fā)屬于Hopf簇發(fā)。同時(shí),圖8(a)和圖8(b)分別給出了周期激勵(lì)增加和減少兩個(gè)方向上激發(fā)態(tài)滯后現(xiàn)象。 圖6 當(dāng)q=1時(shí)的時(shí)間歷程圖Fig.6 The time history for q=1 圖7 當(dāng)q=1時(shí),轉(zhuǎn)換相圖和分岔圖的疊加Fig.7 Overlapping of transition portrait with bifurcation diagram for q=1 圖8 當(dāng)q=1時(shí),簇發(fā)現(xiàn)象的產(chǎn)生機(jī)理Fig.8 The generation mechanism of the periodic bursting oscillation for q=1 下面討論分?jǐn)?shù)階階次對(duì)簇發(fā)振動(dòng)的影響,固定參數(shù)γ=0.8,ω=0.03。圖9和圖10分別給出了q=0.75和q=1.15時(shí)的時(shí)間歷程圖。 圖9 當(dāng)q=0.75時(shí)的時(shí)間歷程圖Fig.9 The time history for q=0.75 圖10 當(dāng)q=1.15時(shí)的時(shí)間歷程圖Fig.10 The time history for q=1.15 通過(guò)對(duì)比圖9、圖6和圖10可以看出,系統(tǒng)的激發(fā)態(tài)持續(xù)時(shí)間會(huì)隨著分?jǐn)?shù)階階次的增大而變長(zhǎng),相對(duì)應(yīng)沉寂態(tài)時(shí)間會(huì)縮短。該結(jié)論也可通過(guò)圖1(a)中快子系統(tǒng)的臨界階次和慢變參數(shù)之間的關(guān)系圖得到,即階次越小,兩個(gè)Hopf分岔點(diǎn)之間的距離越小,從而導(dǎo)致激發(fā)態(tài)的持續(xù)時(shí)間變短。同時(shí),系統(tǒng)幅值也會(huì)隨著分?jǐn)?shù)階階次減小而降低。 為了研究周期擾動(dòng)幅值對(duì)簇發(fā)振動(dòng)的影響,固定q=0.75,ω=0.03,通過(guò)改變?chǔ)玫闹祦?lái)具體研究幅值的影響。由圖1(a)可知,當(dāng)q=0.75時(shí),兩個(gè)Hopf分岔點(diǎn)A,B的參數(shù)值分別為FA=-2.33,F(xiàn)B=0.33,這說(shuō)明,當(dāng)FA 當(dāng)擾動(dòng)幅值γ>2.33,F(xiàn)的變化范圍包含[-2.33,0.33],存在三種吸引子:位于中間的穩(wěn)定極限環(huán)和位于兩邊的穩(wěn)定焦點(diǎn)吸引子,這將導(dǎo)致系統(tǒng)的軌跡產(chǎn)生變化。其中,中間與左右兩端的吸引子的雙穩(wěn)性導(dǎo)致系統(tǒng)軌線在激發(fā)態(tài)與沉寂態(tài)之間相互轉(zhuǎn)遷,這樣會(huì)出現(xiàn)兩種激發(fā)態(tài)和兩種沉寂態(tài)的雙Hopf簇發(fā)現(xiàn)象,圖11給出了當(dāng)γ=2.5時(shí)的時(shí)間歷程圖。 圖11 當(dāng)γ=2.5時(shí)的時(shí)間歷程圖Fig.11 The time history for γ=2.5 如果擾動(dòng)幅值0.33<γ<2.33,整個(gè)系統(tǒng)只涉及快子系統(tǒng)的兩種吸引子,即右側(cè)的穩(wěn)定焦點(diǎn)吸引子和中間的極限環(huán)吸引子。此時(shí),系統(tǒng)只存在單側(cè)雙穩(wěn)性。如圖9所示,為單Hopf簇發(fā)。 如果擾動(dòng)幅值γ<0.33,此時(shí)系統(tǒng)軌線只涉及到穩(wěn)定極限環(huán)吸引子,雙穩(wěn)性演變?yōu)閱畏€(wěn)性,沉寂態(tài)消失,兩個(gè)激發(fā)態(tài)連接在一起,圖12給出了此種情況γ=0.3時(shí)的時(shí)間歷程圖。 圖12 當(dāng)γ=0.3時(shí)的時(shí)間歷程圖Fig.12 The time history for γ=0.3 由此可見(jiàn),在快慢耦合的系統(tǒng)中,如果周期擾動(dòng)因素為慢變過(guò)程,則擾動(dòng)幅值的大小可以調(diào)節(jié)系統(tǒng)的振動(dòng)模式,其機(jī)理是涉及系統(tǒng)吸引子的種類發(fā)生了改變。 本文研究了具有外部周期擾動(dòng)的分?jǐn)?shù)階Brusselator振子,當(dāng)外部周期擾動(dòng)頻率遠(yuǎn)小于系統(tǒng)固有頻率時(shí),系統(tǒng)涉及兩個(gè)不同的時(shí)間尺度,表現(xiàn)出明顯的快慢效應(yīng),存在典型的簇發(fā)振動(dòng)。通過(guò)系統(tǒng)的穩(wěn)定性分析發(fā)現(xiàn),系統(tǒng)存在一條穩(wěn)定的奇線,并利用中心流形定理和數(shù)值模擬證明了其穩(wěn)定性。在非奇線處,根據(jù)分?jǐn)?shù)階系統(tǒng)穩(wěn)定性理論討論了系統(tǒng)的穩(wěn)定性,并給出了發(fā)生Hopf分岔的充分條件。 分?jǐn)?shù)階階次和幅值對(duì)簇發(fā)振動(dòng)具有明顯的影響。隨著分?jǐn)?shù)階階次的減小,激發(fā)態(tài)逐漸減弱,其原因是由于隨著參數(shù)變化,兩個(gè)Hopf分岔點(diǎn)之間的距離變小,導(dǎo)致系統(tǒng)處于激發(fā)態(tài)的時(shí)間縮短。當(dāng)擾動(dòng)幅值較大時(shí),系統(tǒng)存在雙Hopf周期簇發(fā)振動(dòng)現(xiàn)象,其產(chǎn)生機(jī)理是系統(tǒng)涉及快子系統(tǒng)的兩個(gè)Hopf分岔,導(dǎo)致了激發(fā)態(tài)與沉寂態(tài)之間相互轉(zhuǎn)遷的分岔行為;隨著擾動(dòng)幅值的減小,雙Hopf簇發(fā)逐漸演變?yōu)閱蜨opf簇發(fā);當(dāng)幅值繼續(xù)減小時(shí),系統(tǒng)不再涉及Hopf分岔點(diǎn),且僅涉及穩(wěn)定極限環(huán)吸引子,單穩(wěn)性會(huì)使系統(tǒng)一直處于激發(fā)態(tài),沉寂態(tài)消失,但是整個(gè)周期振動(dòng)是兩個(gè)頻率的耦合。 以上結(jié)論為研究分?jǐn)?shù)階Brusselator振子的振動(dòng)與控制提供了重要的理論依據(jù)和技術(shù)服務(wù),同時(shí)整數(shù)階與分?jǐn)?shù)階系統(tǒng)行為的本質(zhì)差異問(wèn)題需要進(jìn)一步深入研究。1.1 情況一:A+F≠0
1.2 情況二:A+F=0
2 簇發(fā)振動(dòng)行為及參數(shù)的影響
2.1 簇發(fā)振動(dòng)及其機(jī)理
2.2 分?jǐn)?shù)階階次對(duì)簇發(fā)振動(dòng)的影響
2.3 周期擾動(dòng)幅值對(duì)簇發(fā)振動(dòng)的影響
3 結(jié) 論