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

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

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

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

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

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

1 改進的六階緊致有限差分格式

二維Helmholtz方程有如下形式

(1)

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

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

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

(2)

其中

利用 Taylor 展開可得該格式的截斷誤差為：

O(h8)。

(3)

易知,R滿足

(4)

(5)

其中

證明首先,易知

(6)

(7)

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

(8)

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

(9)

(10)

(11)

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

(12)

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

(13)

代入式(6),可得

(14)

(15)

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

(16)

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

(17)

(18)

(19)

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

(20)

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

代入式(14),可得

(21)

(22)

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

(23)

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

(24)

(25)

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

將截斷誤差(5)代入式(4)，可得:

(26)

其中

(27)

(28)

2 差分格式的收斂性分析

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

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

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

(29)

對任意nN,定義Wn為矩陣W的第n個列向量,且有

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

(30)

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

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

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

(31)

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

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

(32)

其中:

(33)

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

(34)

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

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

(35)

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

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

(36)

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

因此,可以得到

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

(37)

(38)

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

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

定理1假設(shè)k2不是-Δ的特征值。當(dāng)kh充分小時,差分格式(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說明了誤差方程(30)可以轉(zhuǎn)化為N-1個一維問題(32)。由引理2可知,當(dāng)kh充分小時,一維問題(32)的解存在且唯一。因此,當(dāng)kh充分小時,誤差方程(30)的解也是存在唯一的。所以,當(dāng)kh充分小時,差分格式(28)的解存在且唯一。

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

(39)

(40)

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

(41)

定理得證。

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

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

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

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

(

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∶=

接下來,對差分格式(28)進行頻散分析。為此,令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]。在實際應(yīng)用中,常用標(biāo)準(zhǔn)化的數(shù)值相速度來作為衡量數(shù)值頻散的標(biāo)準(zhǔn)。對于有限差分法,當(dāng)其標(biāo)準(zhǔn)化的數(shù)值相速度越接近1時,其數(shù)值頻散越小,相應(yīng)的數(shù)值精度越高。由文獻[20]可知，要得到標(biāo)準(zhǔn)化數(shù)值相速度kN/k,需要獲得數(shù)值波數(shù)kN。然而,對于差分格式(28),很難獲得其數(shù)值波數(shù)kN。因此,借鑒文獻[16]的方法,直接從頻散方程出發(fā)來獲得差分格式的優(yōu)化參數(shù)。

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

(44)

加細選擇策略:

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

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

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

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

4 數(shù)值算例

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

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

問題1考慮二維 Helmholtz 方程

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

其邊界條件為

u(x0,y0)=0,(x0,y0)Γ∶=?Ω。

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

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

a=1,k=24;a=3,k=24 和a=3,k=48 時,不同網(wǎng)格剖分下所對應(yīng)的 2 種差分格式的離散L2- 范數(shù)誤差 (表2至表4)表明:與OC6相比,IC6的數(shù)值精度有較大的提高;在參數(shù)a相同的情況下,IC6的數(shù)值精度對波數(shù)k的依賴性弱于OC6的,與前面的理論分析結(jié)果一致。a=1,N=256時,不同波數(shù)k所對應(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的直線近似平行，說明2種差分格式均為六階格式。

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

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

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

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

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

問題2考慮二維 Helmholtz 方程

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

其邊界條件為

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

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

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

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

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

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

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

5 結(jié)論

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