二維非均勻空間中氣體放電模型的動力學(xué)分析

姜 雷，彭 靜，袁曉輝

(1.哈爾濱師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，黑龍江 哈爾濱150025;2.東北農(nóng)業(yè)大學(xué) 電氣與信息工程學(xué)院，黑龍江 哈爾濱150030;3.中國科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所，黑龍江 哈爾濱150040)

孤立子是自然界中一種常見的現(xiàn)象，在氣體放電試驗(yàn)［1］、Pt 表面的CO 催化氧化反應(yīng)等許多物理、化學(xué)試驗(yàn)中都可以觀察到孤立子運(yùn)動。從20 世紀(jì)60 年代開始，孤立子就引起了研究者的極大興趣。為解釋其背后的數(shù)學(xué)機(jī)理，Gas -Discharge、Gray -Scott 等模型被用來研究孤立子相互碰撞后產(chǎn)生的復(fù)雜現(xiàn)象，如反彈、融合和消失等。通過對Gray -Scott 模型孤立子解碰撞的數(shù)值解析，NISHIURA［2］等發(fā)現(xiàn)分水嶺解是控制這些復(fù)雜現(xiàn)象變遷的主要因素。分水嶺解是一種不穩(wěn)定的解，當(dāng)其受到不同方向的擾動時會有不同的變化，如反彈或者穿過、融合或者分裂等。與孤立子之間的碰撞相比，非均勻空間上孤立子的動力學(xué)特性研究相對較少。YUAN［3］等研究了一維非均勻空間中Gas-Discharge 模型的孤立子解的特性，隨著參數(shù)變化，孤立子穿過非均勻空間時會產(chǎn)生反彈、停止和消失等現(xiàn)象。與均勻空間中孤立子的碰撞類似，一維非均勻空間中也存在著一類分水嶺解控制著孤立子的狀態(tài)變遷。NISHIURA［4-5］等模擬了二維空間中孤立子從不同角度與非均勻邊界進(jìn)行碰撞的過程。在二維空間中隨著自由度的增加，細(xì)微擾動可以使孤立子的軌道發(fā)生遷移，從而觀察到旋轉(zhuǎn)、分裂等新現(xiàn)象。為了解析二維非均勻空間中孤立子的動力學(xué)特性，筆者利用Gas-Discharge 模型模擬了孤立子與非均勻邊界正面碰撞的過程，觀察到了反射、穿過和分裂等現(xiàn)象的變遷，通過對分水嶺解的穩(wěn)定性解析，解釋了這些現(xiàn)象產(chǎn)生的機(jī)理。筆者還利用Pseudo Arclength 方法獲得了參數(shù)空間中解的全局分歧圖，闡明了控制不同狀態(tài)變遷的分水嶺解之間的聯(lián)系。

1 模型與數(shù)值模擬方法

氣體放電模型(Gas -Discharge)如式(1)所示，最初是由PRUWINS 提出的，用來模擬氣體放電實(shí)驗(yàn)中觀察到的斑圖，現(xiàn)在經(jīng)常被用來研究耗散系統(tǒng)的動力學(xué)特性。該模型是一個三變量的反應(yīng)擴(kuò)散方程組，其中，變量u是激勵子，變量v和w是抑制子。與包含一個激勵子和抑制子的雙參數(shù)反應(yīng)擴(kuò)散方程模型(如Gray -Scott 和FitzHugh-Nagumo 模型)不同，Gas - Discharge 模型中第二個抑制子w對二維空間中孤立子的形成起到關(guān)鍵作用。如果在雙參數(shù)的Gray -Scott 模型中加入一個抑制子，也能獲得穩(wěn)定的孤立子解［6］。

式中：Δ 為拉普拉斯項(xiàng);u=u(t;r);v=v(t;r);w=w(t;r);r=(x;y)∈R2，k2=2.0;k3=1.0;k4=8.5;τ= 40。擴(kuò)散系數(shù)Du=0.9 ×10-4;Dv=1.0 ×10-3;Dw=1.0 ×10-2。參數(shù)k1的值在x方向上按照式(2)所示變化，其中代表左邊的k1值，kR1代表右邊的k1值，ε=-。

