第二類Wiener-Hopf積分方程的復(fù)合Nystr?m-Clenshaw-Curtis方法

郭魏麗，林福榮

（汕頭大學(xué)數(shù)學(xué)系，廣東 汕頭 515063）

0 引言

本文考慮第二類Wiener-Hopf積分方程的高精度快速解法，這是一類定義在正半軸上的卷積方程：

其中 k（t－s）是核函數(shù)，f（t）為已知函數(shù)，u（t）是待求函數(shù).這類方程有廣泛的應(yīng)用背景[1]145-146，186-189，其數(shù)值方法吸引了眾多學(xué)者[2-5].

Wiener-Hopf積分方程的數(shù)值解法可分為兩類.一類是對(duì)半無限區(qū)間進(jìn)行截?cái)啵癁榍蠼庥邢迏^(qū)間上的積分方程.截?cái)喾匠痰男问饺缦耓3-5].

Gohberg等提出應(yīng)用預(yù)處理共軛梯度法求解該方程，并引進(jìn)循環(huán)預(yù)處理算子提高算法的效率[3].Lin等則提出用卷積算子作為逆預(yù)處理算子[5].Kang等提出Nystr?m-Clenshaw-Curtis（NCC）求積方法，并應(yīng)用于方程（2）的離散[4].NCC方法是求解具有半光滑核函數(shù)的積分方程的高精度方法，但離散矩陣沒有結(jié)構(gòu)，而且穩(wěn)定性也存在問題.Chen等提出了修正的NCC求積方法，使得NCC求積方法更加穩(wěn)定[6].另一類數(shù)值解法則不作截?cái)啵苯訉?duì)方程（1）進(jìn)行離散[2，7].

本文在NCC方法的基礎(chǔ)上，提出復(fù)合NCC求積公式，并對(duì)離散方程組做適當(dāng)?shù)奶幚恚沟孟禂?shù)矩陣有Toeplitz塊結(jié)構(gòu).然后引進(jìn)循環(huán)塊矩陣作為預(yù)處理矩陣提高預(yù)處理迭代方法的收斂速度.

本文余下部分安排如下.第2節(jié)介紹NCC方法及其主要性質(zhì)，第3節(jié)討論用復(fù)合NCC方法離散Wiener-Hopf的截?cái)喾匠蹋?），第4節(jié)考慮離散方程組的快速解法，第5節(jié)給出若干例子說明本文提出的算法的有效性.

1 Nystr?m-Clenshaw-Curtis求積方法

先介紹Chebyshev網(wǎng)格上的插值函數(shù)的逼近性質(zhì).設(shè)Tj（s）＝cos（j arccos（s）），j＝0，1，2，…，是一列 Chebyshev 多項(xiàng)式，

是n次Chebyshev多項(xiàng)式Tn（s）＝cos（n arccos（s））的根.

引理1.1（Jackson第五定理[8]95-96）設(shè)h∈C[－1，1]且存在正整數(shù)k使得h的k階導(dǎo)數(shù)h(k)連續(xù).用表示用Qn表示次數(shù)不超過（n－1）次的多項(xiàng)式的集合.設(shè)（Bnh）（x）是 h（x）的最佳逼近多項(xiàng)式，即

則當(dāng)n≥k時(shí)，有

Lebesgue函數(shù)和Lebesgue常數(shù)是分析插值多項(xiàng)式的誤差的有用工具.對(duì)于一個(gè)給定的插值網(wǎng)格{s1，…，sn}，Lebesgue函數(shù)定義為

是第k個(gè)Lagrange插值基函數(shù)，Lebesgue常數(shù)定義為對(duì)于Chebyshev網(wǎng)格有用（Tnh）（s）表示 h（s）在 Chebyshev網(wǎng)格上的插值多項(xiàng)式

下面引理給出（Tnh）（s）的誤差估計(jì)[9].

引理 1.2 設(shè) k 為正整數(shù)，h（s）∈Ck[－1，1].如果 n＞k，則

可見，用Chebyshev網(wǎng)格對(duì)光滑函數(shù)進(jìn)行插值非常有效.

下面介紹NCC求積方法.考慮積分方程

