王肇君， 吳亭亭*， 曾泰山

(1. 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 濟(jì)南 250358; 2. 華南師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院, 廣州 510631; 3. 廣東省大數(shù)據(jù)分析與處理重點(diǎn)實(shí)驗(yàn)室, 廣州 510006)

Helmholtz方程是一類重要的物理方程,在電磁學(xué)[1]、聲學(xué)[2]和地球物理[3]等科學(xué)領(lǐng)域有著廣泛的應(yīng)用。眾所周知,Helmholtz問(wèn)題的困難在于高波數(shù)問(wèn)題[4]。對(duì)于高波數(shù)的Helmholtz問(wèn)題,由于其解高度震蕩,建立其高效和高精度的數(shù)值方法是困難的[5]。因此,尋求Helmholtz方程的高效數(shù)值解法具有重要的理論意義和應(yīng)用價(jià)值。

目前,數(shù)值求解Helmholtz方程的方法主要有有限元法[5-8]、有限差分法[3,9-16]和譜元方法[17]等。其中,有限元法的精度高于有限差分法的,易于復(fù)雜邊界條件和非均勻結(jié)構(gòu)化的處理,但其計(jì)算量較大;有限差分法的計(jì)算簡(jiǎn)單、存儲(chǔ)量小,易于實(shí)現(xiàn)?；谟邢薏罘址?文獻(xiàn)[18]針對(duì)高波數(shù)問(wèn)題提出了一種并行預(yù)條件迭代法來(lái)求解 Helmholtz 方程,有效地提高了計(jì)算效率。由于高波數(shù) Helmholtz 方程的解高度震蕩,應(yīng)用有限差分法求解高波數(shù)問(wèn)題時(shí)常常受到非物理震蕩的困擾,即所謂的“數(shù)值頻散”。數(shù)值頻散通常是由計(jì)算中采用粗網(wǎng)格引起的，因此,數(shù)值求解 Helmholtz 方程時(shí),需要考慮數(shù)值頻散和計(jì)算量。

近年來(lái),求解 Helmholtz 問(wèn)題的有限差分法得到了迅速的發(fā)展。為了提高差分法的數(shù)值精度,一方面,學(xué)者們提出了帶參數(shù)的差分格式,并通過(guò)優(yōu)化參數(shù)得到高精度的差分法[3,9];另一方面,學(xué)者們采用高階差分格式來(lái)求解Helmholtz方程[10-13,15-16,19]。基于極小化數(shù)值頻散的思想,文獻(xiàn)[3,9]構(gòu)造了求解Helmholtz方程的多種二階差分格式:文獻(xiàn)[3]利用旋轉(zhuǎn)坐標(biāo)系和借助加權(quán)平均的思想，建立了一種旋轉(zhuǎn)九點(diǎn)差分格式，從而有效地抑制數(shù)值頻散,提高數(shù)值精度;文獻(xiàn)[9]基于點(diǎn)加權(quán)的思想,提出了一種帶PML的Helmholtz方程的差分格式，該格式為二階格式,且采用極小化數(shù)值頻散的思想來(lái)得到差分格式的優(yōu)化參數(shù)。緊致差分格式具有精度高、計(jì)算簡(jiǎn)單和對(duì)網(wǎng)格節(jié)點(diǎn)數(shù)要求不高等優(yōu)點(diǎn),因此受到了很多學(xué)者的關(guān)注。如:提出了多種四階緊致有限差分格式[10-11,13,15,19];構(gòu)造了二維Helmholtz方程的六階緊致有限差分格式[12,16];通過(guò)對(duì)六階緊致差分格式的截?cái)嗾`差進(jìn)行二階緊致逼近,從而得到一種高精度的緊致差分格式[12];采用9個(gè)點(diǎn)(中心節(jié)點(diǎn)及圍繞該節(jié)點(diǎn)的8個(gè)節(jié)點(diǎn))的加權(quán)平均來(lái)近似零階項(xiàng),以抑制數(shù)值頻散[16];提出了三維Helmholtz方程的六階緊致差分格式[14]。

本文主要構(gòu)造了二維Helmholtz方程的一種六階緊致差分格式。首先,給出了文獻(xiàn)[16]的六階緊致差分法的截?cái)嗾`差。然后,為了提高數(shù)值精度、減弱誤差對(duì)波數(shù)k的依賴性,對(duì)所得截?cái)嗾`差的部分項(xiàng)進(jìn)行二階緊致逼近,從而得到一種改進(jìn)的六階緊致差分格式。其次,將改進(jìn)的六階緊致差分格式的誤差方程進(jìn)行轉(zhuǎn)化,即將二維誤差方程轉(zhuǎn)化為若干個(gè)一維問(wèn)題。在一維問(wèn)題的基礎(chǔ)上,證明了該格式在一定條件下的解存在唯一,且為六階收斂的。再次,通過(guò)極小化數(shù)值頻散的思想給出差分格式優(yōu)化參數(shù)的加細(xì)選取策略。最后,給出2個(gè)數(shù)值算例來(lái)驗(yàn)證所提格式的精確性和有效性。

1 改進(jìn)的六階緊致有限差分格式

二維Helmholtz方程有如下形式

(1)

4um,n-2(um,n+1+um,n-1+um+1,n+um-1,n)]。

接下來(lái),在文獻(xiàn)[16]提出的六階緊致差分格式的基礎(chǔ)上,構(gòu)造改進(jìn)的六階緊致有限差分格式。為此,先回顧文獻(xiàn)[16]提出的六階緊致差分格式。具體地,令

其中,c,d,e是滿足c+d+e=1的參數(shù)。事實(shí)上,h(k2u)是k2u的二階近似[20]。基于Equation-based方法[21]以及加權(quán)平均的思想[20],文獻(xiàn)[16]對(duì)二維Helmholtz方程建立了如下六階緊致差分格式

(2)

其中

利用 Taylor 展開(kāi)可得該格式的截?cái)嗾`差為：

