對(duì)數(shù)非線性薛定諤方程基態(tài)解的數(shù)值解法

馮子旭， 何維清， 張世全

(四川大學(xué)數(shù)學(xué)學(xué)院， 成都， 610064)

1 引 言

玻色-愛(ài)因斯坦凝聚態(tài)(BEC)[1]是由微觀粒子的量子特性引起的宏觀現(xiàn)象[2-5].在最近的物理實(shí)驗(yàn)中，物理學(xué)家把不同組分的BEC混合在一起，使這些BEC 結(jié)合在一起形成量子液滴(液滴態(tài)).液滴態(tài)[6]可以用波函數(shù)描述，滿足如下的對(duì)數(shù)Gross-Pitaevskii方程(LogGPE)：

(1)

這里x是空間變量，t是時(shí)間變量，Ω>0是角速度，β,σ是刻畫粒子間相互作用類型的參數(shù)，角動(dòng)量算子Lz和勢(shì)函數(shù)V(x)定義如下：

Lz=-i(x?y-y?x),

這里的d代表空間維度.

對(duì)于該對(duì)數(shù)非線性薛定諤方程，在實(shí)際應(yīng)用中我們主要關(guān)心其基態(tài)解和動(dòng)力學(xué)演化，本文的研究目標(biāo)就是構(gòu)造該方程基態(tài)解的數(shù)值解法.另一方面，對(duì)于經(jīng)典的非線性薛定諤方程基態(tài)解的研究已有豐富結(jié)果，如,對(duì)數(shù)非線性薛定諤方程(LogSE)基態(tài)解的存在性可以參考文獻(xiàn)[7].

在非線性薛定諤方程基態(tài)解的計(jì)算方面，Edwards和Burnett[8]用Runge-Kutta方法求解了一維和三維球?qū)ΨQ基態(tài)解；Chiofalo 等[9]類比動(dòng)力學(xué)計(jì)算方法提出了虛擬時(shí)間演化法；Bao和Tang[10]提出利用有限元法直接極小化帶質(zhì)量約束的能量泛函，進(jìn)而得到基態(tài)解；Bao[11-13]等結(jié)合梯度流和單位球投影方法提出了歸一化梯度流方法，并構(gòu)造了向后歐拉差分和向后歐拉譜方法，在適當(dāng)初值下有效計(jì)算了基態(tài)和激發(fā)態(tài).值得注意的是，這些方法都是針對(duì)非對(duì)數(shù)非線性薛定諤方程的.眾所周知，由于對(duì)數(shù)非線性項(xiàng)的奇異性，這些方法并不能直接應(yīng)用到對(duì)數(shù)非線性薛定諤方程基態(tài)解的計(jì)算.本文將基于正則化能量泛函和歸一化梯度流構(gòu)造一種求解對(duì)數(shù)非線性薛定諤方程基態(tài)的數(shù)值方法.

本文剩余部分安排如下.第2節(jié)介紹方程(1)的基態(tài)定義和物理意義及帶正則化參數(shù)1<ε?0的正則化方法.第3節(jié)給出了具體的數(shù)值方法和計(jì)算過(guò)程.第4節(jié)我們用數(shù)值算例驗(yàn)證方法的可靠性.最后，第5節(jié)給出全文的結(jié)論.

2 預(yù)備知識(shí)

方程(1)有兩個(gè)重要的相關(guān)物理量，分別是能量和質(zhì)量，定義如下：

(2)

(3)

其中

(4)

數(shù)學(xué)上，基態(tài)解φg(x)是能量泛函在質(zhì)量約束條件下的極小值點(diǎn)，即

(5)

由于對(duì)數(shù)非線性項(xiàng)的奇異性，即當(dāng)ρ→0+時(shí),logρ→-∞，我們不能直接使用前文的方法.為了解決這個(gè)難題，我們參考正則化(ERLogSE)方法[14]，對(duì)能量泛函的對(duì)數(shù)項(xiàng)引入正則化參數(shù)ε進(jìn)行修正，修正后的能量泛函為