在Gas-Discharge 模型中，非均勻空間模擬示意圖如圖1 所示，孤立子是一個局部解。在均勻二維空間中，令方程組左邊為0，可以得到u=v=w，也就是說二維空間中孤立子解在無限遠(yuǎn)處時u=v=w。因此，在模擬時可以認(rèn)為空間中除了孤立子附近區(qū)域，其他部分可以看作一個平面。但是在非均勻空間中，k1的值發(fā)生變化，也就導(dǎo)致在非均勻邊界處u、v、w的值發(fā)生變化，形成一個新的解，形狀如圖1(b)中L/2 附近圖1(a)的灰色線條所示，稱之為非均勻解(HIOP)。因此，空間足夠大時，孤立子在非均勻空間上的動力學(xué)分析可以看作是其與HIOP 碰撞的動力學(xué)研究。根據(jù)NISHIURA 的研究結(jié)果，HIOP 解還有其他4種類型［7-8］。

圖1 非均勻空間模擬示意圖

二維空間中Gas-Discharge 模型的數(shù)值計(jì)算量較大，筆者根據(jù)模型的特點(diǎn)編寫了并行分析程序。模型的模擬采用二維Crank -Nicolson 方法。為了獲得孤立子和HIOP 的數(shù)值解，首先給定一個與HIOP 解非常接近的初值，然后利用Newton- Raphson 迭代法獲得精確解。當(dāng)獲得一個HIOP 解后，可以運(yùn)用Pseudo - Arclength 方法繪制HIOP 解的全局分歧圖。由于求解過程中涉及大量矩陣運(yùn)算，可采用Aztec 軟件包將矩陣運(yùn)算并行化，從而提高程序的運(yùn)行效率，最后采用Arpack 軟件包求解HIOP 解的特征值和特征向量。

2 模擬結(jié)果

圖2 孤立子運(yùn)動相圖

運(yùn)用上述數(shù)值模擬方法，在= (- 7. 5，-6.7)，=(-8.0，-6.7)的參數(shù)空間中，孤立子與HIOP 碰撞后可以觀察到反彈(REB)，穿過(PEN)，以及分裂(SPL)這3 種現(xiàn)象(如圖2 所示)。反彈(REB)是指孤立子從左往右的運(yùn)行過程中，與HIOP 碰撞以后，改變方向從右向左運(yùn)動。從圖2 中可以看到有兩個反彈區(qū)域，這兩個區(qū)域上孤立子返回的位置有所不同，在左上角的反彈區(qū)域，孤立子在HIOP 左側(cè)返回;而在右下角的反彈區(qū)域，孤立子在HIOP 的右側(cè)返回。這也暗示著相圖中反彈(REB)與穿過(PEN)的變遷由兩種不同的機(jī)制控制。分裂是孤立子與HIOP碰撞后變?yōu)閮蓚€孤立子沿著HIOP 反向移動。與反彈類似，兩個分裂區(qū)域分別代表在HIOP 左側(cè)和右側(cè)分裂。在一維非均勻空間中，YUAN 發(fā)現(xiàn)分水嶺解控制孤立子從反彈到穿過之間的變遷。在二維空間中，是否也存在這樣的分水嶺解問題，筆者仔細(xì)模擬了邊界上A、B、C、D點(diǎn)附近孤立子的運(yùn)動情況。在A點(diǎn)附近(，)≈(-7.10，-7.42)，當(dāng)= -7.423 55 時，孤立子穿過HIOP 到達(dá)右側(cè)區(qū)域，而當(dāng)= -7.423 56 時，孤立子與HIOP 碰撞后反彈回左側(cè)區(qū)域。從圖3 的時間演化圖上可以看出，在孤立子決定反彈還是穿過之前，會在一個狀態(tài)停留一段時間。這似乎存在一個靜止的不穩(wěn)定解，當(dāng)受到擾動時，其會選擇運(yùn)動方向。利用牛頓迭代法，可以獲得該狀態(tài)的數(shù)值解，如圖3 所示。這個解與孤立子解不同，首先其速度為0，是一個靜止解。其次其與均勻空間中的靜止解不同，由于參數(shù)的空間非一致性，解在x方向上是不對稱的，在y方向上對稱。實(shí)際上，這個解也是NISHIURA 描述的HIOP 解中的一種。計(jì)算該解的特征值可以得到兩個大于零的特征值，λ1=7.652 3 ×10-2，λ2=2.635 4 ×10-2，其對應(yīng)的特征向量如圖3 所示，分別對應(yīng)x方向和y方向。在一維空間中，只能得到一個大于0 的特征值，特征向量也只對應(yīng)一個方向。將特征向量的值作為擾動加到靜止解上，可以看到靜止解馬上沿著x方向繼續(xù)前進(jìn)，表現(xiàn)為穿過，如果將特征向量反向加到靜止解上，則會反向運(yùn)動，表現(xiàn)為反彈。在y方向上特征值作為擾動也會使孤立子沿著y方向正向或反向移動。由于圖1 所示的模擬中，孤立子與HIOP 正向碰撞，碰撞過程中在y方向的變形是對稱的，因此x方向上的變形對孤立子的運(yùn)動狀態(tài)起到主要影響作用。A點(diǎn)附近的HIOP 解在此起到了分水嶺解的作用。