O(h8)。

(3)

易知,R滿足

(4)

(5)

其中

證明首先,易知

(6)

(7)

類似地,對(duì)Helmholtz方程(1)的兩端關(guān)于y求二階偏導(dǎo),可得

(8)

將式(7)與式(8)相結(jié)合,可得

(9)

(10)

(11)

類似地,對(duì)式(8)的兩端關(guān)于x求二階偏導(dǎo),可得

(12)

將式(11)與式(12)相加,可得

(13)

代入式(6),可得

(14)

(15)

類似地,對(duì)式(8)的兩端關(guān)于y求二階偏導(dǎo),可得

(16)

將式(15)與式(16)相加,可得

(17)

(18)

(19)

類似地,對(duì)式(16)兩端關(guān)于x求二階偏導(dǎo),可得

(20)

將式(19)與式(20)相結(jié)合,可得

代入式(14),可得

(21)

(22)

對(duì)式(16)兩端關(guān)于y求二階偏導(dǎo),可得

(23)

將式(22)與式(23)相加,可得

(24)

(25)

將式(6)、(10)、(14)、(18)、(21)、(25)代入截?cái)嗾`差R(3),即可得到命題的結(jié)論。證畢。

將截?cái)嗾`差(5)代入式(4)，可得:

(26)

其中

(27)

(28)

2 差分格式的收斂性分析

本節(jié)研究改進(jìn)的六階緊致有限差分格式(28)的解的存在性、唯一性和收斂性:首先，將二維誤差方程轉(zhuǎn)化為若干個(gè)一維問(wèn)題;然后,證明了當(dāng)k2不是-Δ的特征值且kh充分小時(shí),差分格式(28)的解存在唯一;進(jìn)一步地,證明了差分格式(28)為六階收斂的。

為便于表述,考慮在方形區(qū)域Ω∶=(0,1)×(0,1)內(nèi)進(jìn)行收斂性分析。首先,引入一些記號(hào)。令步長(zhǎng)h∶=1/N,N為正整數(shù),并記N∶={1,2,…,N-1}。對(duì)于向量V=[Vm:mN],令

對(duì)矩陣W=[Wm,n:m,nN],定義其離散L2- 范數(shù)為：

(29)

對(duì)任意nN,定義Wn為矩陣W的第n個(gè)列向量,且有

其次,為了進(jìn)行收斂性分析,給出差分格式(28)的誤差方程。對(duì)任意m,nN,令Em,n∶=Um,n-um,n。將式(28)與式(26)相減,得到差分格式(28)的誤差方程:

(30)

接下來(lái),為了分析誤差方程(30),將誤差方程(30)轉(zhuǎn)化為若干個(gè)一維問(wèn)題。為此,引入矩陣D(N-1)×(N-1):

顯然,三對(duì)角矩陣D是對(duì)稱正定的,且D的特征值如下:

相應(yīng)的特征向量為rn∶=[rm,n:mN]T,其中

(31)

事實(shí)上,特征向量組{rj,jN}是N-1空間的一組正交基。以下引理將誤差方程(30)轉(zhuǎn)化為N-1個(gè)一維問(wèn)題。

引理1基于特征向量組{rj,jN},誤差方程(30)可轉(zhuǎn)化為N-1個(gè)一維問(wèn)題:對(duì)任意jN,有

(32)

其中:

(33)

借助特征向量的分量表達(dá)式(31),由式(33)知,對(duì)任意m,nN,有:

(34)

下面考慮一維問(wèn)題(32)的解的存在唯一性,并給出收斂性分析。

引理2假設(shè)k2不是-Δ的特征值。對(duì)于任意jN,當(dāng)kh充分小時(shí),一維問(wèn)題(32)的解存在且唯一,且有以下估計(jì)式

(35)

其中,C為某一常數(shù)。

其中,Hj是一維問(wèn)題(32)的系數(shù)矩陣。由于Hj是三對(duì)角托普利茨矩陣,由文獻(xiàn)[22]可知

(36)

所以,當(dāng)ξ=0時(shí),有det(Hj)≠0。當(dāng)ξ>0時(shí),由式(36)知det(Hj)≠0。接下來(lái),證明當(dāng)ξ<0且kh充分小時(shí),det(Hj)的值不為零。當(dāng)ξ<0時(shí),由式(36)得:

因此,可以得到

其中,C1為正常數(shù)。因此

(37)

(38)

聯(lián)立式(37)和式(38),即得引理的結(jié)論。證畢。

最后,給出差分格式(28)的解的收斂性分析。

定理1假設(shè)k2不是-Δ的特征值。當(dāng)kh充分小時(shí),差分格式(28)的解U存在且唯一。若在狄利克雷或諾伊曼邊界條件下的Helmholtz方程(1)的解u在求解區(qū)域Ω上滿足u則有

‖|u-U|‖≤Ch6,

其中,U∶=[Um,n:m,nN],u∶=[um,n:m,nN],C是依賴于u的常數(shù)。

證明下面只給出狄利克雷邊界條件下的證明,在諾伊曼條件下的證明類似可得。為證明差分格式(28)的解存在且唯一,首先分析誤差方程(30)。引理1說(shuō)明了誤差方程(30)可以轉(zhuǎn)化為N-1個(gè)一維問(wèn)題(32)。由引理2可知,當(dāng)kh充分小時(shí),一維問(wèn)題(32)的解存在且唯一。因此,當(dāng)kh充分小時(shí),誤差方程(30)的解也是存在唯一的。所以,當(dāng)kh充分小時(shí),差分格式(28)的解存在且唯一。

令E∶=[Em,n:m,nN],m,n:m,nN]。為了對(duì)差分格式(28)進(jìn)行收斂性分析,只需估計(jì)‖|E‖。首先,根據(jù)范數(shù)的定義和內(nèi)積的性質(zhì),有

(39)

(40)

對(duì)式(39)的右端應(yīng)用式(35),并結(jié)合式(40),可得

(41)

定理得證。

