求解Helmholtz方程的無網格重心插值配點法

楊苗苗，葛永斌

(寧夏大學數學統計學院，寧夏銀川750021)

1.引言

Helmholtz方程是一個橢圓型偏微分方程，代表了波動方程與時間無關的解.這個方程模擬了各種各樣的物理現象.其中包括：振動分析、水波傳播、聲波散射、電磁散射、以及雷達散射等問題.隨著Helmholtz方程被廣泛地應用在工程和科學領域的研究中，該方程本身的復雜特性給數值計算帶來了巨大困難，特別是對于變波數和大波數問題的有效數值計算，雖然其數值解法有很多，但都需要進一步提高計算精度和計算效率，以適應實際復雜大規模問題的求解.所以，對Helmholtz方程的高效精確數值求解方法的研究具有重要的理論價值和實際意義.

無網格法[1]是近年來興起的一種新型、高效的數值計算方法.它基于點的近似思想來構造近似函數，適用于分析各類具有大梯度、奇異性等特殊性質的應用問題.目前發展的無網格方法包括無單元Galerkin法[2]和無網格配點法[3]等.對于無網格配點法，重心插值配點法就是其中的一種.重心插值配點法包括重心Lagrange插值法和重心有理插值法，是指未知函數的近似函數在插值節點處的重心型插值，是一種依賴于偏微分方程的強形式的配點法.該數值方法具有易于理解，思想結構簡單，計算精度高，便于編程的優點.近年來，劉婷和馬文濤[4]采用重心Lagrange插值配點法求解了二維電報方程，得到了一種高效的數值方法.王彩珍等[5]利用重心插值配點法，構造了二維非線性橢圓型方程所對應的離散方法，得到了一種計算效率更高的數值解法.

對于Helmholtz方程的求解，近年來研究者們提出了很多無網格方法，如李鵬等[6]直接利用最小二乘法建立系統的變分公式，推導了大波數Helmholtz問題的加權無網格公式.何锃等[7]通過使用加權正交基函數構造了改進的移動最小二乘法，并結合拉格朗日乘子法離散Dirichlet邊界條件，推導了求解該方程的無網格加權最小二乘法.陳黃鑫和邱偉峰[8]提出了一種求解大波數Helmholtz方程的一階系統最小二乘法，并利用一個非平凡分解來解決對偶的一階系統最小二乘法問題.楊子樂等[9]利用無網格節點MIP法的d適應性和移動最小二乘法來構造近似函數，給出了求解Helmholtz方程的兩種計算方法.李瑩等[10]采用無網格局部徑向基點插值法，將徑向基函數耦合多項式基函數作為近似函數，克服了無網格法中場函數在計算中的冗余性.Assari等[11]將局部支撐薄板樣條的離散配點法近似化為一類具有自由形狀參數的徑向基函數法，用一種特殊的非等距Gauss Legendre求積法來計算出現的對數型奇異積分，發展了一種求解具有對數積分方程的無網格法.李美香等[12]將三角函數作為基函數，并在局部支持域內構造了一種形函數，研究了帶有邊界層問題和波傳播問題的Helmholtz方程.張偉和丁睿[13]通過引入全局徑向基函數和正定緊支徑向基函數構造了求解Helmholtz方程的Galerkin型的無網格方法，使得求解方法有效且具有高精度.屈文鎮等[14]發展了一種新的局部基本解方法，與傳統的基本解法相比，該方法根據節點將計算域劃分為若干重疊子域后得到了稀疏的帶狀矩陣，可以求解大波數Helmholtz問題.陳林沖和李小林[15]基于無網格和無邊界分析，利用Burton-Miller公式提出了一種無邊界方法，該無網格方法對所有的波數都能求出唯一解，并且只涉及到邊界節點，還可以處理特大波數下的Helmholtz問題.

本文針對如下一維Helmholtz方程:

其中[a，b]是求解域，u(x)是關于x的未知函數.k(x)表示波數，可以是常數，也可以是關于x的函數.f(x)是源項.定義Dirichlet邊界條件為

借助Chebyshev插值節點和三種測試節點，運用重心Lagrange插值基函數和重心有理插值基函數推導了求解該類方程的兩種無網格配點法.該重心插值方法理論簡單，需要的插值節點少，最后得到的離散矩陣處理方便，不僅可以計算大波數問題，還可以計算變波數問題.具有推導過程簡單，計算結果精度高和穩定性好等的優點.主要內容如下:第二節推導一維重心Lagrange插值和重心有理插值的計算公式和一些簡化矩陣的理論知識.第三節推導計算一維Helmholtz方程的兩種插值配點法及其離散矩陣.第四節給出一些測試問題的計算結果，并與一些文[19-21]中結果進行對比，最后給出本文的結論.

2.重心插值公式

在區間[a，b]上插入n+1個節點xi，則有a=x0＜x1＜x2＜···＜xn=b，設未知函數u(x)在節點處的函數值ui=u(xi)，i=0，1，···，n，則Lagrange多項式插值公式可表示為:

其中Li(x)為Lagrange插值基函數，具體形式為:

令

定義重心Lagrange插值函數的權重為:

此時，插值基函數(2.2)式可變為:

將(2.5)代入(2.1)式，得:

則上式被稱作改進的Lagrange插值公式.相比Lagrange插值公式，上式在插值節點的數目增加時，計算量減少很多，由O(n2)變為O(n).

利用(2.6)式插值常數1，即給定一組節點為(xi，1)，i=0，1，···，n，得:.

化簡后可得:

將(2.7)代入(2.6)式中可得:

則上式被稱為重心Lagrange插值公式[16]，相比式(2.6)，式(2.8)在保持相同計算量O(n)的同時，還能有效地克服Runge現象.

重心有理插值和重心Lagrange插值相類似，重心有理插值即是將Lagrange函數換為有理函數進行插值.下面我們構造重心有理插值基函數.

首先，選擇一個正整數d，且0≤d≤n.對于節點xi，i=0，1，···，n-d，選擇d+1個節點(xi，ui)，(xi+1，ui+1)，···，(xi+d，ui+d)后，我們利用重心Lagrange插值多項式的形式，構造

然后，根據pi(x)構造一個權函數

最后，利用pi(x)和λi(x)構造出重心有理插值函數

利用常數1的Lagrange插值公式，有

由此可得

將(2.11)代入(2.9)中可得重心有理插值公式：

則上式被稱為重心有理插值公式[17].需要說明的一點是，文[17]中已經證明式(2.13)構造的有理插值的誤差階為O(hd+1).其中h是等距節點網格步長.因此只要選擇合適的參數d，式(2.13)就可達到更高的插值精度，從而得到更小的誤差.

從重心Lagrange插值權函數的表達式(2.4)和重心有理插值權函數的表達式(2.10)能夠看出，插值權的選取跟節點的分布有關.因此我們考慮了三種不同的節點來進行函數插值和數據測試，它們分別是隨機節點、等距節點和Chebyshev節點.由于在進行具體的程序計算時，隨機節點和等距節點可直接調用命令，因此，這部分只給出Chebyshev節點的計算公式:

第一類Chebyshev節點:(Chebyshev I)

第二類Chebyshev節點:(Chebyshev II)

需要說明的是，對于兩類Chebyshev節點的公式，其定義區間為[-1，1]，當針對一般的區間[a，b]時，利用區間的坐標變換公式進行轉化即可.

微分矩陣最初是在研究Chebyshev擬譜方法[18]時提出的，重心插值微分矩陣能夠直接獲得Helmholtz方程中未知函數在計算節點處的導函數，因此，微分矩陣是重心插值配點法求解里非常重要的一部分.

由式(2.8)，函數的重心Lagrange插值公式可表示為：

