求解Schr dinger方程的一種新能量守恒格式

摘要：針對一維非線性Schrdinger方程，通過使用多標(biāo)量輔助變量方法將原方程轉(zhuǎn)化為具有能量守恒性質(zhì)的等價系統(tǒng).對等價系統(tǒng)構(gòu)建帶有懲罰項的數(shù)值格式，該格式保持修正的能量守恒，并且在每一個時間步長只需求解帶有常系數(shù)矩陣的解耦系統(tǒng).最后通過數(shù)值實驗驗證了該數(shù)值格式的有效性.

關(guān)鍵詞：Schrdinger方程;多標(biāo)量輔助變量;能量守恒

中圖分類號：O241.82"" 文獻(xiàn)標(biāo)志碼：A

A New Energy Conservation Scheme for Solving Schrdinger Equations

WANG Ya-hui

（School of Science， Zhejiang Sci-Tech University， Hangzhou 310018， China）

Abstract：For the one-dimensional nonlinear Schrdinger equations， the multiple scalar auxiliary variable method is used to transform the original equation into an equivalent system with energy conservation property. A numerical scheme with penalty term is constructed for solving the equivalent system. The scheme preserves the modified energy conservation， and only needs to solve the decoupling system with constant coefficient matrix at each time step. Finally， the validity of the numerical scheme is verified by numerical experiments.

Key words：Schrdinger equation; multiple scalar auxiliary variable; energy conservation

0 引言

考慮如下一維非線性Schrdinger方程的初邊值問題：

iut+uxx+βu²u=0，

x∈0，L，t∈0，T，

ux，0=φx，x∈0，L，

u0，t=uL，t，t∈0，T，（1）

其中：β為給定實常數(shù);i為虛數(shù)單位;φ（x）為已知光滑復(fù)值函數(shù);u（x，t）為未知復(fù)值函數(shù).

孫志忠^[1]對非線性Schrdinger方程構(gòu)建一般的兩層非線性差分格式，該格式需要通過迭代算法求解，從而產(chǎn)生較高的計算代價.SHEN J等^[2]首次提出了標(biāo)量輔助變量（SAV）方法，該方法主要用于能量守恒或耗散的非線性系統(tǒng)數(shù)值求解.基于SAV方法構(gòu)建的數(shù)值格式在求解時，只需要在每一個時間步長中求解解耦的常系數(shù)線性方程組，這使得計算效率顯著提高，在偏微分方程的求解中有著廣泛的應(yīng)用.如JIANG C L等^[3]針對非線性Klein-Gordon方程，結(jié)合SAV方法的思想，構(gòu)建線性隱式指數(shù)積分保能格式，與一般的指數(shù)積分格式相比，所提出的新格式只需求解解耦的常系數(shù)線性系統(tǒng);ANTOINE X等^[4]對非線性Schrdinger/Gross-Pitaevskii方程構(gòu)建了線性隱式擬譜格式，該格式可以保持修正的能量，并且質(zhì)量在時間方向是近似三階的;CHENG Q等^[5]針對相場囊泡膜模型應(yīng)用多標(biāo)量輔助變量（MSAV）方法，構(gòu)造無條件能量穩(wěn)定和時間方向有二階精度的數(shù)值格式，更多內(nèi)容可參考文獻(xiàn)[6-9].在本文中，針對非線性Schrdinger方程使用MSAV方法，提出了可以

保持修正的能量守恒的數(shù)值格式，并且該數(shù)值格式只需要在每一個時間步長求解解耦的常系數(shù)線性方程組，與兩層非線性差分格式相比計算效率更高.

本文結(jié)構(gòu)如下，在第一節(jié)中，對Schrdinger方程使用MSAV方法后的等價方程建立了差分格式，并給出了該格式的能量守恒證明;在第二節(jié)中，

給出有效求解所得差分格式的具體算法流程；在第三節(jié)中，采用數(shù)值算例驗證理論結(jié)果的有效性；在第四節(jié)中，給出文章的總結(jié).

