基于分裂格式有限點集法對孤立波二維非線性問題的模擬*

任金蓮 任恒飛 陸偉剛 蔣濤

(揚州大學(xué),數(shù)學(xué)科學(xué)學(xué)院,水利與能源動力學(xué)院,揚州 225002)

在提出一種基于時間分裂格式的純無網(wǎng)格有限點集(split-step finite pointset method,SS-FPM)法的基礎(chǔ)上,數(shù)值模擬了含孤立波的二維非線性薛定諤 (nonlinear Schr?dinger,NLS)/ (Gross-Pitaevskii,GP)方程.SS-FPM的構(gòu)造過程為: 1)基于時間分裂的思想將非線性薛定諤方程分成線性導(dǎo)數(shù)項和非線性項; 2)采用基于Taylor展開和加權(quán)最小二乘法的有限點集法,借助Wendland權(quán)函數(shù),對線性導(dǎo)數(shù)項進(jìn)行數(shù)值離散.隨后,模擬了帶有Dirichlet和周期性邊界條件的NLS方程,將所得結(jié)果與解析解做對比.數(shù)值結(jié)果表明: 給出的SS-FPM粒子法的優(yōu)點是在粒子分布非均勻情況下仍具有近似二階精度,且較網(wǎng)格類有限差分算法實施容易,較已有改進(jìn)的光滑粒子動力學(xué)方法計算誤差小.最后,運用SS-FPM對無解析解的二維周期性邊界NLS方程和Dirichlet邊界玻色-愛因斯坦凝聚二分量GP方程進(jìn)行了數(shù)值預(yù)測,并與其他數(shù)值結(jié)果進(jìn)行對比,準(zhǔn)確展現(xiàn)了非線性孤立波奇異性現(xiàn)象和量子化渦旋過程.

1 引 言

含孤立波解的非線性問題常見于非線性光學(xué)和晶體的熱脈沖物理現(xiàn)象[1-3],以及玻色-愛因斯坦凝聚態(tài) (BEC)動力學(xué)特性[3-6]中.該非線性物理現(xiàn)象的研究通常會涉及非線性薛定諤(nonlinear Schr?dinger,NLS)方程或 (Gross-Pitaevskii,GP)方程的求解[5-7],然而NLS/GP方程中含有的非線性項或旋轉(zhuǎn)角動量算子使得多分量或高維情況下的孤立波形解難以用解析手段精確地得到[7,8].目前已提出了多種數(shù)值方法對NLSE/GPE進(jìn)行求解或?qū)?fù)雜孤立波傳播過程進(jìn)行預(yù)測,如有限元法[9]、有限差分法[10,11]、蒙特卡羅法[12,13]、修正的歐拉算法[14]、時間分裂譜方法[11,15]和基于背景網(wǎng)格無網(wǎng)格方法[16-18]等.然而上述基于網(wǎng)格的方法在多維復(fù)雜區(qū)域或非均勻節(jié)點分布情況下編程模擬實現(xiàn)都較復(fù)雜.近些年來,純無網(wǎng)格方法(粒子方法)以其完全不依賴于網(wǎng)格的優(yōu)勢在實數(shù)域偏微分方程的模擬中得到了許多應(yīng)用,如光滑粒子動力學(xué)方法(smoothed particle hydrodynamics,SPH)方法[19-23]和有限點集法 (finite pointset method,FPM)[24-26]等.但上述純無網(wǎng)格方法對復(fù)數(shù)域上含孤立波解非線性問題的模擬研究還處于起步階段,特別是FPM方法對非線性薛定諤方程的模擬在國際上還鮮有研究.

不依賴于背景網(wǎng)格的FPM方法具有粒子方法的特點,其在非線性薛定諤方程的模擬應(yīng)用上還處于試探性研究階段,這與其在數(shù)值穩(wěn)定和數(shù)值精度方面還待進(jìn)一步研究相關(guān).FPM方法精度和穩(wěn)定性的提高,以及非線性薛定諤方程的數(shù)值模擬研究均是國際上的兩個研究熱點.FPM方法模擬復(fù)雜非線性薛定諤方程較網(wǎng)格類方法的主要優(yōu)點在于: 可以任意布點,不受區(qū)域復(fù)雜性限制; 計算空間導(dǎo)數(shù)可采用二階精度顯格式且不依賴于網(wǎng)格; 模擬程序易于實現(xiàn)且便于并行實施.

