牟宴銘，安 靜*

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

0 引言

特征值問題在流體力學(xué)、結(jié)構(gòu)力學(xué)、隨機過程和量子力學(xué)等領(lǐng)域有重要的物理背景和廣泛的應(yīng)用[1-3]，其數(shù)值計算方法一直是很多學(xué)者關(guān)注和研究的課題.到目前為止，已經(jīng)有很多數(shù)值方法求解特征值問題，主要包括有限差分法[4-5]、有限元法[6-8]、譜方法[9-11]等.然而，對于一些特殊的計算區(qū)域，如球域、橢球域、柱體區(qū)域等，直接利用差分法或有限元法求解是比較困難的.主要原因是這些區(qū)域都帶有曲面邊界，區(qū)域剖分比較復(fù)雜，邊界附近單元上的差分逼近或有限元基函數(shù)的構(gòu)造也比較復(fù)雜.另外，由于計算區(qū)域是三維區(qū)域，要獲得較高精度的數(shù)值解需要較高的時間復(fù)雜度和大量的內(nèi)存空間.

近年來，譜方法已經(jīng)成為計算微分方程數(shù)值解的一類重要方法[12-14]，將譜方法應(yīng)用于這些特殊區(qū)域上的特征值問題是一件有意義的事情.文中提出了三維柱體域上二階橢圓特征值問題的一種新的譜逼近方法.首先，將直角坐標系下的二階橢圓特征值問題轉(zhuǎn)化為柱坐標系下的等價形式，再利用變量分離方法將原問題轉(zhuǎn)化為矩形區(qū)域上的一系列二維特征值問題；其次，針對實心圓柱體和空心圓柱體兩種情況，分別引入了兩種Sobolev空間和相應(yīng)的多項式逼近空間，對每個降維的二維特征值問題建立了變分形式和離散格式；然后，利用全連續(xù)算子的譜理論和非一致帶權(quán)Sobolev空間中投影算子的逼近性質(zhì)證明了逼近解的誤差估計.另外，通過構(gòu)造逼近空間的一組有效基函數(shù)，推導(dǎo)了離散格式基于張量積的矩陣形式.最后，給出了一些數(shù)值算例，數(shù)值模擬結(jié)果表明我們的算法是有效的和高精度的.

1 降維格式和極條件

作為一個模型問題，考慮二階橢圓特征值問題

(1)

(2)

其中Ω是實心圓柱體或空心圓柱體區(qū)域，β是一個非負常數(shù).下面將推導(dǎo)問題(1)～(2)基于柱坐標變換的降維格式.

x=rcosθ,y=rsinθ,z=z,

其中R1

首先考慮R1=0的情況，問題(1)～(2)可重述為

(4)

(5)

其中(r,θ,z)∈G,G=(0,R2)×[0,2π)×(-1,1),φ(r,θ,z)=u(rcosθ,rsinθ,z).由于φ在θ方向上是2π周期的，所以由傅里葉基函數(shù)展開可得

(6)

將(6)式代入(3)式，有

由于柱坐標變換在r=0處引入了奇性，為了保證方程(4)的適定性，φm(r,z)需要滿足極條件

m2φm(r,z)|r=0=0,

即

φm(0,z)=0,m≠0.

將(6)式代入(4)式，利用傅里葉基函數(shù)的正交性可以得到一系列等價的二維特征值問題：

(8)

(9)

(10)

其中(r,z)∈(0,R2)×(-1,1).令

則問題(8)～(10)進一步化為下面等價形式：

其中(t,z)∈D,D=(-1,1)×(-1,1)=I2.

類似地，當(dāng)R1>0時，有

(14)

(15)

其中(r,θ,z)∈W,W=(R1,R2)×[0,2π)×(-1,1),θ∈[0,2π),φ(r,θ,z)=u(rcosθ,rsinθ,z).類似于R1=0的推導(dǎo)，可以得到一系列等價的二維特征值問題

(16)

(17)

其中(r,z)∈(R1,R2)×(-1,1).令

則問題(16)～(17)進一步化為如下等價形式：

其中(t,z)∈D.

2 弱形式和離散格式

僅考慮m≥0的情況，m<0的情況可以類似地推導(dǎo).為了推導(dǎo)問題(11)～(13)和問題(18)～(19)的弱形式和離散格式，需要引入下面的Sobolev空間.對于R1=0的情況，定義通常的帶權(quán)Sobolev空間

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

其中ω=1+t是權(quán)函數(shù).引入帶權(quán)Sobolev空間

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

其中

對于R1>0的情況，定義Sobolev空間

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

