陳 悅,安 靜,劉忠敏

(貴州師范大學 數學科學學院,貴州 貴陽 550025)

0 引言

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

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

1 極坐標系下的弱形式及其離散格式

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

(1)

(2)

其中

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

(3)

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

則方程(1)在極坐標系下的等價形式為：

(4)

其內積和范數分別為:

其內積和范數分別為:

(5)

其中

定義逼近空間:

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

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

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

(8)

2 誤差估計

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

證明由邊界條件ψ(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)的特征對，則有：

(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)。

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

令

(14)

證明由算子范數的定義有：

=εNM。

證畢

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

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

(15)

(16)

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

(17)

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

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

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

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

由引理3可以得到

結合(15)可得(16)。

證畢

(18)

相應的內積和范數分別為：

由文獻[23]我們有下面的引理：

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

相應的內積和范數分別為：

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

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

進一步令

(19)

(20)

證明由引理1有

由投影算子ΠNM的性質可得：

又由于

則有

從而有

由引理4有

由引理5有

因此

進一步可得

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

3 算法的有效實現

首先,構造逼近空間中的一組基函數。令

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

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

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

令

我們將尋找

(21)

將(21)代入(8)，讓φNM取遍逼近空間XNM中的所有基函數便可得到如下的線性特征系統(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)。

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

4 非線性特征值問題

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

(22)

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

(23)

(24)

(25)

其中

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

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

(28)

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

(29)

(30)

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

5 數值實驗

為了表明算法的有效性,我們給出了一些數值算例。將在MATLAB 2016a平臺上進行編程計算。

例1：我們在方程(1)中取R=1,V(x,y)=ex+y+1。對于不同的N和M,前4個特征值的數值結果分別在表1和表2被列出。

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

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

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

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

圖2 對于不同的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個數值算例。不失一般性,我們仍然取R=1,V(x,y)=ex+y+1。對于不同的N和M,第一個特征值的數值結果在表3中被列出。

表3 對于不同的N和M，對于的第一個特征值的數值結果Tab.3 Numerical results of the first eigenvalues for different N and M

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

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

6 結論

本文針對圓域上Schr?dinger方程特征值問題提出了一種有效的譜Galerkin 逼近。首先利用極坐標變換將笛卡爾直角坐標系下的二階Schr?dinger方程特征值問題轉化為極坐標系下的一種等價形式。其次,通過引入極點條件和定義適當的帶權Sobolev空間及其逼近空間,建立了極坐標系下二階Schr?dinger方程特征值問題的變分形式及其離散格式,并對逼近解的誤差估計給出了證明,數值算例驗證了算法的有效性和高精度性。