一維帶PML的Helmholtz方程的四階差分法

朱啟華 孫煜然 吳亭亭

( 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，250358，濟(jì)南 )

1 引 言

考慮如下的一維Helmholtz方程:

(1)

對于Helmholtz方程(1), 應(yīng)用PML 技術(shù)將無界域截斷為有界域, 得到帶PML的Helmholtz方程

(2)

其中

其中f0為震源主頻, 一般取f0=25 Hz,LPML為PML 層的厚度,lx為匹配層的點(diǎn)到感興趣區(qū)域與匹配層交界處的距離.a0為一正常數(shù), 可取a0=1.79. 當(dāng)在內(nèi)部區(qū)域時,A≡C≡1, 此時帶PML的Helmholtz方程(2)退化為Helmholtz方程(1).

目前, 數(shù)值求解Helmholtz方程的方法主要有有限元法[4-6]、有限差分法[7-10]等. 有限元法的精度高, 易于復(fù)雜邊界條件和非均勻非結(jié)構(gòu)化網(wǎng)格的處理, 但其計算量相對較大. 有限差分法的計算簡單, 存儲量小, 易于實(shí)現(xiàn), 并且可以通過優(yōu)化差分系數(shù)的方法來提高差分法的數(shù)值精度[1, 7, 9, 10].2013年,Chen等人[1]提出了帶PML的Helmholtz方程的優(yōu)化的9點(diǎn)有限差分法, 該格式為二階格式. 2018年,Dastour等人[2]提出了帶PML的Helmholtz方程的25點(diǎn)優(yōu)化差分格式,該格式為四階格式.文獻(xiàn)[1]與[2]所研究的問題為二維Helmholtz方程. 2019年,Zhou等人[10]提出了一維Helmholtz方程的優(yōu)化差分法, 該格式為二階格式, 沒有考慮吸收邊界條件. 因此, 本文將重點(diǎn)研究一維帶PML 的Helmholtz方程的四階差分法. 首先, 針對一維帶PML 的Helmholtz方程構(gòu)造帶參數(shù)的四階差分格式, 并對差分格式進(jìn)行相容性分析. 其次, 基于極小化數(shù)值頻散的思想, 提出差分格式優(yōu)化系數(shù)的整體選取策略和加細(xì)選取策略. 最后, 數(shù)值算例表明帶加細(xì)參數(shù)的四階差分格式提高了精度, 有效地抑制了數(shù)值頻散, 在計算大波數(shù)問題時具有優(yōu)勢.

2 一維帶PML的Helmholtz方程的四階差分法

本節(jié)建立一維帶PML的Helmholtz方程的四階差分格式, 并對差分格式進(jìn)行相容性分析與頻散分析,提出差分格式的優(yōu)化參數(shù)的整體選擇策略和加細(xì)選擇策略.

2.1一種相容的五點(diǎn)差分格式本小節(jié)建立與一維帶PML的Helmholtz方程相容的五點(diǎn)差分格式.

(3)

(4)

接下來, 建立零階項(xiàng)k2Cu的四階逼近. 借助于Taylor公式, 用五個點(diǎn)來建立k2Cu的四階逼近. 為此, 令

lh(k2Cu)m:=c(k2Cu)m+dlh,0(k2Cu)m,

其中c+d=1. 然后, 用lh(k2Cu)來逼近k2Cu, 即

k2Cu|x=xm≈lh(k2Cu)|x=xm.

(5)

最后, 聯(lián)立方程(4)與(5)得到方程(2)的一種相容的五點(diǎn)差分格式為

(6)

下面的定理1分析了差分格式(6)與帶PML的Helmholtz 方程(2)的相容性.

定理1若d∈(0,1], 且滿足c+d=1，則五點(diǎn)差分格式(6)與帶PML的Helmholtz 方程(2)是相容的, 且為一種四階差分格式.

證假設(shè)xm≤x≤xm+1, 則由Taylor展開式得

(7)

lh(k2Cu)|x=xm=k2Cu-ξ2h4+O(h6)，

(8)

由(7)式與(8)式可知定理1結(jié)論成立.

為便于分析, 下面給出常波數(shù)情況下差分格式(6)在內(nèi)部區(qū)域的具體形式. 由于在內(nèi)部區(qū)域中A≡1,C≡1， 此時, 差分格式(6)可以寫為

LhUm:=A1Um-2+A2Um-1+A3Um+A4Um+1+A5Um+2=gm，

(9)

其中

2.2數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差分析在本小節(jié), 給出差分格式(9)的頻散方程, 并分析數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差.

接下來, 將Ui:=eikx代入差分格式(9), 取gm=0, 并利用歐拉公式得頻散方程為

2A1(2P2-1)+2A2P+A3=0.

(10)

在方程(10)中,將A1,A2,A3中的k用數(shù)值波數(shù)kN替換得[9]