本文針對非線性薛定諤方程的特點,為提高直接推廣FPM方法模擬帶非線性項問題的精度和穩(wěn)定性,首先將非線性薛定諤方程進(jìn)行時間分裂分為非線性項和線性導(dǎo)數(shù)項; 其次,引入Wendland權(quán)函數(shù)[26],采用基于Taylor展開和移動最小二乘思想的顯式FPM格式對線性導(dǎo)數(shù)部分進(jìn)行離散;從而得到一種能夠準(zhǔn)確、高效地求解非線性薛定諤方程的基于時間分裂有限點集法 (split-step finite pointset method,SS-FPM).所得數(shù)值結(jié)果表明,提出的SS-FPM能夠有效地求解帶有Dirichlet和周期性邊界條件的二維非線性薛定諤方程且具有二階精度,并能夠準(zhǔn)確預(yù)測出周期邊界條件下孤立波變化奇異現(xiàn)象和Dirichlet邊界條件下帶外旋轉(zhuǎn)BEC二分量的孤立波值隨時間變化的量子化渦旋過程.

2 Gross-Pitaevskii方程

Gross-Pitaevskii方程常被用來描述量子力學(xué)中BEC動力學(xué)特性[1,11],本文考慮如下帶角動量旋轉(zhuǎn)項的無量綱化GP方程:

初值條件

邊值條件

或周期邊值條件

其中空間變量x=(x,y),或(x,y,z);ψ(x,t)是d維空間中的復(fù)值波函數(shù);Δ是d維空間的拉普拉斯算子,是虛數(shù)單位;φ(x)是一個給定初值的復(fù)值函數(shù).無量綱實數(shù)βd通常用來描述三次非線性項的相互強(qiáng)度|ψ|2ψ,Lz=-i(x?y-y?x)是無量綱旋轉(zhuǎn)角速度θ下的角動量算子的z分量,勢Vd(x,t)是一個實值函數(shù),其形式為

其中γx,γy,γz為實常數(shù);Vd還可以是周期性的函數(shù).

3 SS-FPM

本文針對超冷原子BEC動力學(xué)特性問題的粒子法模擬,將時間分裂格式與FPM進(jìn)行耦合,得到能夠準(zhǔn)確模擬GP方程的SS-FPM.本文提出的SS-FPM方法較直接拓展應(yīng)用FPM方法[24,25]模擬GP方程具有較高精度和更好的長時間模擬穩(wěn)定性,這是由于長時間模擬NLS/GP方程時,方程中的非線性項對算法精度和穩(wěn)定性的要求較高.本文提出的SS-FPM方法的思想是: 首先,應(yīng)用時間分裂法將方程分解成帶非線性和線性導(dǎo)數(shù)項的兩個方程; 其次,采用具有較好穩(wěn)定性的Wendland權(quán)函數(shù)[23]拓展應(yīng)用FPM法對帶線性導(dǎo)數(shù)方程進(jìn)行二階顯格式離散; 最后,采用二階時間分裂偽譜法對兩個方程進(jìn)行交替求解.

3.1 時間分裂格式

引入文獻(xiàn)[7,10,15]中時間分裂(time-split,TS)法,GP方程(1)可被寫為

其中A=-(1/2)Δ-θLz是線性微分算子,B=Vd(x,t)+ βd|ψ(x,t)|2是非線性算子.于是上述方程可以被分解為如下的線性方程

和非線性方程

本文采用的二階分裂步法[7]: 首先求解(8)式,其次以(8)式的解作為初始條件求解(7)式,最后以(7)式的解作為初始條件求解(8)式,從而得到下一時間層的解.

3.2 FPM離散格式

