姚夢(mèng)麗, 田朝薇, 翁智峰

(華僑大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 福建 泉州 362021)

Schr?dinger方程是量子力學(xué)領(lǐng)域的基本方程,對(duì)物理領(lǐng)域的研究具有深遠(yuǎn)的意義.非線性Schr?dinger方程(NLS方程)在量子力學(xué)、非線性光學(xué)、電磁學(xué)、等離子理論、固體物理等領(lǐng)域中有著廣泛的應(yīng)用.考慮以下NLS方程的初邊值問題,即

(1)

式(1)中：i為虛數(shù)單位;a,b,ρ,β均為常數(shù);α(x),v(x)均為已知函數(shù);u為未知函數(shù).

眾所周知,NLS方程滿足質(zhì)量守恒和能量守恒,有

(2)

(3)

式(2),(3)中：M(t)為質(zhì)量函數(shù);E(t)為能量函數(shù).

Schr?dinger方程作為一類重要的量子力學(xué)方程,在實(shí)際復(fù)雜的系統(tǒng)中,由于包含復(fù)數(shù)且為耦合的問題,其精確解往往不容易求得.因此,對(duì)其進(jìn)行數(shù)值求解具有重要意義.許多研究者設(shè)計(jì)了有效的數(shù)值算法求解非線性Schr?dinger方程.孫傳志等[1]研究NLS方程的幾種差分格式.王廷春等[2]用緊致有限差分格式求解非線性耗散Schr?dinger方程.Gong等[3]采用Fourier擬譜方法構(gòu)造二維NLS方程的能量和質(zhì)量守恒格式.Cui等[4]用標(biāo)量輔助變量(SAV)方法構(gòu)造NLS方程的任意高階保結(jié)構(gòu)指數(shù)Runge-Kutta方法.Feng 等[5]研究NLS方程的高階守恒SAV-Gauss配置有限元格式.張法勇等[6]研究帶五次項(xiàng)的NLS方程的守恒差分格式.Wang等[7]用兩層網(wǎng)格有限元方法對(duì)NLS方程進(jìn)行超收斂性分析.Wang等[8]采用伽遼金有限元法求解阻尼NLS方程.Hu等[9]采用牛頓迭代Crank-Nicolson有限元方法求解NLS方程.Chen等[10]采用兩層網(wǎng)格有限體積元法求解NLS方程.以上方法都是基于網(wǎng)格剖分思想求解NLS方程數(shù)值解問題.最近,一種新型的無網(wǎng)格方法重心插值配點(diǎn)法廣泛應(yīng)用于求解微分方程初邊值問題.重心插值配點(diǎn)法具有計(jì)算格式簡單、精度高、程序?qū)嵤┓奖恪⒐?jié)點(diǎn)適應(yīng)性好等特點(diǎn),特別使用Chebyshev節(jié)點(diǎn)的重心Lagrange插值公式可以有效克服Runge現(xiàn)象.目前,重心插值配點(diǎn)法已被推廣到求解高維Fredholm積分方程[11]、耦合Burgers方程[12]、粘彈性波方程[13]、Allen-Cahn方程[14-16]、Cahn-Hilliard方程[17]、Black-Scholes方程[18]、分?jǐn)?shù)階電報(bào)方程[19]等.

算子分裂法[20-22]將復(fù)雜的問題分解成幾個(gè)簡單的子問題,并根據(jù)子問題的性質(zhì)采用不同的求解方法,用算子分裂格式將子問題結(jié)合起來得到原問題的解,可降低原問題的求解規(guī)模,此方法廣泛應(yīng)用于復(fù)雜非線性模型的數(shù)值求解中.基于此,本文提出一種非線性Schr?dinger方程的算子分裂配點(diǎn)格式.

1 重心Lagrange插值配點(diǎn)格式

1.1 重心Lagrange插值