1 格式建立

為了保證問題（1）可以滿足能量守恒，引入兩個標(biāo)量輔助變量（MSAV）

q：=qt= u⁴，1，

p：=p（t）=||u||²-||u⁰||²，

其中 u，v=∫L0uvdx 且 u= u，u.

將問題（1）改寫成如下的等價系統(tǒng)

ut=iuxx+iβu²u u⁴，1q+i1εpu，

qt=12 u⁴，1

2ut，u²u+2u²u，ut=

2qtReu²u，ut，

pt=2Reu，ut， （2）

其中，ε是一個較小的參數(shù)，將其作為懲罰項的系數(shù)，以此近似質(zhì)量守恒.

定理1 系統(tǒng)（2）保持修正的能量守恒

E（t）：=12||ux||²-β4q²-14εp²≡E（0）.

證明 將系統(tǒng)（2）中的第一個方程與ut作內(nèi)積，可得

ut，ut=iuxx，ut+

iβu²u u⁴，1q，ut+i1εpu，ut，（3）

第二個方程與q（t）相乘，根據(jù)q（t）的定義得

12ddtq²=qtqt=2Reu²u，ut.（4）

同樣的，第三個方程與p（t）相乘，可得

12ddtp²=ptpt=2Repu，ut.（5）

對方程（3）兩邊同時取虛部，并結(jié)合方程（4）和（5）得

0=-12ddtux²+β4ddtq²+14εddtp2，

則定理得證.

在給出系統(tǒng)（2）的數(shù)值格式前，首先引入一些記號.令h=L/M，τ=T/N，xi=ih，0≤i≤M，tn=nτ，0≤n≤N，其中M和N均為正整數(shù).記Ωhτ={（xi，tn）|0≤i≤M，0≤n≤N}和Vh={v|v={vⁿi}是定義在Ωhτ上的網(wǎng)格函數(shù)}.對任意網(wǎng)格函數(shù)u，v∈Vh，記

δtu^n+1/2i=uⁿ⁺¹i-uⁿi/τ，

δ²xuⁿi=uⁿi+1-2uⁿi+uⁿi-1/h²，

u^n+1/2i=（uⁿ⁺¹i+uⁿi）/2，

u^*，n+1/2i=3uⁿi/2-u^n-1i/2，

lt;u，vgt;=h∑M-1i=0uivi，

v= h∑^M-1i=0vi²，vSymboleB@=max0≤i≤M-1vi，

其中vi為vi的共軛.

在空間方向使用中心差分近似，時間方向線性部分使用Crank-Nicolson格式，非線性部分使用Adams-Bashforth外推近似，從而得到如下離散系統(tǒng)：

uⁿ⁺¹-uⁿ=iτAu^n+1/2+

iβτu^*，n+1/22u^*，n+1/2 lt;u^*，n+1/24，1gt;q^n+1/2+

iτεp^n+1/2u^*，n+1/2，1≤n≤N-1，

qⁿ⁺¹-qⁿ=

2Relt;u^*，n+1/22u^*，n+1/2，uⁿ⁺¹-uⁿgt; lt;u^*，n+1/24，1gt;，

1≤n≤N-1，

pⁿ⁺¹-pⁿ=

2Relt;u^*，n+1/2，uⁿ⁺¹-uⁿgt;，1≤n≤N-1，（6）

其中，A為中心差分算子所對應(yīng)的矩陣.

注1 u¹ ，q¹和p¹可由u⁰計算得到，即：

u¹-u⁰=iτAu¹+u⁰2+

iβτu⁰²u⁰ lt;u⁰⁴，1gt;q¹+q⁰2+

iτεp¹+p⁰2u⁰，

q¹-q⁰=2Relt;u⁰²u⁰，u¹-u⁰gt; lt;u⁰⁴，1gt;，

p¹-p⁰=2Relt;u⁰，u¹-u⁰gt;.

