孫煜然 朱啟華 吳亭亭

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

1 引 言

考慮如下一維Helmholtz方程:

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

本節(jié)首先對(duì)Helmholtz方程(1)建立帶參數(shù)的緊致差分格式. 其次, 對(duì)其進(jìn)行收斂性分析和頻散分析. 最后,在極小化數(shù)值頻散的基礎(chǔ)上, 提出差分系數(shù)的整體選取策略和加細(xì)選取策略.

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

首先, 對(duì)區(qū)域[0,1]進(jìn)行N等分,N為一正整數(shù). 空間步長(zhǎng)為h:=1/N, 記xm:=mh,m=0,1,2,…,N, 并記Um為u(xm)的差分近似. 其次, 用二階中心差分Dxxu近似uxx. 由Taylor展開可以得到:

(4)

為建立緊致格式, 需對(duì)(4)中的uxxxx進(jìn)行改寫. 為此, 對(duì)Helmholtz方程(1)關(guān)于x求二階導(dǎo)數(shù)得

uxxxx=fxx-k2uxx.

將上式代入(4), 可以得到uxx的四階逼近格式如下:

(5)

接下來, 建立零階項(xiàng)k2u的四階逼近. 令

(6)

其中c+d=1, 且lh(k2u)是k2u的二階逼近[8]：

(7)

將(5)式代入(7)式,得到k2u的四階逼近為

(8)

最后, 聯(lián)立(5)式和(8)式得到Helmholtz方程(1)的四階緊致差分格式為

(9)

其中

fm:=f(xm).

類似于文獻(xiàn)[8]中定理1的證法, 可得差分方程(9)的解的收斂性分析, 即有如下定理1成立.

定理1假定u是方程(1)-(3)的解, 函數(shù)f是足夠光滑的, 且kh足夠小, 則差分方程(9)的解U存在唯一, 且有如下估計(jì)

||U-u||≤Ch4.

上述定理1表明差分格式(9)為四階格式. 易知, 差分格式(9)為三點(diǎn)差分格式. 因此, 差分格式(9)為帶參數(shù)的四階緊致差分格式.

2.2數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差分析在本小節(jié)中,將分析差分格式(9)的數(shù)值波數(shù)與真實(shí)波數(shù)之間的誤差. 首先, 通過經(jīng)典的頻散分析得到緊致差分格式(9)的頻散方程, 該分析基于波數(shù)k為常數(shù)和源項(xiàng)f為零的基礎(chǔ)上. 為便于分析,將緊致差分格式(9)改寫為

As(Um-1+Um+1)+A0Um=B0fm+Bs(fm+1+fm-1),

(10)

其中

接下來, 將Ui:=eikxi代入差分方程(10), 取fi=0, 并利用歐拉公式, 得頻散方程:

2AsP+A0=0,

(11)

其中P:=cos(kh). 用數(shù)值波數(shù)kN替代方程(11)式AS與A0中的k得到關(guān)于數(shù)值波數(shù)與真實(shí)波數(shù)的方程[4]. 基于該方程,下面定理2給出數(shù)值波數(shù)kN和真實(shí)波數(shù)k之間的誤差.

定理2對(duì)于緊致差分格式(9), 有

(12)

證令τ:=kh. 又方程(11)中P依賴于τ, 即P(τ)=cos (τ). 定義X：=(kNh)2,由(11)式得

a(τ)X2+b(τ)X+c(τ)=0,

(13)

(14)

(15)

由于kN是對(duì)k的數(shù)值逼近, 故當(dāng)τ→0時(shí), 應(yīng)有X→0. 結(jié)合方程(13)、(14)和(15),得方程(13)有唯一解

從而, 數(shù)值波數(shù)kN為

(16)

接下來, 將方程(14)和(15)代入方程(16)得

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

(17)

一般的, 可選取IG:=[2.5,400][10]. 下面提出整體的參數(shù)選取策略來選擇d.

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

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

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

(18)

