周魯川 吳亭亭

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

1 引 言

考慮如下一維Helmholtz方程

一維Helmholtz方程在地球物理和醫(yī)學(xué)成像等領(lǐng)域有著廣泛應(yīng)用. 因此, 研究一維Helmholtz方程的高效數(shù)值算法具有重要的理論意義和應(yīng)用價(jià)值. 目前, 數(shù)值求解Helmholtz方程的主要方法有:有限元法[1]、差分法[2]等. 有限元法的精確性高, 易于邊界條件的處理, 但其計(jì)算量相對較大, 不利于大規(guī)模的計(jì)算.差分法計(jì)算簡單, 存儲量小, 易于實(shí)現(xiàn), 并且可以通過優(yōu)化差分系數(shù)的方法來提高差分法的數(shù)值精度[3].因此, 本文將重點(diǎn)討論一維Helmholtz方程的優(yōu)化差分法. 具體地, 我們提出了一維Helmholtz方程的帶參數(shù)的差分格式, 證明了該格式為二階格式. 基于極小化數(shù)值頻散的思想, 提出了差分格式的優(yōu)化系數(shù)的整體選取法和加細(xì)選取法. 最后通過數(shù)值試驗(yàn)比較了不同差分格式的精度和數(shù)值頻散情況, 數(shù)值結(jié)果表明了本文所提格式提高了數(shù)值精度, 有效地抑制了數(shù)值頻散.

2 Helmholtz方程的優(yōu)化差分法

本節(jié)將對方程(1)建立帶參數(shù)的差分格式, 并對其進(jìn)行收斂性分析和頻散分析, 最終給出差分格式的優(yōu)化系數(shù)的整體選取法和加細(xì)選取法.

2.1帶參數(shù)的差分格式的建立本小節(jié)建立一維Helmholtz方程的帶參數(shù)的差分格式.

(4)

故令uxx(x)的近似為

(5)

其次, 由(4)式知

u(xi+1)+u(xi-1)≈2u(xi),

則可令k2u(x)的近似為

(6)

其中,c+d=1.

最后, 聯(lián)立(4)和(5)得到方程(1)的差分格式為

(7)

這里fi=f(xi).

2.2收斂性分析本小節(jié)將對差分方程(7)的解進(jìn)行誤差分析, 有如下定理1成立.

定理1 假定u是方程(1)-(3)的解,f是足夠光滑的, 且kh足夠小, 則差分格式(7)得到的差分解U是唯一的, 且有如下估計(jì)

||U-u||≤Ch2.

證為簡化證明過程, 假設(shè)u的邊值條件是Dirichlet邊值條件, 帶有Neumann邊值條件的情況可類似證明. 首先, 將差分格式(7)表示為矩陣形式:

DU=b,

(8)

其中,

(9)

b=(f1,f2,…,fN)T.

因kh是足夠小的, 所以矩陣(9)為嚴(yán)格對角占優(yōu)矩陣. 因此, 方程組(8)得到的差分解U是存在唯一的.

下面, 建立差分解和真解之間的誤差方程. 為此, 令ei=u(xi)-Ui,則有

將上述誤差方程表示為矩陣形式

DE=T,

(10)

其中,E=(e1,e2,…,eN)T,T=(T1,T2,…,TN)T,Ti=O(h2).

接下來, 進(jìn)行誤差分析. 當(dāng)kh→0, 矩陣D趨于如下矩陣

(11)

又矩陣(11)的特征值[4,5]為

相應(yīng)的特征向量為

當(dāng)h是足夠小時(shí), 易知

(12)

(13)

其中,aj是元素ζj的系數(shù).

并由(10)得到:

(DE,E)=(T,E).

(14)

當(dāng)kh→0, 應(yīng)用(12)和(13)以及Cauchy不等式, 有如下結(jié)論:

且又因

(T,E)≤||T||||E||,

根據(jù)(14)并聯(lián)立上述兩個(gè)不等式得到結(jié)論.

從上述定理中, 可以看出差分格式(7)在c+d=1的條件下, 是一個(gè)二階格式.

2.3數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差分析對于大波數(shù)問題,Helmholtz方程的解具有較強(qiáng)的振蕩性. 事實(shí)上, 隨著波數(shù)k的增加, 數(shù)值解的精度在降低—即所謂的“污染效應(yīng)”. “污染”的結(jié)果就導(dǎo)致數(shù)值解的波數(shù)不同于真實(shí)的波數(shù), 這就是所謂的“數(shù)值頻散”[1,6]. 在本小節(jié)中, 我們將分析差分格式(7)的數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差.

接下來, 將Ui:=e-ikxi代入(7)式, 取fi=0, 并利用歐拉公式, 得頻散方程:

2AScos(kh)+A0=0,

(15)

(16)

在如下定理2中, 我們將給出數(shù)值波數(shù)kN和真實(shí)波數(shù)k之間的誤差分析.

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

(17)

證令τ:=kh. 又因方程(16)中P依賴于τ, 故P(τ)=cos(τ). 引入記號:

f1(τ)=2-2cos(τ),f2(τ)=dcos(τ)+c.

(18)

(19)

結(jié)合方程(16)、(18)、(19), 可得到

上述定理表明了kN以二階近似k. 此外, 與k3h2有關(guān)的項(xiàng)稱為污染項(xiàng), 它依賴于波數(shù)k和差分格式(7)中的參數(shù).