通過對比我們可以看出，兩種重心插值公式的根本區別是權函數的不同.不難發現，重心有理插值除了必要的計算外，還需選取參數d.除此之外，兩種重心插值法在形式上是一致的.因此，利用(2.16)和(2.17)式，函數在節點xj(j=1，2，···，n)處的m階導數可以分別表示為:

寫成矩陣的形式為

3.Helmholtz方程的離散

方法I重心Lagrange插值方法

將重心Lagrange插值公式(2.16)代入到式(1.1)中，可得:

設xj，j=1，2，···，n，為給定的插值節點，令上式方程在所有xj上都成立，則有

注意到Li(xj)=δij，利用重心插值微分矩陣，則方程(3.2)可寫為矩陣形式：

進一步可寫為：

將式(2.16)代入到式(1.2)中，可得:

其中，分別表示n階單位矩陣的第一行和第n行，g1，g2分別為具體的值.聯立求解方程(3.4)和(3.5)，即可求得在重心Lagrange插值公式下插值節點上的函數值ui，回代入式(2.16)中并利用測試節點進行計算可得本文方法I的數值解u(x).

方法II重心有理插值方法

同理，將重心有理插值公式(2.17)代入到式(1.1)中，可得:

4.數值算例

在本節中，我們將兩種重心插值配點法應用于幾個測試問題中以證明本文方法的有效性和準確性.我們給出了包括大波數和變波數Helmholtz方程的數值算例，采用本文方法進行計算并與文[19-21]中所得的數值結果進行比較.其中L∞范數誤差的定義如下給出：

上式中ui和uei分別代表數值解和精確解.Nt是檢測點個數.

問題1[19－20]:

考慮如下非齊次大波數Helmholtz方程

其中邊界條件為:

該問題的精確解為:

首先，由于方法II在進行計算中需要優先選擇參數d，并且0≤d≤N.其中N是插值點個數.因此表1給出當Nt=30，k=10時不同d和N下問題1中方法II的L∞范數誤差.其中Nt是檢測點個數，并且插值節點的類型為第二類切比雪夫點(Chebyshev II)，測試節點也為第二類切比雪夫點.由數值結果我們可以看出，隨著插值點數N的增大，當d越接近于N，方法II的計算結果精度越高，誤差越小.因此在本文剩余算例關于方法II的計算中，我們統一選取d=N-1.這樣既保證了方法II具有更高的精度，又得到了更小的計算誤差.

表1 當Nt=30，k=10時不同d和N下問題1中本文方法II的L∞范數誤差

其次，表2給出了當Nt=N，k=10時四種插值節點下問題1的L∞范數誤差.由數值結果可以得出對于四種插值節點:隨機節點、等距節點、Chebyshev I和Chebyshev II節點來說，隨機插值節點的數值結果很差;等距節點的結果一般;Chebyshev I節點的結果遠不如Chebyshev II節點，并且Chebyshev II節點的精度最高、穩定性最好、計算誤差最小.因此對于插值節點的選擇，我們對于本文剩余算例中的表格和圖都默認選擇Chebyshev II節點，簡記為Chebyshev節點.此外對于兩種插值方法:重心Lagrange插值法和重心有理插值法來說，所得的計算結果相差不大，即精度一樣，穩定性一樣，誤差量級一樣并且L∞范數誤差也近似相等.

表2 當Nt=N，k=10時四種插值節點下問題1的范數誤差

然后，表3給出了當Nt=50，k=10時三種測試節點下的L∞范數誤差.由數值結果我們可以得出對于三種測試節點：隨機節點、等距節點和Chebyshev節點來說，所得的誤差量級一樣，并且誤差也相差不大.此外對于兩種插值方法來說，所得的計算結果也相差不大.

表3 當Nt=50，k=10時三種測試節點下問題1的L∞范數誤差

