李金，張曉蕾，桑瑜，張宇鑫

(華北理工大學(xué) 理學(xué)院，河北 唐山 063210)

引言

柯西主值積分是邊界元法中的重要研究內(nèi)容，被廣泛地應(yīng)用于數(shù)學(xué)物理、電磁力學(xué)、流體力學(xué)、斷裂力學(xué)、化學(xué)等領(lǐng)域[1,2]?？紤]如下形式的柯西主值積分：

(1)

奇異積分的近似計算作為數(shù)值計算領(lǐng)域的研究熱點。目前，已經(jīng)有矩形求積公式[3,4]、梯形求積公式[5,6]、辛普森求積公式[7]、高斯求積公式[8]等方法來近似計算奇異積分。LI用矩形求積公式近似計算了二維奇異積分[9]，得到收斂結(jié)果為O(h2)。文獻(xiàn)[10]中討論了計算柯西主值積分的復(fù)合矩形求積公式，逼近密度函數(shù)得到誤差展開式，設(shè)計外推算法提高收斂階獲得后驗誤差估計。隨后有學(xué)者研究了多維柯西主值積分的復(fù)合矩形求積公式[11]，并給出外推算法。矩形求積公式還可以用來近似計算Hadamard有限部分積分[12]。LI等學(xué)者利用復(fù)合梯形公式近似計算一維、二維及圓上的柯西主值積分[13-15]，得到相應(yīng)的超收斂誤差估計。Monegato給出了近似計算二維柯西主值積分的一類插值求積公式[16]。Kim構(gòu)造了一個基于三角變換的求積公式[17]來近似計算柯西主值積分。

外推法作為一種提高計算精度的方法，近年來被廣泛應(yīng)用于數(shù)值計算[18]。在近似計算奇異積分時，根據(jù)得到的誤差展開式設(shè)計外推算法，從而提高誤差收斂階。在等距節(jié)點的情況下，通常采用把區(qū)間逐次對分的方法計算積分值，這樣，前一次劃分得到的函數(shù)值在區(qū)間重新劃分后仍可被利用，且易于編程，這是外推法實現(xiàn)的前提?；谝陨纤悸?，Navot[19]構(gòu)造了基于Fourier級數(shù)展開的弱奇異積分的Euler-Maclaurin展開式。Lyness研究了二維柯西主值積分的Euler-Maclaurin展開方法[20]。余德浩等學(xué)者根據(jù)梯形求積公式近似計算二階超奇異積分得到誤差漸近展開式，構(gòu)造廣義的外推算法來提高計算精度[21]。

該研究對復(fù)合矩形求積公式近似計算二維柯西主值積分進(jìn)行分析。對積分區(qū)域進(jìn)行均勻剖分，逼近密度函數(shù)與核函數(shù)，通過矩形公式近似計算得到誤差展開式。構(gòu)造序列來逼近奇異點，并對網(wǎng)格進(jìn)行均勻加密，提出了一種外推算法。通過該方法獲得更高的收斂階和后驗誤差估計。該項研究的第一部分給出了矩形求積公式的基本定義及其近似計算二維柯西主值積分獲得誤差展開式的相關(guān)定理；第二部分證明主要定理，并根據(jù)誤差展開式得到收斂階；第三部分在計算二維柯西主值積分誤差展開式的基礎(chǔ)上，設(shè)計外推算法提高收斂階，并給出后驗誤差估計；第四部分用數(shù)值算例來驗證理論分析的正確性。

1主要定理

令a=x0

定理1[22]f(x,y)∈C∞[a,b]×[c,d]，θ1,θ2∈[0,1]。定義

(2)

和

(3)

然后有

(4)

其中ck是與h無關(guān)的常數(shù)。

接下來定義fc(x,y)為f(x,y)的常數(shù)插值

fc(x,y)=f(xi-1,yj-1),(x,y)∈[xi-1,xi]×[yj-1,yj],i,j=1,2,…,n

(5)

定義線性變換

(6)

(7)

(8)

其中ωij定義為科特斯系數(shù)，

(9)

定義φ0(x,y)，φk1(x,y)，φk2(x,y)，

(10)

(11)

(12)

當(dāng)k=1時，

(13)

(14)

現(xiàn)在給出如下定理。

定理2設(shè)f(x,y)∈C∞(Ω)=C∞[a,b]×[c,d]，矩形求積公式In(f,t,s)定義為式(8)，存在一個與h無關(guān)的常數(shù)ck,使得

(15)

其中(t,s)∈(xg-1,xg)×(yl-1,yl)，

(16)

(17)

(18)

(19)

(20)

當(dāng)k=1時，

(21)

(22)

當(dāng)式(18)為0時出現(xiàn)超收斂現(xiàn)象，收斂階為O(h)。

2定理的證明

(23)

(24)

(25)

證明:根據(jù)二維柯西主值積分的定義和線性變換式(6)和式(7)，有

(26)

