含阻尼效應(yīng)的非線性薛定諤方程的共形分裂高階緊致差分格式

羅奕楊,王 蘭*,萬 隆,孔令華

(1.江西師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,江西 南昌 330022;2.豫章師范學(xué)院小學(xué)教育學(xué)院,江西 南昌 330103)

0 引言

非線性薛定諤方程是在數(shù)學(xué)物理中非常重要的模型,它在材料科學(xué)、等離子體物理學(xué)、量子物理學(xué)、非線性光學(xué)、雙分子動(dòng)力學(xué)等領(lǐng)域中得到了廣泛應(yīng)用.含有阻尼效應(yīng)的非線性薛定諤方程描述了在非線性介質(zhì)中的幾種諧振現(xiàn)象,特別是在垂直振蕩水槽中的非線性法拉第諧振[1]和相敏放大器對光纖中孤子的影響[2].本文考慮如下帶有阻尼效應(yīng)的非線性薛定諤(DNLS)方程:

iψt+ψxx+β|ψ|2ψ+iαψ=0,x∈Ω,0

(1)

的周期初邊值問題,取初始條件為

ψ(x,0)=ψ0(x),x∈Ω,

(2)

其中ψ(x,t)是定義在Ω×[0,T]上的未知復(fù)函數(shù),i是虛數(shù)單位,β是不為0的常數(shù)(β>0和β<0分別表示方程是聚焦和散焦的情況),α(>0)表示系統(tǒng)耗散率.若α=0,則方程變?yōu)橐话愕难Χㄖ@方程,其質(zhì)量與能量以及動(dòng)量是守恒的.

根據(jù)R. McLachlan等[3]的研究,初邊值問題(1)～(2)保持共形質(zhì)量守恒律

以及共形動(dòng)量守恒律

近年來,科研工作者對DNLS方程進(jìn)行了一些研究,構(gòu)造了一些數(shù)值格式,如Bao Weizhu等[4]提出了2階時(shí)間分裂譜方法,Jiang Chaolong等[5-6]利用傅里葉偽譜法構(gòu)造了非線性隱式與線性隱式2種方法,Cai Jiaxiang等[7]提出了緊致隱式積分算子方法,Hu Weipeng等[8]提出了DNLS方程的保辛格式等.上述工作雖然有很多優(yōu)點(diǎn),但是部分方法依然存在精度不高、計(jì)算效率較低或共形質(zhì)量守恒律不保持等問題.

基于此,本文提出了一種高精度、高效率且保持共形質(zhì)量守恒律的有限差分方法,該方法結(jié)合了時(shí)間分裂方法[9-10]與高階緊致方法[11](這2種方法已經(jīng)得到廣泛運(yùn)用[12-15]).其中時(shí)間分裂方法可以將復(fù)雜的原問題分裂成幾個(gè)易于求解的子問題,便于構(gòu)造簡便的數(shù)值格式;且高階緊致方法是一種利用了更少的模板節(jié)點(diǎn)達(dá)到了更高精度的方法,它適合于構(gòu)造高精度、高效率的數(shù)值格式.

1 共形分裂高階緊致格式的構(gòu)造

本部分先對Ω×[0,T]進(jìn)行網(wǎng)格剖分.設(shè)Ω=[a,b],將區(qū)間[a,b]作J等分,將區(qū)間[0,T]作N等分,記空間步長與時(shí)間步長分別為h、τ,則h=(b-a)/J,τ=T/N,可以得到時(shí)空區(qū)域Ω×[0,T]的一致網(wǎng)格剖分

Ωhτ={(xj,tn)|xj=a+jh,tn=nτ,j=0,1,…,J,n=0,1,…,N}.

為構(gòu)造DNLS方程(1)的高階精度格式,首先回顧一下高階緊致方法.對于具有一定光滑性的函數(shù)u(x),其2階導(dǎo)數(shù)u″(x)可以由u(xj+1)、u(xj)、u(xj-1) 3點(diǎn)近似[16],即

(u(xj+1)-2u(xj)+u(xj-1))/h2=u″(xj)+h2u(4)(xj)/12+O(h4)=u″(xj)+O(h2),

若略去高階項(xiàng)O(h2),則得到的格式具有2階精度.若再將上式中的4階導(dǎo)數(shù)d4u(x)/dx4寫成d2(d2u(x)/

dx2)/dx2,對前一2階微分采用2階中心差商離散,后一2階微分保持不變,則得到如下離散格式:

(u″(xj+1)+10u″(xj)+u″(xj-1))/12=(u(xj+1)-2u(xj)+u(xj-1))/h2+O(h4).

忽略高階項(xiàng)O(h4),得到截?cái)嗾`差為4階的緊致格式:

(3)

此格式只用了3個(gè)節(jié)點(diǎn)xj-1、xj、xj+1就得到了4階精度的差分格式,若采用一般的差分離散,則這是不可能達(dá)到的.在周期邊界條件下,將其寫成矩陣形式u″=P-1Qu,其中

為構(gòu)造DNLS方程(1)高效的分裂步高階緊致格式,先把它分裂成如下形式:

iψt+ψxx=0,

(4)

iψt+β|ψ|2ψ=0,

(5)

iψt+iαψ=0,

(6)

由用Crank-Nicolson格式對式(4)進(jìn)行時(shí)間離散和對高階緊致格式(3)進(jìn)行空間離散可得

(7)

這說明各個(gè)點(diǎn)在各個(gè)時(shí)刻處的質(zhì)量|ψ(x,t)|2保持不變,即

|ψ(x,tn+1)|2=|ψ(x,tn)|2.

這說明方程(5)滿足逐點(diǎn)的質(zhì)量守恒律.這使得非線性方程(5)可以精確求解.將其代入式(5),則可得方程(5)的精確解

ψ(xj,tn+1)=eiβ|ψ(xj,tn)| 2τψ(xj,tn).

(8)

對于式(6),由觀察可知,此方程為線性常系數(shù)常微分方程,可以用分離變量法得其精確解

ψ(xj,tn+1)=e-ατψ(xj,tn).

(9)

最后,將方程(4)～(6)的解(7)～(9)利用Strang組合,得到在時(shí)間上2階精度的格式:

(10)

其中j=1,2,…,J,n=0,1,2,…,N-1.

格式(10)在時(shí)間上具有2階精度、在空間上具有4階精度;且在計(jì)算過程中,格式(10)在每次時(shí)間的推進(jìn)時(shí)只需要求解線性方程組以及簡單的代數(shù)運(yùn)算,不需要額外的非線性迭代,計(jì)算成本較低,計(jì)算效率非常高,這克服了已有數(shù)值格式的各種不足.不僅如此,格式(10)還保持DNLS方程固有的一些物理性態(tài),如共形質(zhì)量守恒.

2 數(shù)值格式的共形守恒律以及穩(wěn)定性

為分析格式(10)的共形守恒律,需要運(yùn)用到循環(huán)矩陣的一些性質(zhì).

引理1[15]設(shè)X、Y為同階循環(huán)矩陣,則其滿足如下性質(zhì):

(i)XY=YX,X-1、Y-1、XY也都是循環(huán)矩陣;

(ii)對于循環(huán)對稱正定矩陣X,存在唯一的矩陣Z,使得X=ZTZ.

定理1DNLS方程(1)的Strang分裂格式(10)保持如下的離散共形質(zhì)量守恒律,即

證將式(10)中的第1式與第2式兩邊取模并平方,關(guān)于空間指標(biāo)j求和,可得

‖ψ(1)‖2=e-ατ‖ψn‖2,‖ψ(2)‖2=‖ψ(1)‖2,

同理可得‖ψ(4)‖2=‖ψ(3)‖2,‖ψn+1‖2=e-ατ·‖ψ(4)‖2.

用ψ(3)+ψ(4)與式(10)中第3式作內(nèi)積,取虛部并利用引理1可得

‖ψ(3)‖2=‖ψ(2)‖2.

綜上所述有‖ψn+1‖2=e-2ατ‖ψn‖2,即格式(10)滿足離散共形質(zhì)量守恒律Mn+1=e-2ατMn.同時(shí)這表明了格式(10)是無條件穩(wěn)定的,原因是

‖ψn+1‖2=e-2ατ‖ψn‖2=e-2αtn+1‖ψ0‖2≤‖ψ0‖2.

3 數(shù)值實(shí)驗(yàn)

本部分主要通過數(shù)值實(shí)驗(yàn)來檢驗(yàn)所構(gòu)造的格式(10)的有效性,并驗(yàn)證該格式的收斂精度和對共形質(zhì)量守恒律的保持情況.此外,考慮到DNLS方程(1)具有共形動(dòng)量守恒律,定義方程(1)對應(yīng)的離散動(dòng)量表達(dá)式為

其中Dx為1階導(dǎo)數(shù)ψx的離散算子δ2x所對應(yīng)的反對稱矩陣.通過上述公式來模擬格式的動(dòng)量變化情況.

取空間為Ω=[-25,25],DNLS方程(1)中的參數(shù)β=2、α=0.005,取初值為ψ0(x)=sech(x)·e2ix.由于無法給出其精確解,所以為便于分析,將利用格式(10)在相對小的時(shí)間步長與空間步長下得到的數(shù)值解作為“精確解”.這里取h=1/32,τ=1×10-4.為了計(jì)算收斂階,使用公式

oτ=ln(‖e(h,τ1)‖γ/‖e(h,τ2)‖γ)/ln(τ1/τ2),

oh=ln(‖e(h1,τ)‖γ/‖e(h2,τ)‖γ)/ln(h1/h2)

來分別計(jì)算時(shí)間方向與空間方向的收斂階數(shù),其中‖e(h,τ)‖γ表示在空間步長h以及時(shí)間步長τ下誤差的γ范數(shù),γ=2或∞.

在計(jì)算空間方向的收斂階時(shí),可以選取不同的空間步長,并固定一個(gè)相對小的時(shí)間步長來數(shù)值模擬,以此來減少時(shí)間方向?qū)φ`差的影響;在計(jì)算時(shí)間方向的收斂階時(shí)同理.選擇τ=1×10-4,h=1/2、1/4、1/8、1/16來計(jì)算空間方向的收斂階,選擇h=1/32,τ=1/25、1/50、1/100、1/200來計(jì)算時(shí)間方向的收斂階,問題計(jì)算到t=1.從表1和表2容易看出:格式(10)在空間方向上有4階收斂率,在時(shí)間方向上有2階收斂率.這與理論分析結(jié)果相符.