定理2 格式（6）保持無條件能量穩(wěn)定

Eⁿ⁺¹：=-12Relt;Auⁿ⁺¹，uⁿ⁺¹gt;-

β4qⁿ⁺¹²-14εpⁿ⁺¹²≡Eⁿ.

證明 將格式（6）的第一個方程與uⁿ⁺¹-uⁿ做內(nèi)積，并且等式兩邊同時取虛部得

Imlt;uⁿ⁺¹-uⁿ，uⁿ⁺¹-uⁿgt;=

Relt;τAu^n+1/2，uⁿ⁺¹-uⁿgt;+

Relt;βτu^*，n+1/22u^*，n+1/2 lt;u^*，n+1/24，1gt;q^n+1/2，

uⁿ⁺¹-uⁿgt;+Relt;τεpⁿ⁺¹²u^*，n+12，uⁿ⁺¹-uⁿgt;.

類似定理1的證明過程，有

0=τ2Relt;Auⁿ⁺¹，uⁿ⁺¹gt;

-τ2Relt;Auⁿ，uⁿgt;+βτ4［qⁿ⁺¹²-qⁿ²］+

τ4ε［pⁿ⁺¹²-pⁿ²］

成立，則定理得證.

2 數(shù)值實現(xiàn)

由于格式（6）可以化簡成一個解耦的系統(tǒng)進(jìn)行求解，無需進(jìn)行迭代求解，大大減少了運算代價.在該部分，將給出格式（6）的具體求解過程. 令

uⁿ⁺¹=u∧ⁿ⁺¹+qⁿ⁺¹u∨ⁿ⁺¹+pⁿ⁺¹uⁿ⁺¹，（7）

將其帶入到格式（6）的第一個方程中，化簡可得到如下線性解耦格式：

I-12iτAu∧ⁿ⁺¹=I+12iτAuⁿ+

12iβτu^*，n+1/22u^*，n+1/2 lt;u^*，n+1/24，1gt;qⁿ+

iτ2εu^*，n+1/2pⁿ，

I-12iτAu∨ⁿ⁺¹=iβτ2u^*，n+1/22u^*，n+1/2 lt;u^*，n+1/24，1gt;，

I-12iτAuⁿ⁺¹=iτ2εu^*，n+1/2.（8）

根據(jù)格式（8），可求解出u∧ⁿ⁺¹，u∨ⁿ⁺¹，uⁿ⁺¹，然后將式（7）代入到式（6），得到關(guān)于qⁿ⁺¹和pⁿ⁺¹的耦合的線性代數(shù)方程組：

qⁿ⁺¹-qⁿ=

2Relt;u^*，n+1/22u^*，n+1/2，

u∧ⁿ⁺¹+qⁿ⁺¹u∨ⁿ⁺¹+pⁿ⁺¹uⁿ⁺¹-uⁿgt;/

lt;u^*，n+1/24，1gt;，

pⁿ⁺¹-pⁿ=

2Relt;u^*，n+1/2，u∧ⁿ⁺¹+qⁿ⁺¹u∨ⁿ⁺¹+

pⁿ⁺¹uⁿ⁺¹-uⁿgt;.（9）

再通過求解方程組（9），可得qⁿ⁺¹和pⁿ⁺¹.最后，通過式（7）可得uⁿ⁺¹.具體的算法流程如圖1所示.

3 數(shù)值實驗

在該部分，用數(shù)值算例驗證格式（6）的有效性，并將該格式與兩層非線性差分格式^[1]

iδtu^n+1/2i+δ²xu^n+1/2i+

β2uⁿi²+uⁿ⁺¹i²u^n+1/2i=0，

1≤i≤M，1≤n≤N-1，

u⁰i=φxi，1≤i≤M，

uⁿ0=uⁿM，1≤n≤N-1進(jìn)行數(shù)值比較.（10）