進一步,引入Sobolev空間

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

其中

3 誤差估計

本節(jié)證明逼近特征值和特征向量的誤差估計.為了簡潔起見，僅對R1>0的情況給出證明，R1=0的情況可以類似推導(dǎo).用ab表示a≤Mb，其中M為獨立于N且大于零的常數(shù).

證明由于-1

2R1

從而

由Cauchy-Schwarz不等式有

從而有

am(um,um)|um||um.】

類似于定理1的證明，有下面結(jié)論.

定理2bm(um,vm)是定義在L2(D)×L2(D)上的正定連續(xù)的雙線性泛函，即

引理1令Γ和ΓN分別為(24)式和(26)式定義的有界線性算子，則ΓN=ΠNΓ.

取vmN=ΠNΓum-ΓNum，并將其代入(29)式可得

am(ΠNΓum-ΓNum,ΠNΓum-ΓNum)=0.

由定理1有ΓN=ΠNΓ.】

引理2設(shè)(λm,um)和(λmN,umN)分別為問題(21)和(23)的特征對，則

證明由問題(21)和(23)有

兩邊同時除以bm(umN,umN)可得(30)式.】

令

證明根據(jù)算子范數(shù)的定義，可以得到

設(shè)R(λm)和R(λmN)分別是問題(21)和(23)的特征值λm和λmN相應(yīng)的特征函數(shù)空間.

定理4設(shè)(λm,um)和(λmN,umN)分別為問題(21)和(23)的特征對，則下列不等式成立：

||um-umN||am≤M||(Γ-ΓN)|R(λm)||am.(35)

對于um∈R(λm),||um||am=1，我們有

由上面兩個等式及(35)式可得(33)式.

令umN∈R(λmN),||umN||am=1，則

由引理2可得

結(jié)合(33)式可得(34)式.】

對于α≤-1,β≤-1的情況，函數(shù)ωα,β(x)不在L1(I)中，因此它不能用作一般的權(quán)函數(shù).在經(jīng)典的Jacobi多項式中，權(quán)函數(shù)的指數(shù)α,β>-1.為了應(yīng)用方便，將經(jīng)典的Jacobi多項式中權(quán)函數(shù)的指數(shù)α,β中的一個或兩個擴展為整數(shù)是很有必要的.因此，根據(jù)文獻[16]的思想，定義帶整數(shù)指數(shù)對(l,k)的廣義Jacobi多項式

其中n0≤n，且

定義d維Jacobi多項式的張量和權(quán)函數(shù)分別為

定義2維N次多項式空間：

引入非一致帶權(quán)Sobolev空間

其半范數(shù)和范數(shù)分別為

由文獻[17]中的定理8.1和注記8.14，有下面結(jié)果.

證明根據(jù)定理1和(28)式，有

由定理1知，||um-ΠNum||1,D||um-ΠNum||am，所以

從而

||um-ΠNum||am||um-vmN||1,D.

由vmN的任意性，有

由Poincaré不等式可得

又由于

所以

再由引理3可知

結(jié)合文獻[17]中的(3.5.32)式可知(36)式和(37)式成立.】

4 算法的有效實現(xiàn)

為了有效地求解問題(22)和(23)，首先構(gòu)造逼近空間中一組適當(dāng)?shù)幕瘮?shù)，即

其中Li是i次Legendre多項式.因此

令

當(dāng)R1=0,m=0時，令

當(dāng)R1=0,m≥1時，令

(40)

對于R1>0的情形，令

其中

由Legendre多項式的正交性可知，(39),(41)和(43)式中的剛度矩陣和質(zhì)量矩陣都是稀疏矩陣.

5 數(shù)值實驗

為了驗證算法的有效性和高精度，在MATLAB 2016a平臺上進行數(shù)值模擬.

例1取R1=0,R2=1,β=1,m=0,1.對于不同的m和N，表1和表2分別列出了前4個特征值的數(shù)值結(jié)果.

表1 當(dāng)m=0時，對應(yīng)不同N的前4個特征值

表2 當(dāng)m=1時，對應(yīng)不同N的前4個特征值

例2取R1=1,R2=2,β=1,m=0,1.對于不同的m和N，表3和表4列出了前4個特征值的數(shù)值結(jié)果.

總之，數(shù)值模擬實驗結(jié)果說明，文中算法是有效的，并且具有很高的精度.

圖1 數(shù)值解與參考解的誤差

表3 當(dāng)m=0時，對應(yīng)不同N的前4個特征值

表4 當(dāng)m=0時，對應(yīng)不同N的前4個特征值

圖2 數(shù)值解與參考解的誤差