對方程(7)的求解采用FPM法,時間上采用二階龍格-庫塔離散格式,函數(shù)導(dǎo)數(shù)項近似采用顯式FPM離散格式,FPM格式求解函數(shù)一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的基本思想[24,25]是: 第一步,將求解域任意離散成有限點并賦初始值; 第二步,每個點x在權(quán)函數(shù)支持域內(nèi)有相鄰粒子xj(j=1,2,···,n,n為支持域內(nèi)相鄰粒子數(shù)),每個xj在x上進(jìn)行Taylor展開并保留到二階導(dǎo)數(shù)項; 第三步,對Taylor級數(shù)余項進(jìn)行整理,引入最小二乘法得到一個以x處一階/二階導(dǎo)數(shù)為未知函數(shù)的線性方程組; 第四步,求解涉及局部矩陣的線性方程組得到每個點x處函數(shù)導(dǎo)數(shù)近似值; 第五步,采用二階龍格-庫塔時間離散格式得到下一個時間層函數(shù)值.

本文以二維空間上均勻布點為例,采用具有較好穩(wěn)定性Wendland權(quán)函數(shù)[23,27].Wendland權(quán)函數(shù)[23,27]形式如下:

其中 r=|xj-x| 為 支持(域半)徑; ω0是一個正常數(shù),ω0在二維空間中為 7 /64πh2,h為光滑長度,此處取 h ≈ 1.1λ0( λ0為分布粒子初始間距).

對任意點x處函數(shù) ψ (x,t),其支持域內(nèi)相鄰粒子為xj,函數(shù) ψ (xj,t)在x點的Taylor級數(shù)展開

其中ej為Taylor級數(shù)展開時的余項誤差.將(10)式右端第一項放到等式左邊,則(10)式可以寫成

其中

引入加權(quán)最小二乘法[24]可以得到

(12)式可以寫為

其中

根據(jù)極值原理,J取最小可得

通過(13)式包含5 × 5局部系數(shù)矩陣的線性方程組求解可得x處一、二階導(dǎo)數(shù)值.

3.3 SS-FPM離散格式

結(jié)合3.1節(jié)二階時間分裂格式與3.2節(jié)FPM離散格式,對GP方程(1)進(jìn)行離散,可得如下SSFPM離散格式.

對非線性方程(8)通過求解常微分方程的方式進(jìn)行求解[11],

將(14)式求解得到的結(jié)果代入有限點集法中進(jìn)行計算,再將求得的結(jié)果代入線性方程中,得到

最后將方程(15)得到的解代入非線性方程(8)中進(jìn)行求解,得到

其中 j=0,1,2,···,N ,N是區(qū)域上離散化的粒子數(shù); m=0,1,2,···為時間層;(15)式通過3.2節(jié)FPM離散格式進(jìn)行求解得到

時間步長為dt,為保證數(shù)值模擬穩(wěn)定性,時間步長的選擇通常需要滿足限制性條件 (見文獻(xiàn)[19,24,25]),本文選取為粒子初始間距).

初邊值條件施加:在 (14)式和(15)式計算中,初邊值條件的施加也是至關(guān)重要的,初始條件(2)式可以在計算 (14)式前進(jìn)行準(zhǔn)確施加=φ(x).邊值條件(3)式采用文獻(xiàn)[1,7,10]處理方式,將其近似為齊次Dirichlet邊值條件施加( ? Ω 為區(qū)域邊界).周期邊值條件(4)式的施加中為保證邊界上粒子的不足,采用SPH粒子方法的施加方式 (見文獻(xiàn)[19-21]),在物理量值更新前需要施加周期性條件的邊界上取大于等于支持域尺寸2h的粒子數(shù)及物理量賦到相對的邊界上.

4 NLS/GP方程SS-FPM模擬