對于式(24)的第一部分

(27)

相似的可以得到式(25)。i≠g，j≠l的情況，相應(yīng)的黎曼積分可以被同樣證明。當(dāng)k=1時有，

(28)

(29)

以上可以被類似證明。

引理2在定理2的假設(shè)下，有

(30)

證明:通過將fc(x,y)，f(x,y)在奇點(t,s)處泰勒展開，得到

(31)

(32)

結(jié)合式(31)和式(32)，得到引理2結(jié)論。定理2的證明如下

證明:根據(jù)引理2，有

(33)

對于i=g，j=l，有

(34)

結(jié)合式(33)和式(34)可得

(35)

這里d0(φ0)，dk1(φk1)，dk2(φk2)定義為式(18)～式(20)。通過線性變換[xi-1,xi]×[yj-1,yj]映射到[-1,1]×[-1,1]，當(dāng)τ=ξ=0時，d0(φ0)=0。對于最后一部分沒有奇異性，

(36)

(37)

證明完畢。

實際上獲得了矩形公式計算二維柯西主值積分的誤差展開式，因此當(dāng)d0(φ0)=0時，得到超收斂點。接下來將以上結(jié)論推廣到多維的情況。

(38)

這里(t1,t2,...,tq)是奇點，令f(t1,...,tq)∈C∞[ai,bi]q，ai=xi0

fc(x1,x2,...,xq)=f(x1,i1-1,x2,i2-1,...,xq,iq-1)

(39)

然后得到：

(40)

其中科特斯系數(shù)為ωi1,i2,...,iq=hq/[(x1,i1-1-t1)...(xq,iq-1-tq)]。

3外推法

根據(jù)定理2得到矩形公式計算二維柯西主值積分的誤差展開式

(41)

并給出如下外推算法，存在正整數(shù)n10,n20使得m10:=n10(t-a)/(b-a)，m20:=n20(s-c)/(d-c)為常數(shù)。為了簡化計算過程，設(shè)n10=n20，將[a,b]×[c,d]劃分為等子區(qū)域，得到大小為h1=(b-a)/n10=(d-c)/n20的初始網(wǎng)格Π1，細(xì)化網(wǎng)格Π1得到Π2，其大小為h2=h1/2。通過不斷細(xì)化網(wǎng)格得到{Πj}j=1,2,...,這里Πj由Πj-1得到，網(wǎng)格大小為hj，得到如下所示的外推表。

表1 外推表Ti(j)

對于給定的τ,ξ,定義

(42)

(43)

T(hi)=T(hi,hj)=I2i-1n10,2j-1n20(f,ti,sj)

(44)

下面介紹外推算法步驟

第一步:

第二步:

定理3在定理2的誤差展開式基礎(chǔ)上，對于τ,ξ=0，根據(jù)式(42)和(43)有

(45)

后驗誤差估計為

(46)

證明:對于給定的τ,ξ,通過式(41)的誤差展開式有

(47)

根據(jù)柯西主值積分的定義和式(41)，將I(f,ti,sj)在奇點(t,s)處進(jìn)行泰勒展開，有

(48)

f(ti,sj)在奇點(t,s)處進(jìn)行泰勒展開得到

(49)

結(jié)合式(47)～式(49)得到

(50)

這里

(51)

對于給定的τ,ξ和bk(t,s,τ,ξ)均為常數(shù)，由式(50)有

(52)

根據(jù)式(50)和式(52)以及hi=2hi+1,可得

(53)

即

(54)

這里

(55)

繼續(xù)外推過程可以得到收斂階O(h3)，相似的得到收斂階O(h4)。通過這種方法，可以得到更高的收斂階。

4數(shù)值算例

這一部分給出一些數(shù)值例子驗證定理的結(jié)論。

例1：

由表2可得，當(dāng)奇點局部坐標(biāo)為零時收斂階為O(h)，當(dāng)局部坐標(biāo)不為零時則不收斂。因此當(dāng)τ=ξ=0時，出現(xiàn)超收斂現(xiàn)象并驗證了理論分析的正確性。

表2 矩形公式誤差估計

例2：

表3 矩形公式外推誤差估計

表4 后驗誤差估計

例3：

表5 矩形公式外推誤差估計

表6 后驗誤差估計

例4：

表7和表8驗證了當(dāng)τ=ξ=η=0時，得到矩形公式近似計算三維柯西主值積分并進(jìn)行外推的誤差，收斂階可以達(dá)到O(h),O(h2),O(h3),O(h4)，后驗估計的收斂階與表7的誤差收斂階相同，符合理論分析。

表7 矩形公式外推誤差估計

表8 后驗誤差估計

5結(jié)論

通過矩形求積公式近似計算二維柯西主值積分得到誤差展開式。通過構(gòu)造序列來逼近奇點，設(shè)計外推算法獲得更高的計算精度，并得到了后驗誤差估計。這種方法也可以推廣到多維柯西主值積分。