非線性薛定諤方程的幾種差分格式

孫傳志， 汪佳玲

(南京信息工程大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 江蘇 南京 210044)

1926年，物理學(xué)家薛定諤提出薛定諤方程，作為量子力學(xué)領(lǐng)域的基本方程，它對(duì)物理領(lǐng)域的研究具有深遠(yuǎn)意義.隨著社會(huì)的發(fā)展進(jìn)步，單純線性模型已經(jīng)不足以描述這個(gè)世界的所有現(xiàn)象，學(xué)者們將視線逐漸轉(zhuǎn)移到非線性模型的研究.20世紀(jì)70年代，在含非線性項(xiàng)的色散方程的研究中，Hasegawa推導(dǎo)出非線性薛定諤(nonlinear Schr?dinger equation，NLS)方程，廣泛應(yīng)用于應(yīng)用數(shù)學(xué)和量子力學(xué)等物理領(lǐng)域，如光脈沖在色散與非線性介質(zhì)中的傳輸，原子激光產(chǎn)生的Bose-Einstein凝聚效應(yīng)，電磁場中超導(dǎo)電子的運(yùn)動(dòng)等.

NLS方程的一般形式為iut+uxx+β|u|2u=0.其中，β為實(shí)常數(shù)；i為虛數(shù)單位.NLS方程對(duì)分析求解具有一定難度，沒有統(tǒng)一和高效的方法，在對(duì)非線性方程的各類研究中發(fā)現(xiàn)，非線性方程在很多情況下都不具備解析解，只有給出具體的限定條件才能得到確切的結(jié)果.關(guān)于同一個(gè)NLS方程，若給定的初始波函數(shù)不同，波函數(shù)的隨時(shí)演化模式也會(huì)產(chǎn)生變化，為此，利用數(shù)值模擬研究NLS方程.

隨著科學(xué)技術(shù)的創(chuàng)新，非線性薛定諤方程的研究越來越受到國內(nèi)外學(xué)者的重視.張?jiān)品錥1]運(yùn)用算子半群方法證明了NLS方程的解存在唯一性及解的一些性質(zhì)；孟佳[2]利用有限差分法對(duì)具體的NLS方程問題展開數(shù)值求解；Borhanifar等[3]研究薛定諤方程近似解的構(gòu)造，并在實(shí)驗(yàn)中加以驗(yàn)證；Ashyralyev等[4]研究含相關(guān)系數(shù)的高維分?jǐn)?shù)階薛定諤微分方程混合問題的一階及二階精度差分格式，并得到這些差分格式解的穩(wěn)定性分析；Zisowsky等[5]建立并分析求解NLS方程的不同有限差分格式的離散人工邊界條件，借助數(shù)值算例證明其穩(wěn)定性和正確性.

偏微分方程領(lǐng)域中的數(shù)值計(jì)算方法都很經(jīng)典，主要包括算子分裂法，有限元法，有限差分法等[6].有限差分法是求常(偏)微分方程及方程組的定解的方法，它較為靈活簡單，普遍性強(qiáng)，因而方便利用計(jì)算機(jī)進(jìn)行編程求解.在相關(guān)格式數(shù)值研究中，王海等[7]簡要分析二階中心差分的理論基礎(chǔ)，在此基礎(chǔ)上構(gòu)造聲波方程的該格式，建立相關(guān)模型并進(jìn)行數(shù)值模擬；吳宏偉[8]證明緊致差分格式的解的相關(guān)性質(zhì)，并利用該格式數(shù)值求解一類半線性拋物型方程；李華等[9]對(duì)Crank-Nicolson差分格式的誤差、穩(wěn)定性等展開分析證明，并聯(lián)立具體的數(shù)值算例加以論證；文獻(xiàn)[10-11]提出非線性差分格式，證明解存在唯一性和二階收斂L∞范數(shù)，并提出一種迭代算法求解非線性差分格式；Wang[12]構(gòu)造最優(yōu)收斂速度的緊致有限差分格式，即L∞范數(shù)下空間上4階及時(shí)間上2階精度；Patel等[13]針對(duì)具體偏微分方程提出一種無條件穩(wěn)定的緊致有限差分格式，并證明所提差分格式的穩(wěn)定性、一致性和收斂性.

