楊苗苗,葛永斌
(寧夏大學(xué)數(shù)學(xué)統(tǒng)計學(xué)院,寧夏銀川750021)
Helmholtz方程是一個橢圓型偏微分方程,代表了波動方程與時間無關(guān)的解.這個方程模擬了各種各樣的物理現(xiàn)象.其中包括:振動分析、水波傳播、聲波散射、電磁散射、以及雷達散射等問題.隨著Helmholtz方程被廣泛地應(yīng)用在工程和科學(xué)領(lǐng)域的研究中,該方程本身的復(fù)雜特性給數(shù)值計算帶來了巨大困難,特別是對于變波數(shù)和大波數(shù)問題的有效數(shù)值計算,雖然其數(shù)值解法有很多,但都需要進一步提高計算精度和計算效率,以適應(yīng)實際復(fù)雜大規(guī)模問題的求解.所以,對Helmholtz方程的高效精確數(shù)值求解方法的研究具有重要的理論價值和實際意義.
無網(wǎng)格法[1]是近年來興起的一種新型、高效的數(shù)值計算方法.它基于點的近似思想來構(gòu)造近似函數(shù),適用于分析各類具有大梯度、奇異性等特殊性質(zhì)的應(yīng)用問題.目前發(fā)展的無網(wǎng)格方法包括無單元Galerkin法[2]和無網(wǎng)格配點法[3]等.對于無網(wǎng)格配點法,重心插值配點法就是其中的一種.重心插值配點法包括重心Lagrange插值法和重心有理插值法,是指未知函數(shù)的近似函數(shù)在插值節(jié)點處的重心型插值,是一種依賴于偏微分方程的強形式的配點法.該數(shù)值方法具有易于理解,思想結(jié)構(gòu)簡單,計算精度高,便于編程的優(yōu)點.近年來,劉婷和馬文濤[4]采用重心Lagrange插值配點法求解了二維電報方程,得到了一種高效的數(shù)值方法.王彩珍等[5]利用重心插值配點法,構(gòu)造了二維非線性橢圓型方程所對應(yīng)的離散方法,得到了一種計算效率更高的數(shù)值解法.
對于Helmholtz方程的求解,近年來研究者們提出了很多無網(wǎng)格方法,如李鵬等[6]直接利用最小二乘法建立系統(tǒng)的變分公式,推導(dǎo)了大波數(shù)Helmholtz問題的加權(quán)無網(wǎng)格公式.何锃等[7]通過使用加權(quán)正交基函數(shù)構(gòu)造了改進的移動最小二乘法,并結(jié)合拉格朗日乘子法離散Dirichlet邊界條件,推導(dǎo)了求解該方程的無網(wǎng)格加權(quán)最小二乘法.陳黃鑫和邱偉峰[8]提出了一種求解大波數(shù)Helmholtz方程的一階系統(tǒng)最小二乘法,并利用一個非平凡分解來解決對偶的一階系統(tǒng)最小二乘法問題.楊子樂等[9]利用無網(wǎng)格節(jié)點MIP法的d適應(yīng)性和移動最小二乘法來構(gòu)造近似函數(shù),給出了求解Helmholtz方程的兩種計算方法.李瑩等[10]采用無網(wǎng)格局部徑向基點插值法,將徑向基函數(shù)耦合多項式基函數(shù)作為近似函數(shù),克服了無網(wǎng)格法中場函數(shù)在計算中的冗余性.Assari等[11]將局部支撐薄板樣條的離散配點法近似化為一類具有自由形狀參數(shù)的徑向基函數(shù)法,用一種特殊的非等距Gauss Legendre求積法來計算出現(xiàn)的對數(shù)型奇異積分,發(fā)展了一種求解具有對數(shù)積分方程的無網(wǎng)格法.李美香等[12]將三角函數(shù)作為基函數(shù),并在局部支持域內(nèi)構(gòu)造了一種形函數(shù),研究了帶有邊界層問題和波傳播問題的Helmholtz方程.張偉和丁睿[13]通過引入全局徑向基函數(shù)和正定緊支徑向基函數(shù)構(gòu)造了求解Helmholtz方程的Galerkin型的無網(wǎng)格方法,使得求解方法有效且具有高精度.屈文鎮(zhèn)等[14]發(fā)展了一種新的局部基本解方法,與傳統(tǒng)的基本解法相比,該方法根據(jù)節(jié)點將計算域劃分為若干重疊子域后得到了稀疏的帶狀矩陣,可以求解大波數(shù)Helmholtz問題.陳林沖和李小林[15]基于無網(wǎng)格和無邊界分析,利用Burton-Miller公式提出了一種無邊界方法,該無網(wǎng)格方法對所有的波數(shù)都能求出唯一解,并且只涉及到邊界節(jié)點,還可以處理特大波數(shù)下的Helmholtz問題.
本文針對如下一維Helmholtz方程:
其中[a,b]是求解域,u(x)是關(guān)于x的未知函數(shù).k(x)表示波數(shù),可以是常數(shù),也可以是關(guān)于x的函數(shù).f(x)是源項.定義Dirichlet邊界條件為
借助Chebyshev插值節(jié)點和三種測試節(jié)點,運用重心Lagrange插值基函數(shù)和重心有理插值基函數(shù)推導(dǎo)了求解該類方程的兩種無網(wǎng)格配點法.該重心插值方法理論簡單,需要的插值節(jié)點少,最后得到的離散矩陣處理方便,不僅可以計算大波數(shù)問題,還可以計算變波數(shù)問題.具有推導(dǎo)過程簡單,計算結(jié)果精度高和穩(wěn)定性好等的優(yōu)點.主要內(nèi)容如下:第二節(jié)推導(dǎo)一維重心Lagrange插值和重心有理插值的計算公式和一些簡化矩陣的理論知識.第三節(jié)推導(dǎo)計算一維Helmholtz方程的兩種插值配點法及其離散矩陣.第四節(jié)給出一些測試問題的計算結(jié)果,并與一些文[19-21]中結(jié)果進行對比,最后給出本文的結(jié)論.
在區(qū)間[a,b]上插入n+1個節(jié)點xi,則有a=x0<x1<x2<···<xn=b,設(shè)未知函數(shù)u(x)在節(jié)點處的函數(shù)值ui=u(xi),i=0,1,···,n,則Lagrange多項式插值公式可表示為:
其中Li(x)為Lagrange插值基函數(shù),具體形式為:
令
定義重心Lagrange插值函數(shù)的權(quán)重為:
此時,插值基函數(shù)(2.2)式可變?yōu)?
將(2.5)代入(2.1)式,得:
則上式被稱作改進的Lagrange插值公式.相比Lagrange插值公式,上式在插值節(jié)點的數(shù)目增加時,計算量減少很多,由O(n2)變?yōu)镺(n).
利用(2.6)式插值常數(shù)1,即給定一組節(jié)點為(xi,1),i=0,1,···,n,得:.
化簡后可得:
將(2.7)代入(2.6)式中可得:
則上式被稱為重心Lagrange插值公式[16],相比式(2.6),式(2.8)在保持相同計算量O(n)的同時,還能有效地克服Runge現(xiàn)象.
重心有理插值和重心Lagrange插值相類似,重心有理插值即是將Lagrange函數(shù)換為有理函數(shù)進行插值.下面我們構(gòu)造重心有理插值基函數(shù).
首先,選擇一個正整數(shù)d,且0≤d≤n.對于節(jié)點xi,i=0,1,···,n-d,選擇d+1個節(jié)點(xi,ui),(xi+1,ui+1),···,(xi+d,ui+d)后,我們利用重心Lagrange插值多項式的形式,構(gòu)造
然后,根據(jù)pi(x)構(gòu)造一個權(quán)函數(shù)
最后,利用pi(x)和λi(x)構(gòu)造出重心有理插值函數(shù)
利用常數(shù)1的Lagrange插值公式,有
由此可得
將(2.11)代入(2.9)中可得重心有理插值公式:
則上式被稱為重心有理插值公式[17].需要說明的一點是,文[17]中已經(jīng)證明式(2.13)構(gòu)造的有理插值的誤差階為O(hd+1).其中h是等距節(jié)點網(wǎng)格步長.因此只要選擇合適的參數(shù)d,式(2.13)就可達到更高的插值精度,從而得到更小的誤差.
從重心Lagrange插值權(quán)函數(shù)的表達式(2.4)和重心有理插值權(quán)函數(shù)的表達式(2.10)能夠看出,插值權(quán)的選取跟節(jié)點的分布有關(guān).因此我們考慮了三種不同的節(jié)點來進行函數(shù)插值和數(shù)據(jù)測試,它們分別是隨機節(jié)點、等距節(jié)點和Chebyshev節(jié)點.由于在進行具體的程序計算時,隨機節(jié)點和等距節(jié)點可直接調(diào)用命令,因此,這部分只給出Chebyshev節(jié)點的計算公式:
第一類Chebyshev節(jié)點:(Chebyshev I)
第二類Chebyshev節(jié)點:(Chebyshev II)
需要說明的是,對于兩類Chebyshev節(jié)點的公式,其定義區(qū)間為[-1,1],當(dāng)針對一般的區(qū)間[a,b]時,利用區(qū)間的坐標(biāo)變換公式進行轉(zhuǎn)化即可.
微分矩陣最初是在研究Chebyshev擬譜方法[18]時提出的,重心插值微分矩陣能夠直接獲得Helmholtz方程中未知函數(shù)在計算節(jié)點處的導(dǎo)函數(shù),因此,微分矩陣是重心插值配點法求解里非常重要的一部分.
由式(2.8),函數(shù)的重心Lagrange插值公式可表示為:
通過對比我們可以看出,兩種重心插值公式的根本區(qū)別是權(quán)函數(shù)的不同.不難發(fā)現(xiàn),重心有理插值除了必要的計算外,還需選取參數(shù)d.除此之外,兩種重心插值法在形式上是一致的.因此,利用(2.16)和(2.17)式,函數(shù)在節(jié)點xj(j=1,2,···,n)處的m階導(dǎo)數(shù)可以分別表示為:
寫成矩陣的形式為
方法I重心Lagrange插值方法
將重心Lagrange插值公式(2.16)代入到式(1.1)中,可得:
設(shè)xj,j=1,2,···,n,為給定的插值節(jié)點,令上式方程在所有xj上都成立,則有
注意到Li(xj)=δij,利用重心插值微分矩陣,則方程(3.2)可寫為矩陣形式:
進一步可寫為:
將式(2.16)代入到式(1.2)中,可得:
其中,分別表示n階單位矩陣的第一行和第n行,g1,g2分別為具體的值.聯(lián)立求解方程(3.4)和(3.5),即可求得在重心Lagrange插值公式下插值節(jié)點上的函數(shù)值ui,回代入式(2.16)中并利用測試節(jié)點進行計算可得本文方法I的數(shù)值解u(x).
方法II重心有理插值方法
同理,將重心有理插值公式(2.17)代入到式(1.1)中,可得:
在本節(jié)中,我們將兩種重心插值配點法應(yīng)用于幾個測試問題中以證明本文方法的有效性和準(zhǔn)確性.我們給出了包括大波數(shù)和變波數(shù)Helmholtz方程的數(shù)值算例,采用本文方法進行計算并與文[19-21]中所得的數(shù)值結(jié)果進行比較.其中L∞范數(shù)誤差的定義如下給出:
上式中ui和uei分別代表數(shù)值解和精確解.Nt是檢測點個數(shù).
問題1[19-20]:
考慮如下非齊次大波數(shù)Helmholtz方程
其中邊界條件為:
該問題的精確解為:
首先,由于方法II在進行計算中需要優(yōu)先選擇參數(shù)d,并且0≤d≤N.其中N是插值點個數(shù).因此表1給出當(dāng)Nt=30,k=10時不同d和N下問題1中方法II的L∞范數(shù)誤差.其中Nt是檢測點個數(shù),并且插值節(jié)點的類型為第二類切比雪夫點(Chebyshev II),測試節(jié)點也為第二類切比雪夫點.由數(shù)值結(jié)果我們可以看出,隨著插值點數(shù)N的增大,當(dāng)d越接近于N,方法II的計算結(jié)果精度越高,誤差越小.因此在本文剩余算例關(guān)于方法II的計算中,我們統(tǒng)一選取d=N-1.這樣既保證了方法II具有更高的精度,又得到了更小的計算誤差.
表1 當(dāng)Nt=30,k=10時不同d和N下問題1中本文方法II的L∞范數(shù)誤差
其次,表2給出了當(dāng)Nt=N,k=10時四種插值節(jié)點下問題1的L∞范數(shù)誤差.由數(shù)值結(jié)果可以得出對于四種插值節(jié)點:隨機節(jié)點、等距節(jié)點、Chebyshev I和Chebyshev II節(jié)點來說,隨機插值節(jié)點的數(shù)值結(jié)果很差;等距節(jié)點的結(jié)果一般;Chebyshev I節(jié)點的結(jié)果遠不如Chebyshev II節(jié)點,并且Chebyshev II節(jié)點的精度最高、穩(wěn)定性最好、計算誤差最小.因此對于插值節(jié)點的選擇,我們對于本文剩余算例中的表格和圖都默認選擇Chebyshev II節(jié)點,簡記為Chebyshev節(jié)點.此外對于兩種插值方法:重心Lagrange插值法和重心有理插值法來說,所得的計算結(jié)果相差不大,即精度一樣,穩(wěn)定性一樣,誤差量級一樣并且L∞范數(shù)誤差也近似相等.
表2 當(dāng)Nt=N,k=10時四種插值節(jié)點下問題1的范數(shù)誤差
然后,表3給出了當(dāng)Nt=50,k=10時三種測試節(jié)點下的L∞范數(shù)誤差.由數(shù)值結(jié)果我們可以得出對于三種測試節(jié)點:隨機節(jié)點、等距節(jié)點和Chebyshev節(jié)點來說,所得的誤差量級一樣,并且誤差也相差不大.此外對于兩種插值方法來說,所得的計算結(jié)果也相差不大.
表3 當(dāng)Nt=50,k=10時三種測試節(jié)點下問題1的L∞范數(shù)誤差
接著,表4給出了當(dāng)N=24,k=10時三種測試節(jié)點下的L∞范數(shù)誤差.由數(shù)值結(jié)果我們可以得出兩種重心插值方法隨著檢測點數(shù)Nt的變化,所得的誤差量級不變,并且L∞范數(shù)誤差也相差不大.也就是說檢測點的個數(shù)對于本文所提方法的誤差影響不大,因此在剩余算例的計算中,我們可以取任意合理的檢測點個數(shù).
表4 當(dāng)N=24,k=10時三種測試節(jié)點下問題1的L∞范數(shù)誤差
最后,表5給出了當(dāng)Nt=200,k=1時本文兩種無網(wǎng)格配點法與文[19]和文[20]中方法計算L∞范數(shù)誤差的比較.結(jié)果表明當(dāng)N=20時,本文所提兩種方法的誤差量級均為O(10-15),而文[19-20]中的誤差量級僅僅只達到O(10-5).而當(dāng)N=160時本文所提方法的誤差量級均為O(10-13),而文[19-20]中的誤差量級分別為O(10-8)和O(10-12).另外,需要說明是,盡管本文所提方法的計算精度會隨著插值點數(shù)N的增加而減小,但本文方法的誤差量級仍然要高于文[19-20]中所提方法.此外,與有限差分法相比,無網(wǎng)格重心插值配點法并不意味著插值節(jié)點數(shù)越多,所得的計算結(jié)果會越好,有時候我們僅僅可能只用很少的插值節(jié)點便可得到更好的數(shù)值結(jié)果.表6給出了當(dāng)Nt=200,k=1000時本文兩種無網(wǎng)格配點法與文[20]中方法L∞范數(shù)誤差的比較.結(jié)果表明本文所提兩種無網(wǎng)格方法的計算結(jié)果要遠遠優(yōu)于文[20]中所提方法.因此本文方法更適合求解大波數(shù)的Helmholtz問題.
表5 當(dāng)Nt=200,k=1時三種測試節(jié)點下問題1的L∞范數(shù)誤差
表6 當(dāng)Nt=200,k=1000時不同測試節(jié)點下問題1的L∞范數(shù)誤差
問題1在Nt=24,Nt=N,k=10時隨機點(a),等距點(b),Chebyshev I點(c)及Chebyshev II點(d)的四種插值節(jié)點圖可參見圖1.由圖可得:隨機點不規(guī)則的分布在區(qū)間內(nèi),等距點均勻分布在區(qū)間內(nèi),兩類Chebyshev點在端點處分布點數(shù)較多,中間分布點數(shù)相對較少.隨后,在Nt=200,k=1000,N=160時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖參見圖2.由圖可得三種測試點下的精確解和數(shù)值解吻合的很好.
圖1 問題1在N=24,Nt=N,k=10時隨機點(a),等距點(b),Chebyshev I點(c)及Chebyshev II點(d)的插值節(jié)點圖
圖2 問題1在Nt=200,k=1000,N=160時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖
問題2[20]:
考慮如下非齊次變波數(shù)Helmholtz方程
其中邊界條件為:
該問題的精確解為:
表7給出了當(dāng)Nt=100,k(x)=時本文兩種無網(wǎng)格配點法與文[20]中方法L∞范數(shù)誤差的比較.結(jié)果表明當(dāng)N=8時本文所提兩種方法的誤差量級均為O(10-6),而文[20]中的誤差量級為O(10-4).而當(dāng)N=64時本文所提方法的誤差量級為O(10-14)或O(10-15),而文[20]中的誤差量級為O(10-10).表8給出了當(dāng)Nt=200,k(x)=時本文兩種無網(wǎng)格配點法與文[20]中方法L∞范數(shù)誤差的比較.結(jié)果表明本文所提兩種無網(wǎng)格方法的計算結(jié)果要優(yōu)于文[20]中所提方法.因此本文方法更適合求解變波數(shù)的Helmholtz問題.
表7 當(dāng)Nt=100,k(x)=時不同測試節(jié)點下問題2的L∞范數(shù)誤差
表7 當(dāng)Nt=100,k(x)=時不同測試節(jié)點下問題2的L∞范數(shù)誤差
N文[20]隨機節(jié)點等距節(jié)點Chebyshev節(jié)點方法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 當(dāng)Nt=200,k(x)=時不同測試節(jié)點下問題2的L∞范數(shù)誤差
表8 當(dāng)Nt=200,k(x)=時不同測試節(jié)點下問題2的L∞范數(shù)誤差
N文[20]隨機節(jié)點等距節(jié)點Chebyshev節(jié)點方法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)下精確解與數(shù)值解圖參見圖3.由圖可得三種測試點下的精確解和數(shù)值解吻合得非常好.
圖3 問題2在Nt=200,k(x)==64時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖
問題3[21]:
考慮如下齊次Helmholtz方程
其中邊界條件為:
該問題的精確解為:
表9給出了當(dāng)Nt=100,b=0.64時不同測試節(jié)點及不同波數(shù)k=0.5N下本文兩種無網(wǎng)格配點法的L∞范數(shù)誤差.需要說明的是,文[21]中只給出了一些數(shù)值解和精確解的吻合圖而沒有給出相應(yīng)的數(shù)據(jù),因此表9只給了本文所提兩種無網(wǎng)格配點法的計算結(jié)果而沒有和文[21]的對比數(shù)據(jù).計算結(jié)果表明,即使k隨著N的變化而變化,但是本文的兩種無網(wǎng)格配點法仍然精度很高,并且誤差值很小,也就是說本文重心無網(wǎng)格方法更逼近于精確解.
表9 當(dāng)Nt=100,b=0.64時不同測試節(jié)點及不同波數(shù)k=0.5N下問題3的L∞范數(shù)誤差
問題3在N=40,k=0.5N,Nt=100,b=0.64時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖參見圖4.由圖可得三種測試點下的精確解和數(shù)值解吻合的非常好.
圖4 問題3在N=40,k=0.5N,Nt=100,b=0.64時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖
問題4:
考慮如下齊次Helmholtz方程
其中邊界條件為:
該問題的精確解為:
表10給出了當(dāng)Nt=100,k=24時不同測試節(jié)點下本文兩種無網(wǎng)格配點法的L∞范數(shù)誤差.由數(shù)值結(jié)果可以看出,該區(qū)間上本文所提兩種無網(wǎng)格方法的誤差量級最高為O(10-6),盡管誤差量級與前面幾個算例比起來不高,但是穩(wěn)定性仍然很好.
表10 當(dāng)Nt=100,k=24時不同測試節(jié)點下問題4的L∞范數(shù)誤差
問題4在Nt=100,k=24,N=56時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖參見圖5.由圖可得三種測試點下的精確解和數(shù)值解吻合的非常好.
圖5 問題4在Nt=100,k=24,N=56時方法I與方法II在隨機點(a),等距點(b)和Chebyshev測試點(c)下精確解與數(shù)值解圖
可以得出以下結(jié)論:
1)兩種重心插值配點法的計算精度相同,誤差也幾乎相同,因此對于一維Helmholtz方程來說,兩種方法的結(jié)果相差不大.
2)在插值節(jié)點的選取上,數(shù)值結(jié)果表明:隨機節(jié)點很差,等距節(jié)點一般,兩種Chebyshev節(jié)點結(jié)果很好,并且Chebyshev II要優(yōu)于Chebyshev I.在測試節(jié)點的選取上,三種測試節(jié)點的數(shù)值結(jié)果幾乎沒有差異.此外,隨著測試節(jié)點數(shù)量的變化,本文方法的誤差幾乎沒有變化.也就是說,測試節(jié)點的數(shù)量對我們的數(shù)值結(jié)果影響不大.
3)與文[19-21]中所提方法相比,本文所提兩種無網(wǎng)格配點法使用較少的節(jié)點即可獲得高精度的數(shù)值結(jié)果,并保持良好的穩(wěn)定性.并且數(shù)值實驗研究結(jié)果進一步表明,本文方法具有所需配點數(shù)少、精度高、誤差小和效率高等優(yōu)點.此外,本文方法不僅可以計算大波數(shù)問題,還可以計算變波數(shù)問題,并且其數(shù)值結(jié)果要優(yōu)于文[19-21]中方法結(jié)果.
4)該重心插值方法可推廣到二維和三維Helmholtz方程的求解中.現(xiàn)正在進行此項研究.