為驗證提出的SS-FPM法模擬NLS方程的準(zhǔn)確性和預(yù)測無解析解NLS/GP問題孤立波隨時間演化特性的可靠性,本節(jié)首先通過對帶兩種邊值條件的NLS方程的求解,與解析解做比較,對SSFPM法數(shù)值精度和收斂速度進(jìn)行分析; 其次運用提出的粒子方法對周期邊界條件下孤立波奇異特性和BEC中孤立波量子化渦旋過程進(jìn)行數(shù)值預(yù)測,并與其他數(shù)值結(jié)果(SS-FDM (split-step finite difference method)法[4,7]和SS-ICPSPH (split-step implicit corrected parallel smoothed particle hydrodynamics)法)[20]進(jìn)行對比.為分析所提方法的精度和收斂性,本文定義如下誤差 (精確解與數(shù)值解最大誤差范數(shù))和收斂階為:

這里的 λ01和 λ02分別表示不同的粒子初始間距.

4.1 SS-FPM方法精度和收斂性分析

運用SS-FPM對兩個不同邊值條件下NLS方程進(jìn)行求解,分析了所提方法模擬NLS方程的精度和收斂速度,討論了粒子分布非均勻情況下的數(shù)值誤差.

4.1.1 Dirichlet邊值NLS方程

考慮正方形區(qū)域 Ω :[0,2π]× [0,2π]的Dirichlet邊值條件的NLS方程[7,23],其對應(yīng)的方程和初邊值條件分別為

其中V(x,y)=1-sin2xsin2y.

初值條件為

邊值條件為

對應(yīng)該問題的解析解為

圖1 幾個不同時刻處沿 y=0.5π 的 | ψ| 變化曲線 (a)粒子均勻分布; (b)粒子非均勻分布Fig.1.The change curve of | ψ| along y=0.5π at different time: (a)Uniform mode; (b)non-uniform mode.

圖2 兩種不同的粒子分布 (a)均勻粒子分布; (b)非均勻粒子分布Fig.2.Two kinds of particle distribution: (a)Uniform mode; (b)non-uniform mode.

該算例模擬中采用粒子初始間距為 π /64 ,時間步長為dt=10—4.圖1給出了幾個時刻兩種粒子分布情況下SS-FPM法得到沿 y=0.5π 處波函數(shù)變化曲線,并與解析解進(jìn)行對比.圖2展示了本算例模擬中采用的均勻和非均勻兩種粒子分布情況,其中圖2(a)的粒子是均勻分布方式,分別沿x,y方向每隔 π /64 的距離分布一個粒子; 圖2(b)中的粒子是非均勻分布方式,以 ( π,π)為圓心,靠近圓心第一層布6個粒子,第二層布12個粒子,逐步成等差數(shù)列方式向外擴(kuò)展沿圓形分布粒子,相鄰兩個圓形層之間距離為相等均為 π /64 ,但邊界上仍采用均勻布點方式.由圖1可知,隨時間演化提出的粒子方法得到的波函數(shù)與解析解吻合,即使粒子分布非均勻下得到SS-FPM結(jié)果仍與解析解一致.

為進(jìn)一步體現(xiàn)SS-FPM方法模擬Dirichlet邊值NLS方程的數(shù)值精度和收斂性,表1和表2分別給出了提出方法得到數(shù)值結(jié)果的誤差和收斂階.通過觀察表1,給出的粒子方法在粒子分布均勻和非均勻情況下得到的數(shù)值誤差差距不大,從而表明粒子方法模擬方程時無論區(qū)域是否規(guī)則都可以方便處理且具有較高的數(shù)值精度,較網(wǎng)格類方法具更好的靈活推廣應(yīng)用性.由表2可知,SS-FPM法具有二階收斂速度,與文獻(xiàn)[23]中SS-ICPSPH法和文獻(xiàn)[7]中網(wǎng)格類SS-FDM法的收斂階基本一致,但本文提出的粒子方法得到的數(shù)值誤差較SSICPSPH法和SS-FDM法的誤差小.SS-FPM方法相較于直接采用FPM方法,前者具有較小誤差和較快的收斂速度,這是因SS-FPM方法對非線性薛定諤方程中非線性項采用了較準(zhǔn)確的二階精度分裂格式.通過粒子法SS-FPM與網(wǎng)格類SS-FDM的構(gòu)造過程以及粒子分布非均勻情況下圖2和表1結(jié)果可知,本文粒子方法不受網(wǎng)格限制,可以在模擬區(qū)域任意布點情況下對問題進(jìn)行純無網(wǎng)格方法的模擬實現(xiàn),較網(wǎng)格類方法具有更好的靈活應(yīng)用推廣性,易被推廣應(yīng)用于復(fù)雜非規(guī)則區(qū)域上薛定諤問題的模擬.