NLS方程格式守恒性一直備受學(xué)者們的重視[14-17]，張魯明等[18-19]針對(duì)NLS方程構(gòu)造新的守恒差分格式，并對(duì)該守恒格式的收斂穩(wěn)定性展開證明；Lü等[20]推導(dǎo)出NLS方程的一種新求解方法，并利用實(shí)驗(yàn)證明該方法能準(zhǔn)確保持電荷守恒、能量守恒；Ismail等[21]提出求解NLS方程的一種線性隱式守恒方法，并在數(shù)值實(shí)驗(yàn)中表明該方法在時(shí)間和空間上均具備二階精度且能精確地保持能量守恒；Wang[22]提出一種求解非線性耦合薛定諤方程的能量守恒算法，并分析所提算法的可解性、穩(wěn)定性和誤差估計(jì)；Wang[23]提出并分析一個(gè)關(guān)于非線性耦合薛定諤方程的線性守恒差分格式；文獻(xiàn)[24-26]針對(duì)帶周期邊界條件的非線性薛定諤方程，分別提出緊差分的能量守恒格式，并給出無網(wǎng)格比的誤差估計(jì)；Gong等[27]利用Fourier擬譜方法構(gòu)造求解2維非線性薛定諤方程的一個(gè)能量和質(zhì)量守恒格式并加以分析；He[28]構(gòu)造并分析了非線性薛定諤方程的一個(gè)質(zhì)量且能量守恒的局部非連續(xù)Galerkin方法；Cui等[29]結(jié)合SAV方法構(gòu)造非線性薛定諤方程的任意高階的保結(jié)構(gòu)指數(shù)Runge-Kutta方法；Wang等[30]研究兩個(gè)高維薛定諤方程的傅里葉擬譜格式及方程的爆破行為.本文在一定初值、邊值條件下，結(jié)合不同的差分格式對(duì)NLS方程進(jìn)行數(shù)值求解.

1 數(shù)值格式的構(gòu)造

考慮如下NLS方程的初值、邊值問題，即

iut+uxx+β|u|2u=0，x∈R，t>0，

(1)

初值條件為

u(x，0)=u0(x)，x∈R，

(2)

邊值條件為

u(x，t)→0， |x|→∞，t>0.

(3)

式(1)～(3)中：u0(x)為已知光滑函數(shù)，函數(shù)隨著|x|→∞迅速向0衰減.

當(dāng)|x|≥1時(shí)，方程的解飛快地向0衰減，所以可在一個(gè)有限的區(qū)域Ω=(xL，xR)上對(duì)式(1)～(2)展開數(shù)值求解，其中，xL，xR≥1，得到如下初值、邊值問題，即

(4)

初值條件為

(5)

邊值條件為

u(x，t)=0，x=xL或x=xR，x∈?Ω.

(6)

式(4)～(6)的質(zhì)量守恒律為

(7)

能量守恒律為

(8)

2 幾種差分格式

2.1 向前Euler格式

(9)

由二階中心差分算子和向前差分算子，有

(10)

(11)

將式(10)，(11)代入式(9)，有

(12)

結(jié)合式(5)～(6)，有

(13)

(14)

(15)

(16)

(17)

以上格式即為NLS方程的向前Euler格式，是一個(gè)非線性顯式格式，記

(18)

向前Euler格式(15)～(17)的數(shù)值解在時(shí)間和空間方向分別具有1階、2階精度.

2.2 Crank-Nicolson格式

(19)

式(19)中：0≤n≤N-1；1≤j≤J-1.

結(jié)合式(10)及如下式(20)～(22)，即

(20)

(21)

(22)

并利用向前Euler格式類似的推導(dǎo)過程，有

(23)

(24)

(25)

(26)

與向前Euler格式不同，Crank-Nicolson格式是一個(gè)非線性隱式格式，Crank-Nicolson格式(24)～(26)的數(shù)值解在時(shí)間和空間方向上都具有2階精度.

2.3 緊致差分格式

(27)

式(27)中：0≤j≤J；0≤n≤N-1.

結(jié)合式(20)～(22)，有

(28)

