姜 雷，彭 靜，袁曉輝

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

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

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

氣體放電模型(Gas -Discharge)如式(1)所示，最初是由PRUWINS 提出的，用來(lái)模擬氣體放電實(shí)驗(yàn)中觀察到的斑圖，現(xiàn)在經(jīng)常被用來(lái)研究耗散系統(tǒng)的動(dòng)力學(xué)特性。該模型是一個(gè)三變量的反應(yīng)擴(kuò)散方程組，其中，變量u是激勵(lì)子，變量v和w是抑制子。與包含一個(gè)激勵(lì)子和抑制子的雙參數(shù)反應(yīng)擴(kuò)散方程模型(如Gray -Scott 和FitzHugh-Nagumo 模型)不同，Gas - Discharge 模型中第二個(gè)抑制子w對(duì)二維空間中孤立子的形成起到關(guān)鍵作用。如果在雙參數(shù)的Gray -Scott 模型中加入一個(gè)抑制子，也能獲得穩(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 所示，孤立子是一個(gè)局部解。在均勻二維空間中，令方程組左邊為0，可以得到u=v=w，也就是說(shuō)二維空間中孤立子解在無(wú)限遠(yuǎn)處時(shí)u=v=w。因此，在模擬時(shí)可以認(rèn)為空間中除了孤立子附近區(qū)域，其他部分可以看作一個(gè)平面。但是在非均勻空間中，k1的值發(fā)生變化，也就導(dǎo)致在非均勻邊界處u、v、w的值發(fā)生變化，形成一個(gè)新的解，形狀如圖1(b)中L/2 附近圖1(a)的灰色線條所示，稱之為非均勻解(HIOP)。因此，空間足夠大時(shí)，孤立子在非均勻空間上的動(dòng)力學(xué)分析可以看作是其與HIOP 碰撞的動(dòng)力學(xué)研究。根據(jù)NISHIURA 的研究結(jié)果，HIOP 解還有其他4種類型［7-8］。

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

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

2 模擬結(jié)果

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

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

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

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

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

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

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

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

3 結(jié)論

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