圖3 控制REB 和PEN 的分水嶺解

在B附近點(diǎn)可以觀察到與A點(diǎn)類似的轉(zhuǎn)變，當(dāng)= -6.985 04時，孤立子穿過HIOP 到達(dá)右側(cè)區(qū)域，而當(dāng)=-6.985 03時，孤立子與HIOP 碰撞后反彈回左側(cè)區(qū)域。用牛頓迭代法得到一個不穩(wěn)定解，計(jì)算得到該解有兩個正的特征值，λ1=9.563 4 ×10-2，λ2=3.678 2 ×10-2，分別對應(yīng)兩個x，y方向的特征向量φ。因此，從上述數(shù)值計(jì)算可以知道，靜止的非穩(wěn)定解是控制穿過與反彈現(xiàn)象的分水嶺解［9-10］。

相圖中有兩個分裂區(qū)域，分裂現(xiàn)象是一維空間沒有觀察到的新現(xiàn)象。仔細(xì)觀察穿過與分裂邊界上C點(diǎn)附近孤立子的演化發(fā)現(xiàn)，如圖4 所示，孤立子與HIOP 碰撞后繼續(xù)前進(jìn)，碰撞后的孤立子變?yōu)橐粋€具有雙峰形狀的解，而雙峰解并不能穩(wěn)定地繼續(xù)前進(jìn)，當(dāng)= -6.757 170 311 時，雙峰解變?yōu)閮蓚€沿著y方向反向運(yùn)動的孤立子，當(dāng)= -6.757 170 312時，雙峰解變?yōu)橐粋€沿著x方向繼續(xù)前進(jìn)的孤立子。利用牛頓迭代法可以得到一個不穩(wěn)定的雙峰解，有一個大于0 的特征值λ =0.656 5。將其對應(yīng)的特征值加到不穩(wěn)定解雙峰解上可以觀察到雙峰解融合為一個孤立子，將特制值反向加到雙峰解上會分裂為兩個孤立子。因此在穿過與分裂的變化中，不穩(wěn)定的雙峰解充當(dāng)了分水嶺解的角色。在D點(diǎn)附近孤立子從REB 變遷到SPL，同樣也是由一個不穩(wěn)定的雙峰解控制。C、D點(diǎn)附近分水嶺解的區(qū)別在于：C點(diǎn)的分水嶺解位于HIOP 的左側(cè)，而D點(diǎn)的分水嶺解位于HIOP 的右側(cè)。

圖4 C 點(diǎn)附近孤立子的時間演化圖