令{uⁿi（h，τ）}表示空間步長為h，時間步長為τ時的數(shù)值解，{u²ⁿi（h，τ/2）}表示空間步長為h，時間步長為τ/2時的數(shù)值解.定義lSymboleB@范數(shù)意義下的數(shù)值誤差與收斂階為

Ordh=log2ESymboleB@h，τESymboleB@h/2，τ，

其中

ESymboleB@h，τ=

max0≤i≤M，0≤n≤Nuⁿi（h，τ）-uⁿ2i（h/2，τ），

Ordτ=log2ESymboleB@h，τESymboleB@h，τ/2，

其中

ESymboleB@h，τ=

max0≤i≤M，0≤n≤Nuⁿi（h，τ）-u²ⁿi（h，τ/2）.

例1 考慮如下問題

iut+uxx+2u²u=0，

x∈－10，10，t∈0，1，

ux，0=sech（x）exp（2ix），x∈－10，10，

u－10，t=u10，t，t∈0，1.

利用離散格式（6）和兩層非線性差分格式（10）求解上述問題，相關(guān)結(jié)果展示在表1-4和圖2中. 從表1、表2和圖2中可以看出，格式（6）（取參數(shù)ε=0.1）在空間和時間方向均為二階精度，且保持修正的能量守恒. 從表3和表4中可以看出兩層非線性差分格式（10）在時間方向和空間方向上均是二階收斂，并且通過對比表2和表4的CPU運行時間，表明格式（6）的計算效率更高.

4 結(jié)語

通過使用多標(biāo)量輔助變量（MSAV）方法，將一維非線性Schrdinger方程轉(zhuǎn)化為具有能量守恒性質(zhì)的等價系統(tǒng)，結(jié)合空間方向的中心差分近似，時間方向的Crank-Nicolson格式和Adams-Bashforth外推近似，構(gòu)建帶有懲罰項的數(shù)值格式，該格式保持修正的能量守恒，并且給出了求解該系統(tǒng)的算法流程圖，在每一個時間步長只需求解帶有常系數(shù)矩陣的解耦系統(tǒng)，與一般的非線性差分格式相比，可以有效提高計算效率.

參考文獻(xiàn)：

[1] 孫志忠.偏微分方程數(shù)值解法[M].北京：科學(xué)出版社，2005.

[2] SHEN J，XU J，YANG J.The scalar auxiliary variable （SAV） approach for gradient flows[J].Journal of Computational Physics，2018，353：407-416.

[3] JIANG C L，WANG Y S，CAI W J.A linearly implicit energy-preserving exponential integrator for the nonlinear Klein-Gordon equation[J].Journal of Computational Physics，2020，419：109690.

[4] ANTOINE X，SHEN J，TANG Q L.Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrdinger/Gross-Pitaevskii equations[J].Journal of Computational Physics，2021，437：110328.

[5] CHENG Q，SHEN J.Multiple scalar auxiliary variable （MSAV） approach and its application to the phase-field vesicle membrane model[J].SIAM Journal on Scientific Computing，2018，40（6）：A3982-A4006.

[6] SHEN J，XU J，YANG J.A new class of efficient and robust energy stable schemes for gradient flows[J].SIAM Review，2019，61（3）：474-506.

[7] FU Y Y，HU D D，WANG Y S.High-order structure-preserving algorithms for the multi-dimensional fractional nonlinear Schrdinger equation based on the SAV approach[J].Mathematics and Computers in Simulation，2021，185：238-255.

[8] WANG M，HUANG Q M，WANG C.A second order accurate scalar auxiliary variable （SAV） numerical method for the square phase field crystal equation[J].Journal of Scientific Computing，2021，88（2）：33.

[9] LI D F，LI X X，SUN H W.Optimal error estimates of SAV Crank-Nicolson finite element method for the coupled nonlinear Schrdinger equation[J].Journal of Scientific Computing，2023，97（3）：71.

［責(zé)任編輯：趙慧霞］