設(shè)有插值節(jié)點(diǎn)xj和一組對(duì)應(yīng)的實(shí)數(shù)gj(j=0,1,…,m),采用多項(xiàng)式插值,則在次數(shù)不超過m的多項(xiàng)式空間可以找到唯一的插值多項(xiàng)式p(x),滿足p(xj)=gj.將p(x)寫成Lagrange插值公式,即

(4)

式(4)中：Lj(x)為Lagrange插值基函數(shù),滿足

(5)

以及

(6)

令

l(x)=(x-x0)(x-x1)…(x-xm),

(7)

定義重心插值權(quán)wj為

(8)

則Lagrange插值基函數(shù)Lj(x)可表示為

(9)

將式(9)代入式(4),可得

(10)

當(dāng)p(x)=1時(shí),有恒等式

(11)

用式(10)除以式(11),可得重心Lagrange插值公式為

(12)

由式(8)可知：重心插值權(quán)僅依賴插值節(jié)點(diǎn)的分布,因此,對(duì)于一些特殊的節(jié)點(diǎn)分布可以簡化重心插值權(quán).特別地,當(dāng)選擇Chebyshev點(diǎn)族時(shí),可保證重心Lagrange插值具有良好的數(shù)值穩(wěn)定性.

利用重心插值逼近未知函數(shù)時(shí)均采用第二類Chebyshev節(jié)點(diǎn),即

(13)

其重心Lagrange插值的重心插值權(quán)為

(14)

1.2 重心Lagrange插值微分矩陣

對(duì)區(qū)間[a,b]上的節(jié)點(diǎn)a=x1

(15)

由式(15),未知函數(shù)g(x)在節(jié)點(diǎn)x0,x1,…,xm處的p階導(dǎo)數(shù)表示為

(16)

將式(16)寫成矩陣形式,即

g(p)=D(p)g.

(17)

由文獻(xiàn)[23]可知,一階和二階微分矩陣元素的計(jì)算公式為

(18)

通過數(shù)學(xué)歸納法,可得p階微分矩陣的遞推公式為

(19)

2 基于重心Lagrange插值的NLS方程算子分裂格式

2.1 Strang算子分裂法

將式(1)改寫為

iut=L(u)+N(u).

(20)

式(20)中：L(u)=-ρΔu;N(u)=-v(x)u-β|u|2u.

算子分裂法的本質(zhì)是將方程(20)分為兩個(gè)子問題,即

線性部分SA：iut=L(u)=-ρΔu,

(21)

非線性部分SB：iut=N(u)=-v(x)u-β|u|2u.

(22)

式(21),(22)中：SA和SB分別為方程(21)和方程(22)的精確解,可得Strang算子分裂格式為

(23)

2.2 SA和SB的數(shù)值逼近

設(shè)空間節(jié)點(diǎn)xj(j=0,1,…,m)為Chebyshev節(jié)點(diǎn),時(shí)間節(jié)點(diǎn)tk=kτ(k=0,1,2,…,n).

首先,在空間方向上采用重心Lagrange插值配點(diǎn)法,可得半離散格式為

(24)

將式(24)寫成矩陣形式,即

i(uh(t))t=-ρD(2)uh(t),

(25)

式(25)中：uh(t)=[u0(t),u1(t),…,um(t)]′.

(26)

因此,有

(27)

方程(22)采用解析求解,即

(28)

則算子分裂格式為

(29)

3 相容性分析

考慮線性子問題的空間半離散格式(24)的相容性分析.設(shè)p(x)為u(x)的Lagrange插值函數(shù),滿足p(xj)=u(xj),j=0,1,…,m.定義誤差函數(shù)r(x)為

r(x)=u(x)-p(x).

(30)

引理1[19]假設(shè)u∈Cm+1(a,b),則有

式中：lx表示區(qū)間[a,b]的一半.類似地,可得

為了給出線性子問題SA的空間半離散格式(24)的相容性分析,首先,定義算子Γ為

(31)

定理1如果u∈C1([0,T])∪Cm+1([a,b]),有

證明：由式(21)和式(24),可得

