劉仲云，李莉

(長沙理工大學 數(shù)學與統(tǒng)計學院，湖南 長沙，410114)

文中考慮大型稀疏線性方程組

Ax=b

(1)

的迭代解法，其中x，b∈Rn，A∈Rn×n為含Toeplitz塊的全正塊三對角隨機Toeplitz矩陣，即A有如下結(jié)構：

其中:Ti(i=-1，0，1)為Toeplitz矩陣[1]。

一個非奇異矩陣稱為全正矩陣[2]，如果它所有的子式都是非負的。

一個矩陣是一個隨機矩陣，如果它是非負矩陣且每一行的行和都為1。

這類線性方程組頻繁出現(xiàn)在三次均勻B樣條曲面擬合及二階微分方程數(shù)值解等問題中。

1 精化迭代法

求解(1)式的精化迭代法的迭代格式為

x(k+1)=x(k)+d(k)

(2)

Ad(k)=r(k)，r(k)=b-Ax(k)

(3)

其中:x(0)為初始迭代向量，k=0，1，…

當d(k)為(3)的精確解時，迭代格式(2-3)退化為單步迭代精化。不難看出，該迭代格式可以提高解x(k)的精度[3]。另外，倘若用高精度方法求解(3)，則迭代格式(2-3)為迭代精化[4-5]，把該迭代格式稱為精化迭代法。

引理1 假設不考慮舍入誤差，精確計算殘差r(k)，d(k)為(3)的精確解，則x(k+1)為(1)的精確解。

證明 由(3)式知Ad(k)=r(k)，r(k)=b-Ax(k)，從而Ax(k+1)=Ax(k)+Ad(k)=b。這意味著x(k+1)是(1)的精確解，證畢。

由于A的正定性，因此，用廣義極小殘差法(GMRES)非精確求解(3)，這樣就得到精化迭代法的基本框架。

設d(k，l)為用GMRES方法求解(3)得到的第l個近似解，其中，初始迭代向量d(k，0)=0，終止檢驗公式為

(4)

其中：τk是終止判別。

倘若d(k，Lk)滿足終止檢驗公式(4)，從而

x(k+1)=x(k)+d(k，Lk)

(5)

通常，為了確保計算效率，(5)的終止檢驗公式為

‖r(k+1)‖≥θ‖r(k)‖，θ∈[0，1]。

(6)

為方便算法實現(xiàn)，將精化迭代法每次迭代的具體步驟總結(jié)為下面的流程圖(見下圖1)。

輸入初始迭代向量x(0),終止判別τk和精度θ。 Step 1.計算殘差r(k); Step 2.利用GMERES方法求解Ad(k)=r(k),尋找滿足終止檢驗公式(4)的d(k,Lk); Step 3.更新x(k+1)=x(k)+d(k,Lk); Step 4.若滿足終止檢驗公式(6),則輸出近似解x(k+1),計算結(jié)束;否則,令 x(k)=x(k+1),k=k+1,轉(zhuǎn)Step 1。

圖1精化迭代法每次迭代的具體步驟
Fig.1Specific steps of each iteration of the refined iterative method

顯然，對于不同的τk值，d(k，Lk)不同，精化迭代法的收斂速度也不同。當τk等于0時，可以把算法1得到的解認為是公式(1)的精確解。因此，當選取的τk值較小時，迭代序列有較好的收斂速度。

定理2 令x*為線性方程組(1)的精確解。對第k步迭代，用GMRES方法求解(3)的初始迭代向量均為d(k，0)=0，相應的迭代序列為{d(k，l)}Lkl=0且d(k，Lk)滿足終止檢驗公式(4)，則迭代序列{x(k)}k=0滿足

‖x(k+1)-x*‖≤τkκ(A)‖x(k)-x*‖，

其中:κ(A)=‖A‖·‖A-1‖。

特別地，當

ηmax=τmaxκ(A)<1，

迭代序列{x(k)}k=1收斂到x*，其中

證明 由(4)式知

x(k+1)-x*=x(k)-x*+d(k，Lk)=
A-1[Ad(k，Lk)-rk]，

從而

‖x(k+1)-x*‖≤‖A-1‖‖Ad(k，Lk)-rk‖≤τk‖A-1‖‖rk‖≤
τk‖A-1‖·‖A‖·‖x(k)-x*‖=τkκ(A)‖x(k)-x*‖≤ηmax‖x(k)-x*‖。

當ηmax<1時，迭代序列{x(k)}k=1收斂，證畢。

2 數(shù)值實驗

下面給出精化迭代法應用于曲面擬合的兩個例子。采用三次均勻B樣條曲面插值相應的數(shù)據(jù)點，利用精化迭代法求相應的控制點，從而得到插值給定點集的曲面。第k步迭代得到的三次均勻B樣條曲面的擬合誤差計算式為

右端向量b等于數(shù)據(jù)點{Pij}(i=0，1，…，m；j=0，1，…，n)，系數(shù)矩陣均為A=B1?B2，符號?表示Kronecker積，其中B1∈Rm×m，B2∈Rn×n為配置矩陣，形如

例1 利用三次均勻B樣條曲面逼近peak函數(shù)

中的點

例2 利用三次均勻B樣條曲面逼近

中的點

利用三次均勻B樣條曲面迭代逼近例1、2中的點。所有的數(shù)值實驗是在Matlab7.4.0.287(R2014a)環(huán)境中實現(xiàn)。在數(shù)值實驗中，取τk=10-2，ε=10-8。為了比較，還測試了文獻[6]中加權漸進迭代逼近法(WPIA)。表1、2給出了系數(shù)矩陣的階數(shù)不同時，三次均勻B樣條曲面的精化迭代法和WPIA逼近例1、2中的點所需迭代次數(shù)和計算時間。

表1 精化迭代法與加權漸進迭代逼近法求解例1

Table 1 Comparisons between the refined iterative method and WPIA for Case 1

矩陣階數(shù)/n加權漸進迭代逼近迭代次數(shù)/次Time/s精化迭代法迭代次數(shù)/次Time/s相對效率/%900591.24E-01323.93E-0268.31 600616.85E-01292.12E-0169.02 500621.71E+00294.94E-0171.13 600623.67E+00227.40E-0179.84 900636.63E+00201.29E+0080.5

表2 精化迭代法與加權漸進迭代逼近法求解例2

Table 2 Comparisons between the refined iterative method and WPIA for Case 2

矩陣階數(shù)/n加權漸進迭代逼近迭代次數(shù)/次Time/s精化迭代法迭代次數(shù)/次Time/s相對效率/%900511.32E-01313.92E-0270.21 600546.69E-01211.55E-0176.92 500561.55E+00213.28E-0178.93 600563.25E+00206.67E-0179.54 900565.81E+00201.37E+0076.5

在表1、2中，n為系數(shù)矩陣A的階數(shù)；Time為單獨計算100次的平均計算時間，相對效率計算公式為

其中：tw為加權漸進迭代逼近法的計算時間；tr為精化迭代法的計算時間。

由表1、2可看出，在迭代次數(shù)和計算時間上，精化迭代法都要優(yōu)于加權漸進迭代逼近法，且矩陣規(guī)模越大，優(yōu)勢就愈明顯。