2.4優(yōu)化參數(shù)選取策略在本小節(jié)中, 基于極小化數(shù)值頻散的思想, 我們給出差分格式(7)的優(yōu)化參數(shù)的兩種選取策略. 具體地, 極小化數(shù)值頻散就是要極小化數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差. 研究表明, 差分格式的數(shù)值頻散越小, 其數(shù)值精度越高[7,8].

(20)

其次, 為選取合適的參數(shù)c,d來極小化數(shù)值頻散, 令

其中d∈R,G∈IG. 通常, 選取IG:=[Gmin,Gmax]=[4,400][7], 其中Gmin為G的最小值,Gmax為G的最大值, 并且Gmin≥2[7].

通過(20)可以看出, 極小化數(shù)值波數(shù)kN和真實(shí)波數(shù)k之間的誤差等同于極小化范數(shù)||J(d,G)||, 在此總結(jié)為如下選擇策略來選取參數(shù)d.

1) 整體的參數(shù)選取方法.

在給定的條件IG:=[4,400]下, 選取d∈(0,1]以滿足

d=argmin{||J(d,G)||,d∈R}.

(21)

下面,利用最小二乘法極小化目標(biāo)函數(shù)(21)的方法[7]來實(shí)施選擇策略1.若令

J(d,G)=0,

通過整理得方程:

2π2(1-P)d=2π2+G2(P-1).

(22)

故可得到線性方程組

S1d=S2,

(23)

其中

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

c=0.817 7,d=0.182 3.

(24)

為方便引用, 我們稱帶有參數(shù)(24)的差分格式(7)為一維Helmholtz方程的整體三點(diǎn)差分格式(簡稱為global 3p).

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

2) 加細(xì)的參數(shù)選取方法.

估計(jì)區(qū)間IG:=[Gmin,Gmax], 選取d∈(0,1]以滿足

d=argmin{||J(d,G)||,d∈R}.

加細(xì)的選擇方法與整體的相比, 一個(gè)重要的區(qū)別是區(qū)間IG是根據(jù)實(shí)際情況變化的. 在表1中, 列出了多組加細(xì)的優(yōu)化參數(shù). 我們稱帶加細(xì)優(yōu)化參數(shù)(表1)的差分格式(7)為Helmholtz方程的加細(xì)三點(diǎn)差分格式(簡稱為refined 3p). 另外, 稱帶有參數(shù)c=1,d=0的差分格式(7)為傳統(tǒng)的三點(diǎn)差分格式(conventional 3p).

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

3 數(shù)值算例

考慮問題(1)-(3), 取f(x)=-1. 此時(shí)真解為

圖1 三種差分格式數(shù)值相速度曲線

接下來, 我們需要對邊界條件(3)進(jìn)行二階近似.為此, 由Taylor公式得

又因u(1)(xN)=iku(xN)及u(2)(xN)=-k2u(xN)+1, 故邊界條件(3)的近似為

(25)

基于上述問題,我們將比較三種差分格式的數(shù)值精度： 一種是傳統(tǒng)的三點(diǎn)差分格式(conventional 3p), 一種是整體三點(diǎn)差分格式(global 3p), 一種是加細(xì)三點(diǎn)差分格式(refined 3p). 在計(jì)算中, 誤差范數(shù)采用相對C-范數(shù). 其中, C-范數(shù)具體定義為: 對任意復(fù)向量z=[z1,z2,…,zM],

表2與表3對應(yīng)于k=30,200時(shí), 不同的網(wǎng)格點(diǎn)數(shù)(N)下所對應(yīng)的三種差分格式的相對C-范數(shù)誤差. 由表2與表3看出, refined 3p的數(shù)值精度比conventional 3p和global 3p高. 從表3中還可看出refined 3p在N=513時(shí)便可達(dá)到global 3p在N=1 025時(shí)的數(shù)值精度. 這就意味著, 當(dāng)采用refined 3p時(shí), 可用較少的網(wǎng)格點(diǎn)數(shù)得到具有較高精度的數(shù)值解. 進(jìn)一步, 圖2和圖3表明上述三種格式均為二階格式. 因此, 在計(jì)算中refined 3p相比其他兩種格式有著顯著的優(yōu)勢[9,10].

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

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

圖2 k=30時(shí)相對C-范數(shù)下的誤差

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

為進(jìn)一步比較三種差分格式的計(jì)算效力, 圖4(a)與(b)分別給出了三種差分格式在條件kh=1,kh=0.5下的相對C-范數(shù)誤差. 從圖4可以看出, 當(dāng)kh是一個(gè)常數(shù)時(shí), 即每個(gè)波長內(nèi)取的網(wǎng)格節(jié)點(diǎn)固定時(shí),refined 3p的誤差較之其他兩種差分格式的誤差要小. 而且, 隨著波數(shù)k的增大, refined 3p的誤差增大的更緩慢一些. 這些都表明, 我們所提的格式refined 3p有效地提高了數(shù)值精度, 抑制了數(shù)值頻散[11,12].

圖4 k∈[200,1200]時(shí), 相對C-范數(shù)誤差

4 結(jié) 論

在本文中, 首先基于加權(quán)平均的思想建立了帶參數(shù)的差分格式, 分析了格式的收斂性. 接下來, 給出了數(shù)值波數(shù)和真實(shí)波數(shù)之間的誤差, 并基于極小化數(shù)值頻散的思想提出了兩種參數(shù)選取策略. 最后, 數(shù)值算例表明帶加細(xì)化參數(shù)的差分格式提高了數(shù)值精度, 有效地抑制了數(shù)值頻散.