Γu(x,t)-Γu(xj,t)=iut(x,t)+ρuxx(x,t)-iut(xj,t)-ρuxx(xj,t)=
iut(x,t)-iut(xj,t)+ρuxx(x,t)-ρuxx(xj,t)=R1+R2.

(32)

式(32)中：

R1=iut(x,t)-iut(xj,t),R2=ρuxx(x,t)-ρuxx(xj,t).

(33)

對(duì)于R1,有

R1=iut(x,t)-iut(xj,t)=i(ut(x,t)-ut(xj,t))=irt(xj,t).

(34)

由引理1可得

(35)

類似地,可得

(36)

將式(35),(36)代入式(32)即可,證畢.

4 數(shù)值算例

定義L∞誤差和L2誤差分別為

式中：ui為點(diǎn)xi處的精確解;Ui為點(diǎn)xi處的數(shù)值解.

4.1 算例1

比較重心Lagrange插值配點(diǎn)格式和有限差分格式的數(shù)值結(jié)果,以驗(yàn)證插值配點(diǎn)格式的有效性及高精度.ρ=1,β=1,選取帶右端項(xiàng)的真解u(x,t)=costsin πx,0≤x≤1,右端項(xiàng)為f(x)=-isintsin πx-π2costsin πx+|costsin πx|2costsin πx,初始條件為u(x,0)=sin πx,0≤x≤1, 邊界條件為u(0,t)=u(1,t)=0,0≤t≤1,勢函數(shù)v(x)=0.

固定時(shí)間步長(τ=10-4),改變空間節(jié)點(diǎn)數(shù)(M),重心Lagrange插值配點(diǎn)格式和有限差分格式的誤差,如表1所示.由表1可知：算例1的重心Lagrange插值配點(diǎn)格式選取6個(gè)節(jié)點(diǎn)差值即可達(dá)到10-5量級(jí),而有限差分格式選取256個(gè)節(jié)點(diǎn)差值才能達(dá)到10-5量級(jí),這表明重心Lagrange插值配點(diǎn)格式在空間上用較少的點(diǎn)可達(dá)到更高的精度.

表1 重心Lagrange插值配點(diǎn)格式和有限差分格式的誤差(算例1)Tab.1 Errors of barycentric Lagrange interpolation collocation scheme and finite difference scheme (example 1)

固定空間節(jié)點(diǎn)數(shù)(M=8),改變時(shí)間步長,重心Lagrange插值配點(diǎn)格式的誤差及時(shí)間收斂階,如表2所示.表2中：Rt,1,Rt,2分別為L∞誤差和L2誤差的時(shí)間收斂階.由表2可知：基于重心Lagrange插值配點(diǎn)格式的NLS方程在時(shí)間上是二階精度.

表2 重心Lagrange插值配點(diǎn)格式的誤差及時(shí)間收斂階(算例1)Tab.2 Errors and time convergence orders of barycentric Lagrange interpolation collocation scheme (example 1)

不同數(shù)值格式的空間收斂階,如圖1所示.圖1中：Rs為空間收斂階.由圖1可知：有限差分格式的空間收斂階是二階精度;重心Lagrange插值配點(diǎn)格式的空間收斂階滿足指數(shù)收斂的性質(zhì).

圖1 不同數(shù)值格式的空間收斂階(算例1)Fig.1 Spatial convergence orders for different numerical schemes (example 1)

4.2 算例2

算例2的重心Lagrange插值配點(diǎn)格式和有限差分格式的誤差(τ=10-4),如表3所示.由表3可知：算例2的重心Lagrange插值配點(diǎn)格式選取8個(gè)節(jié)點(diǎn)差值即可達(dá)到10-4量級(jí),而有限差分格式選取128個(gè)節(jié)點(diǎn)差值也只能達(dá)到10-4量級(jí),故重心Lagrange插值配點(diǎn)格式具有較高的精度.

表3 重心Lagrange插值配點(diǎn)格式和有限差分格式的誤差(算例2)Tab.3 Errors of barycentric Lagrange interpolation collocation scheme and finite difference scheme (example 2)