表1 粒子分布均勻/非均勻兩種情況下的最大誤差erTable 1.Maximum error er under uniform/nonuniform particles distribution.

表2 四種不同方法在t=2時的數(shù)值收斂階Table 2.The rate of convergence obtained using four different methods at t=2.

4.1.2 周期邊值NLS方程

為體現(xiàn)提出的SS-FPM方法對帶周期性邊值NLS方程模擬的準(zhǔn)確性,考慮正方形區(qū)域Ω:[0,2π]×[0,2π]的周期邊界條件的非線性薛定諤方程,其對應(yīng)的方程和初值條件[11]為

初值條件

周期邊界條件

對應(yīng)的解析解為

模擬中,選取 k1=k2=1 ,粒子初始間距 π /64 ,時間步長 d t=10-4.圖3給出了兩個不同位置、不同時刻SS-FPM法得到的波函數(shù)實部變化曲線,并與解析解進(jìn)行比較.通過圖3可知,SS-FPM方法準(zhǔn)確地展示了波函數(shù)實部隨時間演化呈現(xiàn)出的周期性波形,且數(shù)值結(jié)果與解析解吻合,從而表明提出的粒子方法可以準(zhǔn)確地模擬帶周期邊值的NLS方程.表3給出了周期邊界條件下當(dāng)t=2時刻三種不同方法的數(shù)值收斂階.由表3可知,本文提出的SS-FPM法與其他兩種數(shù)值方法收斂階接近,且較SS-ICPSPH方法誤差更小; 給出的SS-FPM法模擬周期邊界下非線性薛定諤問題具有二階收斂速度.

圖3 兩個不同位置不同時刻ψ實部變化曲線圖 (a)沿對角線; (b)沿y=πFig.3.The change curve of real part at two positions with different times: (a)Along the diagonal; (b)along y=π .

表3 三種不同方法在t=2時的數(shù)值收斂階Table 3.The rate of convergence obtained using three different particle methods at t=2.

4.2 數(shù)值模擬預(yù)測

為進(jìn)一步體現(xiàn)SS-FPM方法模擬NLS/GP方程的能力,采用SS-FPM對帶周期邊界NLS方程中孤立波奇異特性和BEC中GP方程描述的孤立波量子化渦旋過程進(jìn)行數(shù)值預(yù)測.

4.2.1 具有奇異性周期邊值NLS方程

考慮正方形區(qū)域 Ω :[0,2π]×[0,2π]的周期邊界條件的非線性薛定諤方程,其對應(yīng)的方程和初值條件[11]為

初值條件

其中β=1.

該算例為無解析解的帶周期邊界NLS方程,它將描述隨時間演化孤立波出現(xiàn)奇異特性,常被用來驗證一種數(shù)值方法預(yù)測周期性NLS方程出現(xiàn)奇異值現(xiàn)象的可靠性和穩(wěn)定性[11].本小節(jié)運用SSFPM粒子法對該算例進(jìn)行了模擬,并與文獻(xiàn)[7]中SS-FDM法結(jié)果進(jìn)行對比(見圖4).觀察圖4知,SS-FPM方法得到的帶周期邊界NLS方程的奇異值與SS-FDM結(jié)果吻合,表明提出的粒子方法模擬預(yù)測NLS方程描述的奇異特性是可靠的.

圖4 兩個不同時刻波函數(shù) | ψ| 三維圖和等值線圖 (a1),(a2)t=0; (b1),(b2)t=0.0108Fig.4.The 3D graphs and contour of | ψ| at two different times: (a1),(a2)t=0; (b1),(b2)t=0.0108.

4.2.2 BEC中角動量旋轉(zhuǎn)GP方程