(6)

其中正則化函數(shù)

或

(7)

根據(jù)正則化能量泛函，我們有如下的正則化方程：

fε(ρ)=(Fε(ρ))′

(8)

(9)

對(duì)于上式所得的基態(tài)解能量與原問(wèn)題的基態(tài)解能量的誤差，我們有下面的定理.

(10)

證明 首先，?φ∈S∩H1，有

利用不等式0≤log(1+x)≤x，?x≥0可得

-Cε≤E(φ)-Eε(φ)≤0.

根據(jù)基態(tài)的定義可知，

綜合上述不等式即可得到定理2.1.

3 數(shù)值解法

為求解(5)式的能量泛函極小化問(wèn)題，我們首先對(duì)(1)式使用虛時(shí)演化：t→-it,再利用離散歸一化梯度流(DNGF) 方法.

選擇一個(gè)固定的時(shí)間步長(zhǎng)Δt>0.對(duì)n= 0,1,2, …和一個(gè)時(shí)間序列tn=nΔt，在每個(gè)時(shí)間步t∈[tn,tn+1)，DNGF 方法如下：

?據(jù)墓志，韓顯宗卒于太和二十三年，時(shí)年三十四歲。參見(jiàn)趙超《漢魏南北朝墓志匯編》，天津古籍出版社2008年版，第39頁(yè)。

fε(|φε|2)φε,x∈Rd

(11)

(12)

φε(x,0)=φ0(x),x∈Rd，

(13)

i) 對(duì)無(wú)約束條件的能量泛函使用速降法;

ii) 為了滿足約束條件，在每個(gè)時(shí)間步將解拉回到解空間S.

接下來(lái)我們用向后歐拉傅里葉譜(BFFP) 方法來(lái)離散(11)～(13)式.由勢(shì)函數(shù)的性質(zhì)，(11)～(13)式的解會(huì)在無(wú)窮遠(yuǎn)處呈指數(shù)衰減到零，所以在實(shí)際計(jì)算時(shí)我們可以把計(jì)算區(qū)域截?cái)嗟揭粋€(gè)合適的有界區(qū)域上，并假設(shè)解滿足零邊界條件.不失一般性，我們以二維的BFFP方法為例，對(duì)其它維數(shù)來(lái)說(shuō)也是類似的.

在二維的情況下，對(duì)(11)～(13)式的截?cái)鄦?wèn)題如下：

fε(|φε|2)φε,x∈D,t∈[tn,tn+1)

(14)

(15)

φε(x,t)=0,x∈Γ=?D,t∈[tn,tn+1)

(16)

φε(x,0)=φ0(x),x∈D且有

(17)

這里的計(jì)算區(qū)域D=[a,b]×[c,d]，其中|a|， |b|，|c|和|d|是足夠大的常數(shù).選擇空間網(wǎng)格，其中J和K是兩個(gè)正偶數(shù)，網(wǎng)格點(diǎn)如下：

定義下標(biāo)集

TJK=(j,k)|j=1,2,…,J-1,

k=1,2,…,K-1,

k=0,1,2,…,K.

引入記號(hào)

則求解(14)～(15)式的BFFP方法離散格式為[11，15]：

(18)

(19)

對(duì)Δ，?x和?y,我們可以使用傅里葉譜方法逼近，即

(20)

(21)

(22)

在每一個(gè)時(shí)間步，為了求非線性方程組(18)式的解φε,(1)，我們使用帶穩(wěn)定參數(shù)項(xiàng)的離散傅里葉變換方法，迭代求解下面的線性方程組,直至收斂：

iΩyk([?x]φε,(1),m)jk-iΩxj([?y]φε,(1),m)jk-

(23)

對(duì)(23)式兩邊使用離散傅里葉變換，有

(24)

(25)

的φε,(1),m+1即為當(dāng)前時(shí)間步的解φε,(1).