(11)

基于方程(11), 下面的定理2給出了數(shù)值波數(shù)kN與真實(shí)波數(shù)k之間的誤差分析.

定理2對于差分格式(9), 有

(12)

證令τ:=kh.又因?yàn)榉匠蘌依賴于τ, 故P(τ)=cos(τ). 進(jìn)一步, 引入記號:

f1(τ):=cos2(τ)-8cos(τ)+7,

f2(τ):=-2dcos2(τ)+4dcos(τ)+3c+d.

(13)

(14)

結(jié)合方程(11), (13)和(14)可得

上述定理2表明kN是k的四階逼近. 此外, 與k5h4有關(guān)的項(xiàng)稱為污染項(xiàng), 它依賴于波數(shù)k和差分格式(9)中的參數(shù)d. 因此, 可以適當(dāng)?shù)剡x取d來使得kN更好地逼近k, 下一節(jié)將給出d的選取策略.

2.3五點(diǎn)差分格式的優(yōu)化系數(shù)的選擇策略在本小節(jié), 基于極小化數(shù)值頻散的思想, 給出差分格式(9)的優(yōu)化參數(shù)的選擇策略. 具體地, 極小化數(shù)值頻散就是要極小化數(shù)值波數(shù)kN與真實(shí)波數(shù)k之間的誤差. 研究表明, 差分格式(9)的數(shù)值頻散越小, 其精度越高[1,7]. 如果差分格式具有最優(yōu)的收斂階, 并能夠選擇合適的參數(shù)來抑制數(shù)值頻散, 就可以把它看作是帶PML的Helmholtz方程的優(yōu)化差分格式.

(15)

其中N:=P2-8P+7,D:=-2dP2+4dP+3-2d.

1) 策略1: 整體的參數(shù)選取方法.在給定的條件IG=[2.5,400] 下, 選取d∈(0,1]以滿足

(16)

下面, 利用最小二乘法極小化目標(biāo)函數(shù)的方法來實(shí)現(xiàn)整體的參數(shù)選擇策略[1,7]. 若令J(d,G)=0,

通過整理可得

-8π2(P2-2P+1)d=G2(P2-8P+7)-12π2.

其中

該方程的系數(shù)矩陣有r行和1列, 是一個超定方程組. 選擇r=100, 并用最小二乘法去解方程組得到差分格式(9)的一組優(yōu)化參數(shù)為

c=0.898 5，d=0.101 5.

(17)

稱帶有參數(shù)(17)的差分格式(6)或(9)為Helmholtz 方程的整體五點(diǎn)差分格式(簡稱為global 5p).

整體的參數(shù)選取方法只給出一組優(yōu)化參數(shù), 這種做法比較粗糙. 為進(jìn)一步提高差分格式的數(shù)值精度,下面給出加細(xì)的參數(shù)選取方法.

2) 策略2: 加細(xì)的參數(shù)選取方法.第一步：估計區(qū)間IG：=[Gmin,Gmax]；第二步：選取d∈(0,1]以滿足(16)式.

加細(xì)的選取方法與整體的選取方法相比, 一個重要的區(qū)別在于區(qū)間IG是根據(jù)實(shí)際情況變化的. 在表1 中, 列出了多組加細(xì)的優(yōu)化參數(shù). 稱帶加細(xì)優(yōu)化參數(shù)(以表1中的參數(shù)為加細(xì)優(yōu)化參數(shù))的差分格式(6)或(9)為Helmholtz方程的加細(xì)五點(diǎn)差分格式(簡稱為refined5p). 另外, 稱帶有參數(shù)c=1,d=0的差分格式(6)或(9)為Helmholtz方程的五點(diǎn)差分格式(簡稱為5pwithc=1,d=0).

表1 加細(xì)優(yōu)化參數(shù)

為比較5pwithc=1,d=0,global5p和refined5p三種差分格式的數(shù)值頻散, 圖1和圖2分別給出了這三種差分格式的標(biāo)準(zhǔn)化的數(shù)值相速度曲線圖和標(biāo)準(zhǔn)化的數(shù)值群速度曲線圖.global5p的數(shù)值頻散相比5pwithc=1,d=0的數(shù)值頻散有較大改善, refined5p的數(shù)值頻散最小. 在數(shù)值算例中會進(jìn)行更詳細(xì)地比較.

圖1 三種差分格式的標(biāo)準(zhǔn)化的數(shù)值相速度曲線圖

圖2 三種差分格式的標(biāo)準(zhǔn)化的數(shù)值群速度曲線圖

3 數(shù)值算例

為了驗(yàn)證本文所構(gòu)造格式的有效性, 在本節(jié)中通過兩個算例來進(jìn)行數(shù)值驗(yàn)證.

算例1考慮如下一維Helmholtz方程：

該問題的真解為

