陳 悅,安 靜,劉忠敏

(貴州師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,貴州 貴陽 550025)

0 引言

薛定諤方程特征值問題具有重要的物理背景,它在原子物理、核物理和計(jì)算量子化學(xué)等領(lǐng)域有著廣泛地應(yīng)用,例如,在非線性彈性框架下機(jī)械結(jié)構(gòu)振動(dòng)模的計(jì)算、描述玻色-愛因斯坦凝聚穩(wěn)態(tài)結(jié)構(gòu)的Gross-Pitaevskii方程[1-4]、以及用于計(jì)算量子化學(xué)和材料科學(xué)中分子系統(tǒng)基態(tài)電子結(jié)構(gòu)的Hartree-Fock和Kohn-Sham方程[5-7]。關(guān)于薛定諤方程特征值問題的理論分析和數(shù)值計(jì)算已經(jīng)有很多成果[8-15],但它們主要都是基于一些低階的有限元方法,要獲得一些高精度的數(shù)值結(jié)果,將會(huì)花費(fèi)很多內(nèi)存容量和計(jì)算時(shí)間,尤其是對(duì)一些特殊區(qū)域(如圓域，球域等)上的非線性薛定諤方程特征值問題。眾所周知,譜方法是一類非常重要的數(shù)值方法,由于其具有譜精度的特點(diǎn),我們一般只需要較少的自由度就能獲得較高的精度,但其計(jì)算區(qū)域要求是矩形或立方體區(qū)域。最近，文獻(xiàn)[16-20]提出了圓域上二階、四階方程及其特征值問題有效的譜方法，但這些方法都是基于常系數(shù)或徑向變系數(shù)的情況。另外，文獻(xiàn)[21]提出了無界域上三維薛定諤方程基于降維格式的一種譜方法,該方法也是基于徑向變系數(shù)的情況,由于薛定諤方程特征值問題一般都是指數(shù)衰減的,我們通常把無界域截?cái)酁橐粋€(gè)圓域(二維)或球域(三維),那么如何提出圓域上帶有一般變系數(shù)的薛定諤方程特征值問題的譜方法將是有意義的事情。

因此,本文的目的是提出了圓域上二階變系數(shù)Schr?dinger方程特征值問題的一種有效的譜伽遼金方法。首先利用極坐標(biāo)變換將笛卡爾直角坐標(biāo)系下的二階Schr?dinger方程特征值問題轉(zhuǎn)化為極坐標(biāo)系下的一種等價(jià)形式。其次,極條件被推導(dǎo),克服了極點(diǎn)奇性引入的困難。再結(jié)合特征函數(shù)的邊界條件和在θ方向的周期性,我們定義了帶權(quán)的Sobolev空間及其逼近空間,建立了極坐標(biāo)系下二階Schr?dinger方程特征值問題的一種弱形式和相應(yīng)的離散格式。基于緊算子的譜理論、非一致帶權(quán)Sobolev空間中投影算子的逼近性質(zhì)以及傅里葉基函數(shù)的逼近性質(zhì),我們對(duì)逼近解的誤差估計(jì)給出了證明。最后,我們給出了一些數(shù)值實(shí)驗(yàn),數(shù)值結(jié)果表明我們的算法是有效的和高精度的。

1 極坐標(biāo)系下的弱形式及其離散格式

作為一個(gè)模型問題，我們首先考慮下面的二階變系數(shù)薛定諤方程特征值問題：

(1)

(2)

其中

Δu(x,y)=Δψ(t,θ)

(3)

為了使(3)有意義,ψ(t,θ)需要滿足下面的本質(zhì)極條件：

則方程(1)在極坐標(biāo)系下的等價(jià)形式為：

(4)

其內(nèi)積和范數(shù)分別為:

其內(nèi)積和范數(shù)分別為:

(5)

其中

定義逼近空間:

XNM=span{φmN(t)eimθ:φmN(t)∈PmN,-M≤m≤M},

其中PmN={p∈PN:mp(-1)=p(1)=0},PN為次數(shù)不超過N的多項(xiàng)式空間。則(5)的離散格式為：找λNM∈C,ψNM∈XNM,使得

A(ψNM,φNM)=λNMB(ψNM,φNM),?φNM∈XNM

(8)

2 誤差估計(jì)

為了敘述方便,我們用ab表示a≤Cb,其中C為與M,N無關(guān)的正常數(shù)。

證明由邊界條件ψ(1,θ)=0和Cauchy-Schwarz不等式，有：

則有：

將上面不等式兩邊在區(qū)域D上積分可得：