接著，表4給出了當N=24，k=10時三種測試節點下的L∞范數誤差.由數值結果我們可以得出兩種重心插值方法隨著檢測點數Nt的變化，所得的誤差量級不變，并且L∞范數誤差也相差不大.也就是說檢測點的個數對于本文所提方法的誤差影響不大，因此在剩余算例的計算中，我們可以取任意合理的檢測點個數.

表4 當N=24，k=10時三種測試節點下問題1的L∞范數誤差

最后，表5給出了當Nt=200，k=1時本文兩種無網格配點法與文[19]和文[20]中方法計算L∞范數誤差的比較.結果表明當N=20時，本文所提兩種方法的誤差量級均為O(10－15)，而文[19-20]中的誤差量級僅僅只達到O(10－5).而當N=160時本文所提方法的誤差量級均為O(10－13)，而文[19-20]中的誤差量級分別為O(10－8)和O(10－12).另外，需要說明是，盡管本文所提方法的計算精度會隨著插值點數N的增加而減小，但本文方法的誤差量級仍然要高于文[19-20]中所提方法.此外，與有限差分法相比，無網格重心插值配點法并不意味著插值節點數越多，所得的計算結果會越好，有時候我們僅僅可能只用很少的插值節點便可得到更好的數值結果.表6給出了當Nt=200，k=1000時本文兩種無網格配點法與文[20]中方法L∞范數誤差的比較.結果表明本文所提兩種無網格方法的計算結果要遠遠優于文[20]中所提方法.因此本文方法更適合求解大波數的Helmholtz問題.

表5 當Nt=200，k=1時三種測試節點下問題1的L∞范數誤差

表6 當Nt=200，k=1000時不同測試節點下問題1的L∞范數誤差

問題1在Nt=24，Nt=N，k=10時隨機點(a)，等距點(b)，Chebyshev I點(c)及Chebyshev II點(d)的四種插值節點圖可參見圖1.由圖可得：隨機點不規則的分布在區間內，等距點均勻分布在區間內，兩類Chebyshev點在端點處分布點數較多，中間分布點數相對較少.隨后，在Nt=200，k=1000，N=160時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖參見圖2.由圖可得三種測試點下的精確解和數值解吻合的很好.

圖1 問題1在N=24，Nt=N，k=10時隨機點(a)，等距點(b)，Chebyshev I點(c)及Chebyshev II點(d)的插值節點圖

圖2 問題1在Nt=200，k=1000，N=160時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖

問題2[20]:

考慮如下非齊次變波數Helmholtz方程

其中邊界條件為:

該問題的精確解為:

表7給出了當Nt=100，k(x)=時本文兩種無網格配點法與文[20]中方法L∞范數誤差的比較.結果表明當N=8時本文所提兩種方法的誤差量級均為O(10－6)，而文[20]中的誤差量級為O(10－4).而當N=64時本文所提方法的誤差量級為O(10－14)或O(10－15)，而文[20]中的誤差量級為O(10－10).表8給出了當Nt=200，k(x)=時本文兩種無網格配點法與文[20]中方法L∞范數誤差的比較.結果表明本文所提兩種無網格方法的計算結果要優于文[20]中所提方法.因此本文方法更適合求解變波數的Helmholtz問題.

表7 當Nt=100，k(x)=時不同測試節點下問題2的L∞范數誤差

表7 當Nt=100，k(x)=時不同測試節點下問題2的L∞范數誤差

N文[20]隨機節點等距節點Chebyshev節點方法I方法II方法I方法II方法I方法II 81.610(-4)1.175(-6)1.195(-6)1.194(-6)1.194(-6)1.198(-6)1.198(-6)161.695(-6)2.887(-15)3.997(-15)2.887(-15)3.997(-15)2.776(-15)3.997(-15)321.335(-8)6.772(-15)4.996(-15)6.772(-15)4.996(-15)6.772(-15)4.996(-15)64 1.011(-10)4.219(-14)8.573(-15)4.208(-14)8.704(-15)4.186(-14)8.704(-15)