接下來, 需要對邊界條件(20)進(jìn)行四階逼近. 為此, 由Taylor 公式展開得

將上述三個式子整理得

又由(18)式和(20)式可知

u(1)(xN)=iku(xN)，u(2)(xN)=-k2u(xN)+1.

故可得到邊界條件(20)的近似為

對u(x1),u(xN-1)通過內(nèi)插法進(jìn)行五階逼近,可得

表2與表3分別對應(yīng)于k=800和k=1 000時, 不同的網(wǎng)格剖分所對應(yīng)的三種差分格式的相對C-范數(shù)誤差. 從表中可以看出, 對于相同的步長而言,refined5p的誤差比global5p的誤差和5pwithc=1,d=0的誤差都要小,這說明了在這三種差分格式中refined5p的數(shù)值精度最高. 進(jìn)一步, 圖3和圖4表明了上述三種格式均為四階格式. 因此,refined5p相比其他兩種格式有著更明顯的優(yōu)勢, 尤其在計算大波數(shù)問題時.

表2 k=800時相對C-范數(shù)下的誤差

表3 k=1 000時相對C-范數(shù)下的誤差

圖5分別給出了三種差分格式在kh=0.25,kh=0.5下的相對C-范數(shù)誤差. 當(dāng)kh固定時, 即表示每個波長內(nèi)取的網(wǎng)格節(jié)點(diǎn)固定,refined5p比其他兩種差分格式的誤差要小. 不難發(fā)現(xiàn), 當(dāng)波數(shù)k由800增加到1 200時, 三種差分格式的誤差也在增大. 但相對而言,refined5p的誤差增大的比較緩一些, 從而更好地說明了refined5p的數(shù)值頻散比5pwithc=1,d=0和global5p數(shù)值頻散要小. 圖5表明了本文所提出的差分格式refined5p可以更好地提高數(shù)值精度, 抑制數(shù)值頻散.

圖3 k=800時相對C-范數(shù)下的誤差

圖4 k=1 000時相對C-范數(shù)下的誤差

算例2考慮無界區(qū)域上的一維Helmholtz方程

-u″-k2u=δ(x-x0),x0=0,x∈(-,+)，

將無界域截斷為有界域, 可借助于PML技術(shù). 例如下面的方程

(a)kh=0.25

(b)kh=0.5圖5 當(dāng)k∈[800,1 200]時, 相對C-范數(shù)誤差

假設(shè)感興趣的區(qū)域?yàn)閇-5,5], 可以選擇優(yōu)化系數(shù)(加細(xì)參數(shù)選取策略)并用差分格式(6)來求解問題(21)-(22). 在計算中, 采用相對L-范數(shù)和相對L2-范數(shù)來衡量誤差. 具體地, 若設(shè)u為真解,U為數(shù)值解, 則

表4和表5分別給出了k=50,k=100, 取不同的步長時的相對L-誤差與相對L2-誤差以及收斂階, 從表中可以看出, 當(dāng)步長固定時, 相應(yīng)的誤差隨著k的增大而增大. 進(jìn)一步, 圖6更好地說明了在不同相對范數(shù)的意義下, 差分格式均達(dá)到四階.

表4 k=50時 相對L-范數(shù)和相對L2-范數(shù)及誤差階

表4 k=50時 相對L-范數(shù)和相對L2-范數(shù)及誤差階

表5 k=100時相對L-范數(shù)和相對L2-范數(shù)及誤差階

表5 k=100時相對L-范數(shù)和相對L2-范數(shù)及誤差階

(a) 相對L-誤差 (b) 相對 L2-誤差圖6 不同相對范數(shù)意義下的差分格式

(a) k=1的精確解和取時所得數(shù)值解實(shí)部曲線

(b) k=1的精確解和取時所得數(shù)值解的虛部曲線圖7 k=1時精確解與數(shù)值解曲線

(a) k=3的精確解和取時所得數(shù)值解實(shí)部曲線圖8 k=3時精確解與數(shù)值解曲線

(b)k=3的精確解和取時所得數(shù)值解的虛部曲線圖8 k=3時精確解與數(shù)值解曲線

4 結(jié) 語

在本文中,給出了一維帶PML的Helmholtz方程的五點(diǎn)差分格式, 該格式為四階格式. 首先, 建立了帶參數(shù)的五點(diǎn)差分格式, 并給出該格式的相容性分析. 接下來, 分析了數(shù)值波數(shù)和真實(shí)波數(shù)之間的誤差. 然后, 基于極小化數(shù)值頻散的思想, 提出了差分系數(shù)的整體優(yōu)化選取策略和加細(xì)優(yōu)化選取策略. 最后,通過數(shù)值算例驗(yàn)證了所提出的差分格式能夠提高數(shù)值精度, 有效地抑制數(shù)值頻散, 在計算大波數(shù)問題時具有顯著優(yōu)勢.