為驗證SS-FPM方法預(yù)測GP方程描述BEC動力學(xué)特性的有效性,考慮如下二分量帶旋轉(zhuǎn)項GP方程,并與其他數(shù)值結(jié)果進(jìn)行比較.正方形區(qū)域 Ω :[-8,8]×[-8,8]上具有二分量角動量旋轉(zhuǎn)GP方程,及其對應(yīng)的初邊值條件[8]為:

本文模擬中選取w1=w2=1.5x2+0.5y2,Θ=0.7,σ=-100,?=0.8.

圖5給出了兩種數(shù)值方法得到的不同時刻沿x軸變化的波函數(shù)曲線.由圖5中兩個不同時刻波函數(shù)變化曲線可知,BEC動力學(xué)特性在角動量旋轉(zhuǎn)下的變化是復(fù)雜的; SS-FPM方法得到的數(shù)值結(jié)果與SS-FDM法[7]得結(jié)果吻合,表明提出的粒子方法模擬預(yù)測GP方程描述BEC動力學(xué)性質(zhì)是可靠的.圖6給出了三個物理量 R e(ψ),Im(ψ),|ψ| 三維數(shù)值結(jié)果,可以明顯觀察到角動量旋轉(zhuǎn)項影響下隨時間演化的量子化渦旋變化情況.值得注意的是圖5和圖6僅展示了第一分量的數(shù)值結(jié)果,這是因為根據(jù)選取模擬參數(shù)兩個分量變化是一致的.

圖5 兩個不同時刻 | ψ| 沿x軸 (y=0)變化曲線 (a)t=0.05; (b)t=0.25Fig.5.The change curve of | ψ| along x-axis (y=0)at two different times: (a)t=0.05; (b)t=0.25.

圖6 兩個不同時刻 R e(ψ),Im(ψ),|ψ| 的三維數(shù)值結(jié)果 (a1),(a2),(a3)t=0; (b1),(b2),(b3)t=0.25Fig.6.Three-dimensional numerical results of R e(ψ),Im(ψ),|ψ| at two different times: (a1),(a2),(a3)t=0; (b1),(b2),(b3)t=0.25.

5 結(jié) 論

為提高直接推廣有限點集法模擬二維非線性薛定諤方程或Gross-Pitaevskii方程的數(shù)值精度,本文給出了基于分裂格式的有限點集法(SS-FPM),該方法兼顧時間分裂格式和傳統(tǒng)FPM方法的優(yōu)點.數(shù)值算例中,考慮了粒子分布均勻和非均勻情況下帶有不同邊界條件的二維非線性薛定諤方程,并與解析解進(jìn)行對比,對SS-FPM方法的精度和收斂性進(jìn)行了分析,驗證了數(shù)值預(yù)測的準(zhǔn)確性; 隨后運用SS-FPM法對無解析解NLS/GP方程進(jìn)行了數(shù)值預(yù)測,與其他數(shù)值結(jié)果進(jìn)行了比較.數(shù)值結(jié)果表明:

1)SS-FPM方法模擬二維非線性薛定諤方程具有二階收斂速度,較已有分裂有限差分法具有較好的靈活應(yīng)用性,且在粒子分布非均勻情況下仍具有較高數(shù)值精度;

2)SS-FPM法能夠成功地預(yù)測帶周期邊界條件下孤立波傳播過程的奇異現(xiàn)象,且與其他數(shù)值結(jié)果相吻合;

3)SS-FPM法準(zhǔn)確地預(yù)測了帶外旋轉(zhuǎn)BEC二分量的孤立波值量子化渦旋隨時間的演變過程.

目前未見文獻(xiàn)將時間分裂格式與FPM耦合對非線性薛定諤方程進(jìn)行模擬研究,本文針對不同邊界條件下二維非線性薛定諤方程給出的SSFPM法較網(wǎng)格類方法具有更好的靈活推廣應(yīng)用性,為復(fù)雜區(qū)域上含孤立波非線性問題的數(shù)值預(yù)測提供了一種準(zhǔn)確有效的粒子方法.