其中右端函數(shù) g（t）和未知函數(shù) x（t）是光滑函數(shù)，h（t，s）是一個(gè) p－半光滑核函數(shù)（p 為正整數(shù)），即

其中hi（t，s）∈Cp（[－1，1]×[－1，1]），i＝1，2.注意到

由于h1（t，s）和h2（t，s）在[－1，1]×[－1，1]中光滑，h1（t，s）x（s）和h2（t，s）x（s）也在[－1，1]×[－1，1]光滑.NCC求積方法的關(guān)鍵點(diǎn)是對(duì)固定的t，把（4）中的被積函數(shù)看作s的函數(shù)，在整個(gè)區(qū)間[－1，1]進(jìn)行有效逼近.

定義

其中系數(shù) αl，j由插值條件

確定，即

上式中

于是

經(jīng)計(jì)算得

則由（7）和（8）可得

由（5），（6），（11），（12）得

其中W＝CSLC－1，

類似地，設(shè)

其中S由（9）給出，

由（13）和（14）可得（3）的 NCC離散方程組

說明：（1）對(duì)于一般區(qū)間[a，b]上的積分方程

可以通過簡(jiǎn)單的仿射變換化為[-1，1]上的積分方程：

（2）設(shè)h1（t，s）x（s）∈Cp（[－1，1]2），h2（t，s）x（s）∈Cq（[－1，1]2），則當(dāng)NCC方法用于（15）時(shí)，數(shù)值解有如下的誤差界：其中 c 為正常數(shù)，r＝min（p，q）.

2 求解第二類Wiener-Hopf積分方程的NCC方法

考慮截?cái)郬iener-Hopf積分方程（2）的離散.必須指出，當(dāng)積分區(qū)間比較大時(shí)，直接應(yīng)用NCC方法可能導(dǎo)致核函數(shù)的函數(shù)值過大，造成很大的舍入誤差.比如像這樣的函數(shù)，當(dāng)積分區(qū)間較大時(shí)用NCC方法得不到高精度的解.復(fù)合NCC公式可以作為一種解決方法.其基本思想是對(duì)區(qū)間進(jìn)行合理的分段，然后在每個(gè)小區(qū)間上應(yīng)用NCC方法，從而避免出現(xiàn)較大的舍入誤差的現(xiàn)象.通過分析我們發(fā)現(xiàn)：如果把區(qū)間[0，T]劃分為等長(zhǎng)的子區(qū)間，應(yīng)用復(fù)合NCC公式并將求積節(jié)點(diǎn)作適當(dāng)?shù)闹嘏牛瑒t系數(shù)矩陣具有Toeplitz塊結(jié)構(gòu).這樣可以由Toeplitz塊結(jié)構(gòu)和循環(huán)塊結(jié)構(gòu)設(shè)計(jì)快速算法，顯著減少計(jì)算量.

積分方程（2）可以寫成如下積分方程組：

考慮在每個(gè)小區(qū)間[bj－1，bj]，j＝1，2，…，m，應(yīng)用NCC求積公式，每個(gè)小區(qū)間都取n個(gè)節(jié)點(diǎn).令

則（16）的離散方程為

定義如下的mn維向量f和u：

設(shè)A為m×m分塊矩陣

其中 Aij為 n×n 矩陣，其（k，l）元素為

上述方程組可以寫成如下緊湊格式：

當(dāng)mn很大時(shí)，用直接法求解離散方程組的計(jì)算量很大.因此，我們考慮用迭代法求解，以提高解方程的效率.

3 離散方程組的預(yù)處理迭代解法

由（18）定義的分塊矩陣沒有特殊的結(jié)構(gòu).經(jīng)過研究，我們發(fā)現(xiàn)對(duì)求積節(jié)點(diǎn)作適當(dāng)重排后，系數(shù)矩陣具有Toeplitz塊結(jié)構(gòu).由此我們得到求解線性方程組（17）的快速算法.文獻(xiàn)[10]研究了核函數(shù)光滑的情形.如果 Tn＝[ti，j]n×n滿足 ti，j＝ti－j，則稱 Tn是一個(gè) n 階Toeplitz 矩陣.特別地，如果 Cn＝[ci－j]n×n滿足 ci＝ci－n，i＝1，2，…，n－1，則稱 cn是一個(gè) n階循環(huán)矩陣.循環(huán)矩陣可由傅里葉矩陣對(duì)角化：