重心Lagrange插值配點(diǎn)格式的數(shù)值解和精確解(M=64,t=1),如圖2所示.不同數(shù)值格式的空間收斂階,如圖3所示.重心Lagrange插值配點(diǎn)格式守恒量的保持情況,如圖4所示.

(a) 數(shù)值解 (b) 精確解圖2 重心Lagrange插值配點(diǎn)格式的數(shù)值解和精確解Fig.2 Numerical solution and exact solution of barycentric Lagrange interpolation collocation scheme

圖3 不同數(shù)值格式的空間收斂階(算例2)Fig.3 Spatial convergence orders for different numerical schemes (example 2)

(a) 總質(zhì)量 (b) 質(zhì)量誤差

由圖2可知：NLS方程的重心Lagrange插值配點(diǎn)格式是穩(wěn)定的.由圖3可知：有限差分格式是二階精度,重心Lagrange插值配點(diǎn)格式呈指數(shù)遞減.由圖4可知：基于重心Lagrange插值配點(diǎn)格式的NLS方程滿足質(zhì)量守恒和能量守恒.

4.3 算例3

通過算例3進(jìn)行孤立波的演化實(shí)驗(yàn).選取空間節(jié)點(diǎn)M=248,時(shí)間步長τ=10-4,對(duì)應(yīng)的計(jì)算區(qū)間為[-10π,10π]×[0,T],ρ=1,β=2,選取真解u(x,t)=exp(i(2x-3t))sech (x-4t),初始條件為u(x,0)=exp(2ix)sechx,邊界條件為u(a,t)=u(b,t)=0,勢函數(shù)v(x)=0.

算例3的重心Lagrange插值配點(diǎn)格式的誤差及時(shí)間收斂階(M=248),如表4所示.由表4可知：算例3驗(yàn)證了時(shí)間收斂階為二階精度.

表4 重心Lagrange插值配點(diǎn)格式的誤差及時(shí)間收斂階(算例3)Tab.4 Errors and time convergence orders of barycentric Lagrange interpolation collocation scheme (example 3)

初始孤立波波形,如圖5所示.

圖5 初始孤立波波形Fig.5 Initial isolated wave shape

重心Lagrange插值配點(diǎn)格式的孤立波波形,如圖6所示.由圖6可知：算例3的重心Lagrange插值配點(diǎn)格式的孤立波波形隨時(shí)間的變化而變化,符合物理意義,而且每個(gè)時(shí)刻的波形都沒有出現(xiàn)震蕩,這說明該格式是穩(wěn)定的.

(a) t=1 (b) t=3 (c) t=5圖6 重心Lagrange插值配點(diǎn)格式的孤立波波形Fig.6 Isolated wave shapes of barycentric Lagrange interpolation collocation scheme

重心Lagrange插值配點(diǎn)格式守恒量的保持情況,如圖7所示.由圖7可知：算例3的重心Lagrange插值配點(diǎn)格式保持質(zhì)量守恒和能量守恒.

(a) 總質(zhì)量 (b) 總能量圖7 重心Lagrange插值配點(diǎn)格式守恒量的保持情況(算例3)Fig.7 Conservation of conserved quantities of barycentric Lagrange interpolation collocation scheme (example 3)

5 結(jié)束語

將算子分裂格式與重心Lagrange插值配點(diǎn)格式相結(jié)合,提出一種求解NLS方程的有效數(shù)值格式.此外,對(duì)方程的非線性部分采用解析求解,降低非線性迭代的計(jì)算量,提高計(jì)算效率,并給出線性子問題空間半離散格式的相容性分析.最后,通過數(shù)值算例驗(yàn)證格式的高精度和有效性.與經(jīng)典的有限差分格式進(jìn)行比較可知,文中格式可用較少的點(diǎn)達(dá)到較高的精度.今后,可將該方法推廣到其他非線性微分方程,為解決同類問題提供一種高精度數(shù)值求解方法,如耦合的非線性Schr?dinger方程、時(shí)間分?jǐn)?shù)階非線性Schr?dinger方程等.