表8 當Nt=200，k(x)=時不同測試節點下問題2的L∞范數誤差

表8 當Nt=200，k(x)=時不同測試節點下問題2的L∞范數誤差

N文[20]隨機節點等距節點Chebyshev節點方法I方法II方法I方法II方法I方法II 85.563(-4)2.046(-6)2.047(-6)2.048(-6)2.048(-6)2.048(-6)2.048(-6)166.045(-6)1.221(-14)1.366(-14)1.221(-14)1.399(-14)1.221(-14)1.399(-14)324.377(-8)2.676(-14)7.105(-15)2.676(-14)7.105(-15)2.676(-14)7.105(-15)64 3.028(-10)5.596(-14)1.993(-13)5.584(-14)1.993(-13)5.584(-14)1.988(-13)

問題2在Nt=200，k(x)=N=64時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖參見圖3.由圖可得三種測試點下的精確解和數值解吻合得非常好.

圖3 問題2在Nt=200，k(x)==64時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖

問題3[21]:

考慮如下齊次Helmholtz方程

其中邊界條件為:

該問題的精確解為:

表9給出了當Nt=100，b=0.64時不同測試節點及不同波數k=0.5N下本文兩種無網格配點法的L∞范數誤差.需要說明的是，文[21]中只給出了一些數值解和精確解的吻合圖而沒有給出相應的數據，因此表9只給了本文所提兩種無網格配點法的計算結果而沒有和文[21]的對比數據.計算結果表明，即使k隨著N的變化而變化，但是本文的兩種無網格配點法仍然精度很高，并且誤差值很小，也就是說本文重心無網格方法更逼近于精確解.

表9 當Nt=100，b=0.64時不同測試節點及不同波數k=0.5N下問題3的L∞范數誤差

問題3在N=40，k=0.5N，Nt=100，b=0.64時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖參見圖4.由圖可得三種測試點下的精確解和數值解吻合的非常好.

圖4 問題3在N=40，k=0.5N，Nt=100，b=0.64時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖

問題4:

考慮如下齊次Helmholtz方程

其中邊界條件為:

該問題的精確解為:

表10給出了當Nt=100，k=24時不同測試節點下本文兩種無網格配點法的L∞范數誤差.由數值結果可以看出，該區間上本文所提兩種無網格方法的誤差量級最高為O(10－6)，盡管誤差量級與前面幾個算例比起來不高，但是穩定性仍然很好.

表10 當Nt=100，k=24時不同測試節點下問題4的L∞范數誤差

問題4在Nt=100，k=24，N=56時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖參見圖5.由圖可得三種測試點下的精確解和數值解吻合的非常好.

圖5 問題4在Nt=100，k=24，N=56時方法I與方法II在隨機點(a)，等距點(b)和Chebyshev測試點(c)下精確解與數值解圖

5.結論

可以得出以下結論:

1)兩種重心插值配點法的計算精度相同，誤差也幾乎相同，因此對于一維Helmholtz方程來說，兩種方法的結果相差不大.

2)在插值節點的選取上，數值結果表明:隨機節點很差，等距節點一般，兩種Chebyshev節點結果很好，并且Chebyshev II要優于Chebyshev I.在測試節點的選取上，三種測試節點的數值結果幾乎沒有差異.此外，隨著測試節點數量的變化，本文方法的誤差幾乎沒有變化.也就是說，測試節點的數量對我們的數值結果影響不大.

3)與文[19-21]中所提方法相比，本文所提兩種無網格配點法使用較少的節點即可獲得高精度的數值結果，并保持良好的穩定性.并且數值實驗研究結果進一步表明，本文方法具有所需配點數少、精度高、誤差小和效率高等優點.此外，本文方法不僅可以計算大波數問題，還可以計算變波數問題，并且其數值結果要優于文[19-21]中方法結果.

4)該重心插值方法可推廣到二維和三維Helmholtz方程的求解中.現正在進行此項研究.