下面, 利用奇異值分解法極小化目標(biāo)函數(shù)(18)的方法[9,10]來實(shí)施選擇策略1.若令

J(d,G)=0,

通過整理得方程

(19)

S1d=S2,

(20)

其中

方程組(20)的系數(shù)矩陣有r行和1列, 為超定方程組. 選擇r=100, 并用奇異值分解法求解方程組(20),可得緊致差分格式(9)的一組優(yōu)化參數(shù)為.

c=0.902 5,d=0.097 5.

(21)

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

整體的參數(shù)選取方法只給出一組優(yōu)化參數(shù), 這種做法比較粗糙. 因?yàn)閷?duì)不同的頻率和速度等問題,均采用同一組優(yōu)化參數(shù)進(jìn)行計(jì)算,仍可能會(huì)產(chǎn)生較大的數(shù)值頻散. 為進(jìn)一步提高差分格式的數(shù)值精度,下面給出加細(xì)的參數(shù)選取方法.

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

步驟1：估計(jì)區(qū)間IG：=[Gmin,Gmax];

步驟2：選取d∈(0,1] 以滿足

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

同樣可通過奇異值分解法求解相應(yīng)的線性系統(tǒng).

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

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

標(biāo)準(zhǔn)化的數(shù)值相速度是衡量數(shù)值頻散的重要標(biāo)準(zhǔn). 圖1給出了compact 3p, global 3p和 refined 3p的標(biāo)準(zhǔn)化的數(shù)值相速度曲線. 我們發(fā)現(xiàn) global 3p相比 compact 3p 有較大的改善, refined 3p的數(shù)值頻散最小. 在數(shù)值算例部分將給出更加詳細(xì)的比較.

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

3 數(shù)值算例

在本節(jié),通過兩個(gè)算例來展示本文構(gòu)造的差分格式的有效性. 具體地,將比較compact 3p,global 3p,refined 3p這三種差分格式的數(shù)值精度. 在計(jì)算中，誤差范數(shù)采用相對(duì)C-范數(shù).其中, C-范數(shù)具體定義為: 對(duì)任意復(fù)向量z=[z1,z2,…,zM],其范數(shù)為

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

(22)

其中xN+1為虛擬點(diǎn). 基于方程(1), 有uxxx=-k2ux， 將其代入方程(22), 得到ux的四階逼近為

綜上,得到建立邊界條件(3)的四階近似為

(23)

在計(jì)算過程中,把虛擬點(diǎn)xN+1作為未知量處理, 建立個(gè)數(shù)與未知量個(gè)數(shù)相同的方程. 為此, 對(duì)于節(jié)點(diǎn)xm，m=1,2,…,N,應(yīng)用內(nèi)部差分格式(9)得到方程, 在虛擬點(diǎn)xN+1處應(yīng)用(23)式建立方程.

表2、表3為算例1對(duì)應(yīng)于k=200、500時(shí), 不同的網(wǎng)格點(diǎn)數(shù)N所對(duì)應(yīng)的三種差分格式的相對(duì)C-范數(shù)誤差. 由表格看出, 當(dāng)步長(zhǎng)相同時(shí), refined 3p 的誤差比compact 3p和 global 3p 小得多.在這三種差分格式中, refined 3p的數(shù)值精度最高. 對(duì)于大波數(shù)問題, 為達(dá)到相同精度, refined 3p僅需要compact 3p四分之一的節(jié)點(diǎn)數(shù).圖2、圖3進(jìn)一步說明這三種差分格式均為四階格式, 與理論分析一致.

表2 算例1,k=200時(shí)相對(duì)C-范數(shù)下的誤差

表3 算例1,k=500時(shí)相對(duì)C-范數(shù)下的誤差

圖2 算例1,k=200時(shí)相對(duì)C-范數(shù)下的誤差

圖3 算例1,k=500時(shí)相對(duì)C-范數(shù)下的誤差