由截?cái)嗾`差的表達(dá)式(27)知,當(dāng)kh充分小時(shí),差分格式(28)的誤差僅依賴于函數(shù)u,對(duì)波數(shù)k的依賴性較弱。

3 差分格式的參數(shù)優(yōu)化選取策略

本節(jié)給出差分格式(28)的頻散方程。從頻散方程出發(fā),基于極小化數(shù)值頻散的思想,給出差分格式(28)的優(yōu)化參數(shù)的加細(xì)選擇策略。

為了得到優(yōu)化參數(shù)d和e,首先給出差分格式(28)的頻散方程。為此,將差分格式(28)進(jìn)行改寫。為便于分析,引入以下記號(hào):

(

W

)

m,n

∶=

W

m+1,n+1

+

W

m-1,n-1

+

W

m+1,n-1

+

W

m-1,n+1

,

(

W

)

m,n

∶=

W

m+1,n

+

W

m-1,n

+

W

m,n-1

+

W

m,n+1

。

從而將差分格式(28)改寫為:

A0Um,n+A(U)m,n+A

(42)

其中,

A∶=

接下來(lái),對(duì)差分格式(28)進(jìn)行頻散分析。為此,令g=0,并將Um,n∶=e-ik(xmcos θ+ynsin θ)代入差分方程(42),利用歐拉公式得頻散方程:

4APQ+2A(P+Q)+A0=0,

(43)

其中,P∶=cos(khcosθ),Q∶=cos(khsinθ)。

標(biāo)準(zhǔn)化的數(shù)值相速度和標(biāo)準(zhǔn)化的數(shù)值群速度是衡量數(shù)值頻散的重要標(biāo)準(zhǔn)[3]。在實(shí)際應(yīng)用中,常用標(biāo)準(zhǔn)化的數(shù)值相速度來(lái)作為衡量數(shù)值頻散的標(biāo)準(zhǔn)。對(duì)于有限差分法,當(dāng)其標(biāo)準(zhǔn)化的數(shù)值相速度越接近1時(shí),其數(shù)值頻散越小,相應(yīng)的數(shù)值精度越高。由文獻(xiàn)[20]可知，要得到標(biāo)準(zhǔn)化數(shù)值相速度kN/k,需要獲得數(shù)值波數(shù)kN。然而,對(duì)于差分格式(28),很難獲得其數(shù)值波數(shù)kN。因此,借鑒文獻(xiàn)[16]的方法,直接從頻散方程出發(fā)來(lái)獲得差分格式的優(yōu)化參數(shù)。

為了得到差分格式(28)的優(yōu)化參數(shù),再次考慮頻散方程(43)。令v為波傳播的相速度,λ為波長(zhǎng),G∶=λ/h為每個(gè)波長(zhǎng)的網(wǎng)格節(jié)點(diǎn)數(shù)。由λ=2πv/ω,k=ω/v，可得kh=2π/G。在頻散方程(43)兩端同時(shí)乘以h2,并用2π/G替代kh,可得:

(44)

加細(xì)選擇策略:

步驟1： 估計(jì)區(qū)間IG∶=[Gmin,Gmax],其中Gmin代表G的最小值,Gmax代表G的最大值。令

其中，n=1,2,…,s；s為給定的正整數(shù)。

步驟3： 將1/G和θ代入方程(44),得到一個(gè)超定方程組。

步驟4： 應(yīng)用最小二乘法求解得到的方程組,即可得到差分格式(28)的優(yōu)化參數(shù)。

4 數(shù)值算例

下面將通過(guò)2個(gè)數(shù)值算例來(lái)驗(yàn)證本文提出的改進(jìn)六階緊致差分格式的有效性。在本節(jié)中,數(shù)值解與精確解之間的誤差用離散L2-范數(shù)(式(29))衡量。根據(jù)經(jīng)驗(yàn)將區(qū)間[2.5,400]劃分為7個(gè)小區(qū)間，基于加細(xì)選擇策略，對(duì)于部分區(qū)間IG，得到差分格式(28)的多組加細(xì)的優(yōu)化參數(shù)(表1)。在下面的計(jì)算中，均采用表1中的參數(shù)作為差分格式(28)的優(yōu)化參數(shù)。將本文提出的差分格式(28)(簡(jiǎn)記為IC6)與文獻(xiàn)[16]構(gòu)造的帶優(yōu)化參數(shù)的六階緊致差分格式(簡(jiǎn)記為OC6)進(jìn)行了數(shù)值精度的對(duì)比。

表1 加細(xì)的優(yōu)化參數(shù)Table 1 The refined optimal parameters

問(wèn)題1考慮二維 Helmholtz 方程

Δu+k2u=-a2sin(ax)sin(ky),(x,y)Ω∶=(0,π)×(0,π),

其邊界條件為

u(x0,y0)=0,(x0,y0)?！??Ω。

假設(shè)a、k為正整數(shù)，該問(wèn)題的真解為

u(x,y)∶=sin(ax)sin(ky),(x,y)[0,π]×[0,π]。

a=1,k=24;a=3,k=24 和a=3,k=48 時(shí),不同網(wǎng)格剖分下所對(duì)應(yīng)的 2 種差分格式的離散L2- 范數(shù)誤差 (表2至表4)表明:與OC6相比,IC6的數(shù)值精度有較大的提高;在參數(shù)a相同的情況下,IC6的數(shù)值精度對(duì)波數(shù)k的依賴性弱于OC6的,與前面的理論分析結(jié)果一致。a=1,N=256時(shí),不同波數(shù)k所對(duì)應(yīng)的2種差分格式的離散L2-范數(shù)誤差(表5)表明:當(dāng)網(wǎng)格剖分N固定,隨著波數(shù)k的增大,OC6與IC6的離散L2-范數(shù)誤差都有所增大,但是與OC6 相比,IC6 的數(shù)值精度依然有顯著的提高。此外,離散L2-范數(shù)誤差曲線圖(圖1)表明了2種差分格式的誤差曲線與斜率為-6的直線近似平行，說(shuō)明2種差分格式均為六階格式。

表2 a=1,k=24 時(shí)的離散 L2- 范數(shù)誤差Table 2 The error in the discreteL2- norm for a=1,k=24

表3 a=3,k=24 時(shí)的離散 L2- 范數(shù)誤差Table 3 The error in the discreteL2- norm for a=3,k=24

表4 a=3,k=48 時(shí)的離散 L2- 范數(shù)誤差Table 4 The error in the discreteL2- norm for a=3,k=48

表5 a=1,N=256 時(shí)的離散 L2- 范數(shù)誤差Table 5 The error in the discreteL2- norm for a=1,N=256

圖1 問(wèn)題1的離散 L2- 范數(shù)誤差曲線

問(wèn)題2考慮二維 Helmholtz 方程

Δu+k2u=0,(x,y)Ω∶=(0,π)×(0,π),

其邊界條件為

假設(shè)k為常數(shù),則上述問(wèn)題的真解為u=eiβysinx,且有 12+β2=k2。在y=π 處使用 Sommerfeld 型邊界條件uy-iβu=0的原因是:將特征值移到復(fù)平面,從而消除在實(shí)特征值處發(fā)生共振的可能性[23]。因此,對(duì)任意k>0,該問(wèn)題的解存在且唯一。

β=100,k=100.005 0;β=199.997 5,k=200和β=499.999 0,k=500時(shí)不同網(wǎng)格剖分下所對(duì)應(yīng)的2種差分格式的離散L2-范數(shù)誤差(表6至表8)表明:與OC6相比,IC6的誤差明顯減小,說(shuō)明IC6的數(shù)值精度高于OC6的。N=700時(shí),不同波數(shù)k所對(duì)應(yīng)的2種差分格式的離散L2-范數(shù)誤差(表9)表明:隨著波數(shù)k的增大,OC6和IC6的離散L2-范數(shù)誤差都有所增大,但是IC6的數(shù)值精度依然高于OC6的。此外,離散L2-范數(shù)誤差曲線圖(圖2)表明2種差分格式的誤差曲線與斜率為-6的直線近似平行，說(shuō)明2種差分格式均為六階格式。

表6 β=100,k=100.005 0時(shí)的離散 L2- 范數(shù)誤差Table 6 The error in the discreteL2- norm for β=100,k=100.005 0

表7 β=199.997 5,k=200 時(shí)的離散 L2- 范數(shù)誤差Table 7 The error in the discreteL2- norm for β=199.997 5,k=200

表8 β=499.999 0,k=500 時(shí)的離散 L2- 范數(shù)誤差Table 8 The error in the discreteL2- norm for β=499.999 0,k=500

表9 N=700 時(shí)的離散 L2- 范數(shù)誤差Table 9 The error in the discreteL2- norm for N=700

圖2 問(wèn)題 2 的離散 L2- 范數(shù)誤差曲線

5 結(jié)論

本文建立了二維 Helmholtz 方程的一種改進(jìn)的六階緊致差分格式,該格式的特點(diǎn)是其誤差較弱地依賴波數(shù)k,收斂性分析表明該格式是六階收斂的。數(shù)值算例驗(yàn)證了該格式可以有效地提高數(shù)值精度,抑制數(shù)值頻散。