證畢

|A(ψ,φ)|‖ψ‖1,*w‖φ‖1,*w,

=‖ψ‖1,*w‖φ‖1,*w,

類似于定理1的證明，我們有下面的定理：

|B(ψ,φ)|‖ψ‖w‖φ‖w,

(9)

A(TNMψ,φNM)=B(ψ,φNM),?φNM∈XNM

(10)

(11)

引理2 令T和TNM是分別由(9)和(10)定義的有界線性算子，則有

TNM=ΠNMT。

A(ΠNMTψ-TNMψ,φNM)=A(ΠNMTψ-Tψ,φNM)+A(Tψ-TNMψ,φNM)=0

(12)

在(12)中取φNM=ΠNMTψ-TNMψ得到：

A(ΠNMTψ-TNMψ,ΠNMTψ-TNMψ)=0。

由定理1可得：

TNM=ΠNMT。

顯然

TNM|XNM:XNM→XNM。

引理3 令(λ,ψ)和(λNM,ψNM)分別為弱形式(5)和離散格式(8)的特征對(duì)，則有：

(13)

證明由(5)和(8)式我們有

A(ψNM-ψ,ψNM-ψ)-λB(ψNM-ψ,ψNM-ψ)

=A(ψNM,ψNM)-2A(ψNM,ψ)+A(ψ,ψ)-λB(ψNM,ψNM)+2λB(ψNM,ψ)-λB(ψ,ψ)

=λNMB(ψNM,ψNM)-2λB(ψNM,ψ)+λB(ψ,ψ)-λB(ψNM,ψNM)+2λB(ψNM,ψ)-λB(ψ,ψ)

=λNMB(ψNM,ψNM)-λB(ψNM,ψNM)=(λNM-λ)B(ψNM,ψNM)。

將上面等式兩邊同時(shí)除以B(ψNM,ψNM)可得到(13)。

令

(14)

證明由算子范數(shù)的定義有：

=εNM。

證畢

設(shè)S(λ)和S(λNM)分別表示(5)和(8)式中λ和λNM相應(yīng)的特征函數(shù)空間。

定理4令(λ,ψ)和(λNM,ψNM)分別為(5)和(8)式的特征對(duì)，則有

(15)

(16)

‖ψ-ψNM‖A‖(T-TNM)|S(λ)‖A

(17)

對(duì)于ψ∈S(λ),‖ψ‖A=1,有

由上面的兩個(gè)等式和(17)有：

‖ψ-ψNM‖A‖(T-TNM)|S(λ)‖A

對(duì)于任意ψNM=S(λNM),‖ψNM‖A=1,則有

由引理3可以得到

結(jié)合(15)可得(16)。

證畢

(18)

相應(yīng)的內(nèi)積和范數(shù)分別為：

由文獻(xiàn)[23]我們有下面的引理：

‖f(θ)-fM(θ)‖kMk-s|f(θ)|s,

相應(yīng)的內(nèi)積和范數(shù)分別為：

則由文獻(xiàn)[24]中的定理1.8.2，我們有下面的引理：

由于ψ(t,θ)在θ方向上是以2π為周期，則由傅里葉基函數(shù)展開有：

進(jìn)一步令

(19)

(20)

證明由引理1有

由投影算子ΠNM的性質(zhì)可得：

又由于

則有

從而有

由引理4有

由引理5有

因此

進(jìn)一步可得

將上式代入(15)式得(19)，再結(jié)合(19)得(20)，證畢

3 算法的有效實(shí)現(xiàn)

首先,構(gòu)造逼近空間中的一組基函數(shù)。令

φmi(t)=Li(t)-Li+2(t),i=0,1,…,N-2,

其中Li(t)是次數(shù)為i的Legendre 多項(xiàng)式，則逼近空間XNM為：

XNM=span{φmk(t)eimθ:-M≤m≤M,k=0,1,…,N-1-sign(|m|)}。

令

我們將尋找

(21)

將(21)代入(8)，讓?duì)誑M取遍逼近空間XNM中的所有基函數(shù)便可得到如下的線性特征系統(tǒng)：

(A+B+Q)U=λNMCU。

其中：

U=[u-M,0,u-M,1,…,u-M,N-2,…,u00,u01,…,u0,N-1,…,uM,0,uM,1,…,uM,N-2]T,

A=(akjnm),B=(bkjnm),C=(ckjnm),Q=(qkjnm)。

由勒讓德多項(xiàng)式和傅里葉基函數(shù)的正交性質(zhì)可知，矩陣A,B,C都是稀疏的分塊帶狀矩陣，矩陣Q的稀疏性依賴于變系數(shù)V(x,y)的性質(zhì)。

