郭魏麗，林福榮

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

0 引言

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

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

Wiener-Hopf積分方程的數(shù)值解法可分為兩類.一類是對半無限區(qū)間進行截斷，化為求解有限區(qū)間上的積分方程.截斷方程的形式如下[3-5].

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

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

本文余下部分安排如下.第2節(jié)介紹NCC方法及其主要性質(zhì)，第3節(jié)討論用復合NCC方法離散Wiener-Hopf的截斷方程（2），第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 多項式，

是n次Chebyshev多項式Tn（s）＝cos（n arccos（s））的根.

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

則當n≥k時，有

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

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

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

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

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

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

其中右端函數(shù) g（t）和未知函數(shù) x（t）是光滑函數(shù)，h（t，s）是一個 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)鍵點是對固定的t，把（4）中的被積函數(shù)看作s的函數(shù)，在整個區(qū)間[－1，1]進行有效逼近.

定義

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

確定，即

上式中

于是

經(jīng)計算得

則由（7）和（8）可得

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

其中W＝CSLC－1，

類似地，設(shè)

其中S由（9）給出，

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

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

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

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

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

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

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

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

則（16）的離散方程為

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

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

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

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

當mn很大時，用直接法求解離散方程組的計算量很大.因此，我們考慮用迭代法求解，以提高解方程的效率.

3 離散方程組的預處理迭代解法

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

則方程組（17）化為

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

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

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

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

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

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

應用PCGNR方法求解線性方程組（20）的每步迭代的主要工作是計算和求解下面簡要說明快速算法的步驟：

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

4 數(shù)值例子

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

例4.1.考慮

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

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

例4.2.考慮

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

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

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