式(28)中：0≤n≤N-1；0≤j≤J-1；c3是與h和τ無關(guān)的常數(shù).

對(duì)式(28)兩邊同時(shí)乘以Ah，經(jīng)過變換有

(29)

由于

(30)

(31)

將式(31)代入式(29)，有

(32)

(33)

(34)

(35)

緊致差分格式也是一個(gè)非線性隱式格式，利用前面定義的矩陣A和H，可以將式(33)表示為

(36)

緊致差分格式(33)～(35)的數(shù)值解在時(shí)間和空間方向上分別具有2階，4階精度.

3 守恒性

定理1Crank-Nicolson格式(24)～(26)精確保離散質(zhì)量守恒和能量守恒，質(zhì)量守恒和能量守恒分別為

(37)

(38)

式(37)～(38)中：n=0，1，2，…，N.

證明:式(24)的緊形式為

(39)

對(duì)上式取虛部，有

將(un+1-un)與式(39)兩邊同時(shí)作內(nèi)積，并取實(shí)部，有

(40)

定理2緊致差分格式(34)～(36)精確保離散質(zhì)量守恒和能量守恒，質(zhì)量守恒和能量守恒分別為

(41)

(42)

將式(36)的兩邊與(un+1-un)作內(nèi)積并取實(shí)部，得

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

4.1 孤立波的演化實(shí)驗(yàn)

圖1 初始條件下t=0 s時(shí)的孤立波波形

(a) t=1 s (b) t=10 s

Crank-Nicolson格式和緊致差分格式的孤立波波形分別如圖3，4所示.

(a) t=1 s (b) t=10 s

(a) t=1 s (b) t=10 s

由圖2～4可知：當(dāng)t=10 s時(shí)，向前Euler格式所得的孤立波波形有一些振蕩，可能是因?yàn)橄蚯癊uler格式作為顯格式，對(duì)差分步長較敏感.因此，相對(duì)其他兩種格式而言，向前Euler格式稍顯不足.

當(dāng)t=10 s時(shí)，選取時(shí)間步長(τ=0.000 01)繼續(xù)采用向前Euler格式進(jìn)行數(shù)值模擬，孤立波波形，如圖5所示.當(dāng)t=10s時(shí)，孤立波波形隨著時(shí)間步長τ的減小，振蕩越弱即波形更加穩(wěn)定.

圖5 向前Euler格式在τ=0.000 01的孤立波波形

4.2 守恒量的保持實(shí)驗(yàn)

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

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

由圖6，7可知：Crank-Nicolson格式和緊致差分格式均能很好地保持離散質(zhì)量守恒和能量守恒.

4.3 數(shù)值格式的計(jì)算效率實(shí)驗(yàn)

5 結(jié)論

分別利用經(jīng)典的向前差分算子、二階中心差分算子、Crank-Nicolson方法和緊致差分算子構(gòu)造了向前Euler格式、Crank-Nicolson格式和緊致差分格式.引用守恒性相關(guān)理論知識(shí)證明了Crank-Nicolson格式和緊致差分格式精確保持離散質(zhì)量守恒和能量守恒.

通過數(shù)學(xué)軟件MATLAB，分別對(duì)3種格式進(jìn)行了孤立波的數(shù)值實(shí)驗(yàn)和計(jì)算，并研究了Crank-Nicolson格式和緊致差分格式對(duì)守恒量的保持.雖然3種格式均能對(duì)NLS方程進(jìn)行數(shù)值求解，但是具有不同的特性.向前Euler差分格式作為非線性顯式格式，易于理解和編程，并且計(jì)算時(shí)間較少.而Crank-Nicolson格式和緊致差分格式都是非線性隱式格式，通過數(shù)值實(shí)驗(yàn)可以驗(yàn)證它們都能夠精確保持離散質(zhì)量守恒和能量守恒，與理論證明相符合.

這些格式不僅能應(yīng)用于NLS方程，還能應(yīng)用于偏微分方程領(lǐng)域的諸多方程中，如KdV方程、Klein-Gordon方程等.每種差分格式都具有不同的優(yōu)劣性，在具體的實(shí)際情況下，可采用較為理想的差分格式進(jìn)行數(shù)值計(jì)算.