在圖4中,進(jìn)一步比較了三種差分格式的計(jì)算效力. 圖4(a),圖4(b)分別為算例1對(duì)應(yīng)于三種差分格式在條件kh=1,kh=0.5下的相對(duì)C-范數(shù)誤差.kh為一個(gè)常數(shù), 表示每個(gè)波長(zhǎng)內(nèi)取的網(wǎng)格節(jié)點(diǎn)固定. 在圖4中, 波數(shù)k由200增加到1 200,我們看到refined 3p的精度始終比compact 3p和 global 3p的精度高. 而且, 隨著波數(shù)k的增大, compact 3p和 global 3p的誤差越來越大, 而refined 3p的誤差則維持在穩(wěn)定的狀態(tài), 說明refined 3p可以較好的減弱大波數(shù)的污染.

圖4(a) 算例1,當(dāng)k∈[200,1 200]時(shí)相對(duì)C-范數(shù)誤差

圖4(b) 算例1,當(dāng)k∈[200,1 200]時(shí)相對(duì)C-范數(shù)誤差

3.2算例2考慮如下Helmholtz方程：

該問題的精確解為u(x)=eikx. 接下來,用該例子來比較compact 3p, global 3p和refined 3p三種差分格式的數(shù)值精度. 在計(jì)算中, 邊界的處理方式與算例1中對(duì)邊界的處理方式相同.

在表4與表5中, 給出了對(duì)于不同節(jié)點(diǎn)數(shù)N下k=200、500時(shí)的三種差分格式的誤差. 容易發(fā)現(xiàn), refined 3p的精度最高. 為達(dá)到某一精度, refined 3p僅需要compact 3p一半的節(jié)點(diǎn)數(shù)即可. 進(jìn)一步, 圖5與圖6說明這三種差分格式均為四階格式, 與理論分析一致.

表4 算例2,k=200時(shí)相對(duì)C-范數(shù)下的誤差

表5 算例2,k=500時(shí)相對(duì)C-范數(shù)下的誤差

在圖7中,對(duì)compact 3p, global 3p和refined 3p 這三種差分格式進(jìn)行了更深入的比較. 圖7(a),圖7(b)分別給出了三種差分格式在kh=1,kh=0.5時(shí)的相對(duì)C-范數(shù)下的誤差. 在這兩個(gè)圖中, 波數(shù)k由200增加到1 200, 我們看到refined 3p比其他兩種格式的精度要高. 并且, 隨著波數(shù)k的增大, refined 3p的數(shù)值精度的變化相比compact 3p, global 3p的數(shù)值精度的變化要緩慢的多, 這說明refined 3p在計(jì)算大波數(shù)問題時(shí)有優(yōu)勢(shì).

綜上,compact 3p, global 3p, refined 3p這三種差分格式均為四階格式. 在這三種格式中, refined 3p的數(shù)值精度最高, 并可有效抑制數(shù)值頻散, 在計(jì)算大波數(shù)問題時(shí)有優(yōu)勢(shì).

圖5 算例2,k=200時(shí)相對(duì)C-范數(shù)下的誤差

圖6 算例2,k=500時(shí)相對(duì)C-范數(shù)下的誤差

(a)

(b)

圖7 算例2,當(dāng)k∈[200,1 200]時(shí), 相對(duì)C-范數(shù)誤差

4 結(jié) 語(yǔ)

本文給出了一維Helmholtz方程的四階優(yōu)化緊致差分格式. 首先, 建立了帶參數(shù)的四階緊致差分格式并給出收斂性分析. 接下來, 給出了數(shù)值波數(shù)和真實(shí)波數(shù)之間的誤差分析, 并基于極小化數(shù)值頻散的思想提出了兩種參數(shù)選取策略, 同時(shí)畫出數(shù)值相速度曲線, 比較說明加細(xì)的參數(shù)選取策略效果更好. 最后,通過數(shù)值算例驗(yàn)證了所提出的差分格式提高了數(shù)值精度,并且有效地抑制了數(shù)值頻散.