從模擬的結(jié)果可以看到，在二維空間中分水嶺解可以是不穩(wěn)定的單峰解，也可以是不穩(wěn)定的雙峰解。與均勻空間中的分水嶺解不同，非均勻空間中分水嶺解在x方向是非對稱的。分水嶺解是控制孤立子碰撞后狀態(tài)變遷的關(guān)鍵因素，因此，如果能夠獲得分水嶺解的全局分歧圖，根據(jù)分水嶺解的不穩(wěn)定特征，就可以知道在參數(shù)空間內(nèi)將會發(fā)生哪些現(xiàn)象。利用Pseduo -Arclength 方法，從已知的單峰和雙峰分水嶺解出發(fā)可以得到的 全 局 分 歧 圖(圖5)。圖5 中，實(shí)線代表雙峰解(SP)，虛線代表單峰解(SS)。從圖5 中可以看出雙峰解和單峰解并不是單獨(dú)的分支，通過解的拓展可以相互連接起來。在相圖中A-D點(diǎn)附近的4 種分水嶺解都是相互關(guān)聯(lián)的。還可以看到，除了4 種類型的分水嶺解，還存在小的雙峰解和單峰解，這些解也都是不穩(wěn)定解。當(dāng)孤立子與HIOP 正面碰撞以后其形狀會向分水嶺解接近，隨著參數(shù)的變化，接近的軌道不同導(dǎo)致了各種現(xiàn)象的分化。在圖5 中，虛線代表時均勻空間上的雙峰解和單峰解，這也說明非均勻空間的分水嶺解和均勻空間的分水嶺解可以通過解的拓展連接起來。

圖5 kL1 = -7.0 的全局分歧圖

3 結(jié)論

通過對Gas-Discharge 模型在二維非均勻空間上孤立子解的動力學(xué)分析，揭示了孤立子狀態(tài)變遷的機(jī)理。在二維非均勻空間中，Gas - Discharge 模型的非均勻參數(shù)引起了多種形態(tài)的HIOP 解，其中4 種不穩(wěn)定的HIOP 解充當(dāng)了分水嶺解的角色，這4 種解可以通過參數(shù)空間中的解拓展相互聯(lián)系。HIOP 解的全局分歧圖為解釋孤立子的動力學(xué)特性提供了依據(jù)。

［1］PRUWINS S C，OR-GUIL M，BODE M，et al. Interacting pulses in three -component reaction -diffusion systems on two -dimensional domains［J］. Phys Rev Lett，1997(78)：3781 -3784.

［2］NISHIURA Y，TERAMOTO T，UEDA K I. Dynamic transitions through scattors in dissipative systems［J］.Chaos，2003，13(12)：962 -972.

［3］YUAN X，TERAMOTO T，NISHIURA Y. Heterogeneity-induced defect bifurcation and pulse dynamics for a three - component reaction - diffusion system［J］.Phys Rev E，2007，75(3)：622 -630.

［4］NISHIURA Y，TERAMOTO T，YUAN X. Heterogeneity-induced spot dynamics for a three -component reaction - diffusion system［J］. Communications on Pure and Applied Analysis，2012(11)：307 -338.

［5］NISHIURA Y，UEDA K I. A mathematical mechanism for instabilities in stripe formation on growing domains［J］. Physica D，2011(241)：37 -59.

［6］NISHIURA Y，TERAMOTO T，UEDA K I. Scattering and separators in dissipative systems［J］. Phys Rev E，2003，71(5)：621 -630.

［7］NISHIURA Y，TERAMOTO T，YUAN X，et al. Dynamics of traveling pulses in heterogeneous media［J］.Chaos，2007，17(3)：71 -75.

［8］NISHIURA Y，TERAMOTO T，UEDA K I. Scattering of traveling spots in dissipative systems［J］. Chaos，2005，15(4)：75 -84.

［9］TERAMOTO T，SUZUKI K，NISHIURA Y. Rotational motion of traveling spots in dissipative systems［J］.Phys Rev E，2009，77(4)：620 -628.

［10］TERAMOTO T，YUAN X，BAR M，et al. Onset of unidirectional pulse propagation in an excitable medium with asymmetric heterogeneity［J］. Phys Rev E，2009，77(4)：602 -607.