則方程組（17）化為

其中Λij為對(duì)角陣（其對(duì)角元由的第一列經(jīng)快速傅里葉變換得到），有

其中P是一個(gè)mn×mn的置換矩陣，Dk是n×n的矩陣，滿足

由于線性方程組（19）的系數(shù)矩陣可能不是對(duì)稱正定的，我們考慮對(duì)其法方程組應(yīng)用共軛梯度法（CGNR）和預(yù)處理共軛梯度法（PCGNR）.記則（19）的法方程組為

其中B*表示B的轉(zhuǎn)置.

取預(yù)處理矩陣為Q＝M*M，用PCGNR方法求解下列方程組

為完整起見，我們給出求解（20）的PCGNR方法（如果Q改為單位矩陣，則為CGNR方法）.

應(yīng)用PCGNR方法求解線性方程組（20）的每步迭代的主要工作是計(jì)算和求解下面簡(jiǎn)要說明快速算法的步驟：

因此，在求解這類積分方程時(shí)，我們可以選取不太大的n（如n＝16），這樣，每步迭代的計(jì)算量比較小.

4 數(shù)值例子

本節(jié)給出數(shù)值例子說明本文提出的方法的有效性.在CGNR和PCGNR方法中，初始解為零向量，終止條件為剩余向量的相對(duì)誤差小于10－14，即∈＝10－14.所有的計(jì)算在Dell Optiplex 9020臺(tái)式計(jì)算機(jī)上用Matlab運(yùn)行.在下面的表格中，“誤差”表示數(shù)值解的相對(duì)誤差，即其中分別表示數(shù)值解和精確解；分別表示高斯消去法，CGNR方法（矩陣-向量相乘用快速算法），CGNR方法（矩陣-向量相乘不用快速算法），PCGNR方法的計(jì)算時(shí)間，單位為秒；ICGNR，IPCGNR分別表示兩種方法的迭代次數(shù).

例4.1.考慮

取T＝80，n＝16，相關(guān)的結(jié)果如表1所示.從表1可以看出：（1）PCGNR方法的收斂比CGNR方法快得多，但在計(jì)算時(shí)間方面并沒有優(yōu)勢(shì).主要原因在于構(gòu)造預(yù)處理矩陣需要較多的時(shí)間.（2）隨著m的增大，數(shù)值解的精度提高非常快.（3）當(dāng)m較小時(shí)，快速算法的優(yōu)勢(shì)不明顯；隨著m增大，快速算法的優(yōu)勢(shì)越來越大，如m＝128時(shí)，TCGNR約為的1/9，為TGE的1/7.

表1 誤差、計(jì)算時(shí)間和迭代次數(shù)（例題4.1）

例4.2.考慮

取T＝80，n＝16，相關(guān)的結(jié)果如表2所示.我們看到與例題4.1類似的結(jié)果：（1）PCGNR方法的收斂比CGNR方法快得多，但在計(jì)算時(shí)間方面并沒有優(yōu)勢(shì).（2）隨著m的增大，數(shù)值解的精度提高非常快.（3）當(dāng)m較小時(shí)，快速算法的優(yōu)勢(shì)不明顯；隨著m增大，快速算法的優(yōu)勢(shì)越來越大.

表2 誤差、計(jì)算時(shí)間和迭代次數(shù)（例題4.2）

從以上的數(shù)值結(jié)果我們看到：復(fù)合NCC方法是一個(gè)精度非常高的方法，這與理論結(jié)果相符；利用矩陣的Toeplitz塊結(jié)構(gòu)，快速矩陣-向量乘法可以提高算法的效率；引進(jìn)預(yù)處理矩陣雖然能加快迭代方法的收斂速度，但在計(jì)算時(shí)間方面沒有優(yōu)勢(shì).如何更加高效的預(yù)處理矩陣，是一個(gè)值得進(jìn)一步研究的問題.