4 非線性特征值問題

在這一節(jié)，將提出的算法應(yīng)用于非線性特征值問題的計(jì)算,考慮下面的非線性特征值問題：

(22)

類似于(4)的推導(dǎo)，我們可得到方程(22)在極坐標(biāo)系下的等價(jià)形式為：

(23)

(24)

(25)

其中

則(24)-(25)的離散形式為：找(λNM,ψNM)∈C×XNM,使得

A(ψNM,φNM)=λNMB(ψNM,φNM),?φNM∈XNM,

(28)

由于(28)是非線性的,我們需要通過迭代法進(jìn)行求解，我們建立了如下的迭代格式:

(29)

(30)

因此,把非線性格式(28)轉(zhuǎn)化為變系數(shù)的迭代格式(29)和(30),從而可用第3節(jié)提出的算法有效地求解。

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

為了表明算法的有效性,我們給出了一些數(shù)值算例。將在MATLAB 2016a平臺(tái)上進(jìn)行編程計(jì)算。

例1：我們?cè)诜匠?1)中取R=1,V(x,y)=ex+y+1。對(duì)于不同的N和M,前4個(gè)特征值的數(shù)值結(jié)果分別在表1和表2被列出。

表1 N=15和不同的M情況下前4個(gè)逼近特征值的結(jié)果Tab.1 Numerical results of the top four eigenvalues for N=15 and different M

表2 M=8和不同的N情況下前4個(gè)逼近特征值的結(jié)果Tab.2 Numerical results of the top four eigenvalues for M=8 and different N

從表1、表2可知,當(dāng)固定N=15,M≥7和固定M=8,N≥11時(shí),前4個(gè)逼近特征值達(dá)到了至少12 位有效數(shù)字的精度。另外,我們以M=30,N=60時(shí)的數(shù)值解作為參考解,在圖1中畫出了數(shù)值解和參考解之間的誤差曲線,從圖1可以觀察到我們提出的算法是收斂的和高精度的。為了進(jìn)一步直觀地描述逼近特征值的收斂率，我們?cè)趫D2中作出了log-log尺度下逼近特征值的誤差曲線，從圖2中可以觀察到我們的算法是指數(shù)收斂的。

圖1 對(duì)于不同的M(左)和N(右)，數(shù)值解與參考解之間的誤差圖像Fig.1 The error figures between the reference solutions and approximate solutions for different M(left) and different N(right)

圖2 對(duì)于不同的M(左)和N(右)，逼近解與參考解在log-log尺度下的誤差曲線Fig.2 The error figures on the log-log scale between the reference solutions and approximate solutions for different M(left) and different N(right)

例2：我們以Schr?dinger方程(22)作為第2個(gè)數(shù)值算例。不失一般性,我們?nèi)匀蝗=1,V(x,y)=ex+y+1。對(duì)于不同的N和M,第一個(gè)特征值的數(shù)值結(jié)果在表3中被列出。

表3 對(duì)于不同的N和M，對(duì)于的第一個(gè)特征值的數(shù)值結(jié)果Tab.3 Numerical results of the first eigenvalues for different N and M

從表3可知,當(dāng)M≥6,N≥11時(shí),第一個(gè)逼近特征值達(dá)到至少13 位有效數(shù)字的精度。類似地,我們也選擇M=30,N=60時(shí)的數(shù)值解作為參考解,在圖3中給出了數(shù)值解和參考解之間的誤差圖像,從圖3我們可以觀察到我們提出的算法也是收斂的和高精度的。

圖3 對(duì)于不同的M(左)和N(右),數(shù)值解與參考解之間的誤差圖像Fig.3 The error figures between the reference solutions and approximate solutions for different M(left) and different N(right)

6 結(jié)論

本文針對(duì)圓域上Schr?dinger方程特征值問題提出了一種有效的譜Galerkin 逼近。首先利用極坐標(biāo)變換將笛卡爾直角坐標(biāo)系下的二階Schr?dinger方程特征值問題轉(zhuǎn)化為極坐標(biāo)系下的一種等價(jià)形式。其次,通過引入極點(diǎn)條件和定義適當(dāng)?shù)膸?quán)Sobolev空間及其逼近空間,建立了極坐標(biāo)系下二階Schr?dinger方程特征值問題的變分形式及其離散格式,并對(duì)逼近解的誤差估計(jì)給出了證明,數(shù)值算例驗(yàn)證了算法的有效性和高精度性。