極坐標(biāo)下二維非線性薛定諤方程的有限差分方法

楊程程，張榮培

（沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，遼寧 沈陽 110034）

非線性薛定諤方程是量子力學(xué)中最重要的方程之一，在等離子物理、非線性光學(xué)、激光晶體中的自聚焦、晶體中熱脈沖的傳播以及在極低溫度下的Bose-Einstein 凝聚體的動(dòng)力學(xué)等領(lǐng)域內(nèi)有著重要的應(yīng)用［1-4］。

近年來，許多學(xué)者在求解非線性薛定諤方程時(shí)應(yīng)用了許多數(shù)值方法，例如有限差分方法［5］、有限元法［6］、譜方法［7］和緊致積分因子法［8］等等。但這些方法均在直角坐標(biāo)系下求解，而在極坐標(biāo)下求解的非線性薛定諤方程的文章比較少［9］，本文考慮在圓形區(qū)域上求解極坐標(biāo)下的二維非線性薛定諤方程。

式中，u(x,y)為復(fù)函數(shù)；i2=-1 為虛數(shù)單位；Δu=uxx+uyy為拉普拉斯算子。

邊界條件：在r方向的左邊界為Neumann 邊界；在r方向的右邊界為齊次Dirichlet 邊界；在θ方向?yàn)橹芷谶吔纭?/p>

1 數(shù)值方法

在式（1）中，將直角坐標(biāo)系下的x與y轉(zhuǎn)換為極坐標(biāo)形式，即x=rcosθ，y=rsinθ。其中拉普拉斯算子可表示為Δu=urr+r-1ur+r-2uθθ。式（1）在極坐標(biāo)下可改寫為式中，0

圖1 r方向的節(jié)點(diǎn)劃分

根據(jù)泰勒展開式，式（3）中的導(dǎo)數(shù)可以離散為如下形式：

令u在網(wǎng)格節(jié)點(diǎn)(rj,θk)的數(shù)值解為uj,k，略去余項(xiàng)，得到式（3）的二階中心差分格式：

由于式（1）在r方向的左端是齊次Neuman 邊界，在r方向的右端是齊次Dirichlet 邊界，在θ方向是周期邊界，所以不存在邊界的奇性。

為了能夠用Kronecker積的形式來表示式（8），首先將式（8）寫成矩陣的形式，然后再將矩陣的形式轉(zhuǎn)化成Kronecker 積的形式。定義數(shù)值解uj,k的矩陣為

將非線性項(xiàng)寫成向量的形式記為F?？紤]到u在θ方向是以2π 為周期的，因此邊界值為uNr+1,k=0，uj,0=根據(jù)式（8）、式（9）和式（10），可將非線性常微分方程式（11）寫成如下緊致形式：

式中，矩陣的定義如下：

下面將式（12）向量化，將U和F按照列重新排成向量的形式，記U=vec(U)，F(xiàn)(U)=vec(F)，由Kronecker 積的性質(zhì)(BT?A)vec(X)=vec(AXB)［10］可得

采用積分因子方法求解式（13），將時(shí)間區(qū)域剖分為tn=nΔt。其中，Δt為時(shí)間步長(zhǎng)；n=1,2,…,N。令，式（13）寫成

式（14）兩邊同乘積分因子e-Lt，然后從tn到tn+1求積分得

對(duì)式（15）應(yīng)用梯形積分公式，得

式（16）中的指數(shù)矩陣與向量的乘積采用Kroylov 子空間的方法來求解。此方法的主要思想就是利用Kroylov子空間的元素來近似，其中具體實(shí)現(xiàn)方法參考文獻(xiàn)［11］。由于積分的部分采用梯形求積公式進(jìn)行離散，故精度為二階精度。

2 數(shù)值試驗(yàn)

不同時(shí)間下的計(jì)算結(jié)果以及爆破情況如圖3、圖4 和圖5 所示。其中，圖3 表示在t=0 時(shí)該方程的初始解，圖4 表示在t=0.03 時(shí)的解，圖5 表示在t=0.043 時(shí)的解。由圖3、圖4 和圖5 可知，在t=0時(shí)有2 個(gè)環(huán)形解，且隨著時(shí)間的增長(zhǎng)，環(huán)形解的振幅逐漸增大，而環(huán)形解的半徑逐漸減小，出現(xiàn)爆破現(xiàn)象，該結(jié)果與文獻(xiàn)［6］的結(jié)果吻合。數(shù)值結(jié)果表明，本文提出的方法可以在圓形求解區(qū)域有效地捕捉爆破解。

圖2 計(jì)算區(qū)域在r和θ方向上的網(wǎng)格

圖3 t=0時(shí)該方程的初始解

圖4 t=0.03時(shí)的解

圖5 t=0.043時(shí)的解

3 結(jié)語

本文對(duì)極坐標(biāo)下的二維非線性薛定諤方程進(jìn)行求解，將拉普拉斯算子在極坐標(biāo)下表示，在空間上采用中心差分方法將離散后的式子用Kronecker積與矩陣向量化表示，最后采用積分因子方法進(jìn)行求解。在實(shí)現(xiàn)過程中采用Kroylov子空間的方法來求解指數(shù)矩陣與向量的乘積，再利用不動(dòng)點(diǎn)迭代方法在每個(gè)單元上求解非線性方程組。數(shù)值實(shí)驗(yàn)結(jié)果表明，采用Kronecker 積與矩陣向量化來求解在極坐標(biāo)下的二維非線性薛定諤方程是有效的。