4 數(shù)值模擬

表1 對(duì)不同參數(shù)β，兩種正則化方法得出的基態(tài)解能量及能量之差Tab.1 The energy and energy difference of the ground state solution obtained by two regularization methods for different β

由上表可知，在數(shù)值上來(lái)說(shuō)兩種正則化方法的結(jié)果幾乎沒(méi)有區(qū)別.因此，在后面計(jì)算二維情況的基態(tài)解時(shí)我們統(tǒng)一使用第一種正則化方法.

從上圖1可知，左圖的正則化方法能量關(guān)于參數(shù)ε的收斂階高于一階，略低于三階;而右圖的正則化方法能量關(guān)于參數(shù)ε的收斂階高于一階，非常接近二階,與理論分析相符，并且數(shù)值試驗(yàn)的結(jié)果優(yōu)于理論結(jié)果.再由上表2可知，計(jì)算耗時(shí)基本不受正則化參數(shù)大小的影響.所以在求解基態(tài)時(shí)我們可以選擇盡可能小的正則化參數(shù).

表2 對(duì)不同的正則化參數(shù)ε，兩種正則化方法得出的基態(tài)解的運(yùn)行時(shí)間Tab.2 Running time of the ground state solution obtained by the two regularization methods for different ε

下面我們用第一種正則化方法來(lái)研究在正則化參數(shù)ε=10-14時(shí)的質(zhì)量M0和能量E0隨計(jì)算時(shí)間t的變化，結(jié)果如下.

根據(jù)圖2，我們可以看到該數(shù)值方法在此算例中是質(zhì)量守恒，能量穩(wěn)定的.

圖2 質(zhì)量M0隨時(shí)間的變化(a);能量E0隨時(shí)間的變化(b)Fig.2 The change of mass M0over time(a), and the change of energy E0 over time(b)

最后，應(yīng)用前面提出的數(shù)值方法計(jì)算LogSE的基態(tài)解.如文獻(xiàn)[6]，我們?cè)O(shè)置各參數(shù)為：M=1 000,γx=γy=γ=0.04,β=1,σ=1,D=[-30,30]×[-30,30].我們選擇如下函數(shù)的線性組合作為初始函數(shù)[16]：

φ(x)=C|x|ke-γ|x|2/2+ikθ,

θ(x,y)=arctan(y/x),

C是歸一化常數(shù)，使得‖φ0‖2=M.我們研究如下四種情況：

Case 1 Ω/γ=0,φ0(x,y)=Ce-γ|x|2/2;

Case 2 Ω/γ=0.375,

φ0(x,y)=C(x+iy)e-γ|x|2/2;

Case 3 Ω/γ=0.45,

φ0(x,y)=C|x|e-γ|x|2/2+iθ;

Case 4 Ω/γ=0.5125,

x1=(5,0)T,x2=(-5,0)T,

x3=(0,-5)T.

進(jìn)一步，我們?nèi)∵x擇非常小的正則化參數(shù)ε=10-14，時(shí)間步長(zhǎng)Δt=0.001以及網(wǎng)格J=K=28，計(jì)算結(jié)果如下.

圖3 (a～d)依次為Case 1～4的基態(tài)解Fig.3 From (a) to (d) are the ground states of Case 1～4

我們的計(jì)算結(jié)果與文獻(xiàn)[6]中的結(jié)果相符，這說(shuō)明了我們的方法的可靠性.

5 結(jié) 論

本文提出了帶正則化參數(shù)的ERLogSE方法，解決了對(duì)數(shù)非線性項(xiàng)在原點(diǎn)處的奇異性帶來(lái)的困難，并分析了正則化基態(tài)解與原問(wèn)題基態(tài)解的能量誤差估計(jì).基于正則化的能量泛函，本文設(shè)計(jì)了相應(yīng)的求基態(tài)解的數(shù)值格式，并用數(shù)值算例驗(yàn)證了方法的可靠性和理論分析的正確性.