表1 空間收斂階及誤差(τ=1×10-4)

表2 時(shí)間收斂階及誤差(h=1/32)

(a)α=0 (b)α=0.005圖1 在不同時(shí)間上的波形

圖2與圖3分別展示了當(dāng)α=0和α=0.005時(shí)質(zhì)量與動(dòng)量隨時(shí)間的演化情況.圖4描繪了質(zhì)量和動(dòng)量的衰減率隨時(shí)間的演化關(guān)系,其縱坐標(biāo)的數(shù)據(jù)由(ln(Qn/Q0))/2±0.1(其中Qn=Mn或In)計(jì)算得到,中間的線條是參照線y=-αt=-0.005t.從圖2可以看出:當(dāng)α=0時(shí),質(zhì)量與初始值相比是基本保持不變的,符合薛定諤方程的質(zhì)量守恒律,而動(dòng)量與初始值相比具有一定的誤差.由圖3和圖4可以看到:質(zhì)量隨時(shí)間的衰減率完全與理論分析的結(jié)果一樣,且動(dòng)量也隨著時(shí)間而衰減.

(a)質(zhì)量 (b)動(dòng)量圖2 當(dāng)α=0時(shí)質(zhì)量與動(dòng)量隨時(shí)間的演化

(a)質(zhì)量 (b)動(dòng)量圖3 當(dāng)α=0.005時(shí)質(zhì)量與動(dòng)量隨時(shí)間的演化

圖4 質(zhì)量與動(dòng)量隨時(shí)間的衰減率

4 結(jié)論

本文通過采用分裂步方法和高階緊致方法,構(gòu)造了含有阻尼效應(yīng)的非線性薛定諤方程的共形分裂高階緊致方法.通過采用分裂步技術(shù)把原方程分解成更簡單的2個(gè)微分方程,其中非線性問題能夠精確求解,這可以避免求解非線性代數(shù)方程組.對于空間導(dǎo)數(shù)的離散采用高效的高階緊致格式,這樣得到的數(shù)值格式不僅能夠保持原方程的共形結(jié)構(gòu),而且計(jì)算效率還非常高.這種數(shù)值格式的構(gòu)造思路可以推廣到其他微分方程保結(jié)構(gòu)格式的構(gòu)造中,在未來的工作中將深入討論此類格式.