帶BTTB矩陣線性互補問題的塊預(yù)處理模系矩陣分裂迭代方法

吳敏華， 李郴良

(桂林電子科技大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院，廣西 桂林 541004)

彈塑性接觸模型根據(jù)力學(xué)現(xiàn)象可以轉(zhuǎn)化為線性互補問題，接觸面經(jīng)離散后的系數(shù)矩陣往往非常大且具有Toeplitz結(jié)構(gòu)。Zhao等[1]利用多重網(wǎng)格求解，將多重網(wǎng)格方法嵌套到有效集，提出了完全多重網(wǎng)格方法。Vollebregt[2]基于共軛梯度法和快速傅里葉變換，提出了BCCG+FAI方法，同時給出了一個新的預(yù)處理算子。Bai[3]將線性互補問題轉(zhuǎn)化為不動點問題，將矩陣適當分裂，提出了更具一般性的模系矩陣分裂迭代方法。Bai等[4]還利用并行技術(shù)，提出了模系矩陣多分裂方法。1986年Strang[5]引入循環(huán)矩陣作為預(yù)優(yōu)共軛梯度法的預(yù)處理矩陣來求解Toeplitz系統(tǒng)Tx=b，從此，基于預(yù)優(yōu)共軛梯度法的Toeplitz迭代解法在過去的幾十年里有了飛速發(fā)展[6-8]。

用模系矩陣迭代方法求解大型線性互補問題十分高效。為此，將預(yù)優(yōu)共軛梯度法和模系矩陣迭代方法結(jié)合，利用系數(shù)矩陣是正定的二級對稱BTTB矩陣的優(yōu)點，構(gòu)造BCCB塊預(yù)處理算子[8]，用共軛梯度法求解線性互補問題。

1 預(yù)備知識

對于線性互補問題，求一對可行的互補解w,z∈Rn，使

w=Tz+q≥0,z≥0,zTw=0，

(1)

其中T∈Rn×n為對稱正定的矩陣，q=(q1,q2,…,qn)T∈Rn。

引理1[9]對給定的線性互補問題(1)和標量α>0，若αI+T是非奇異的，則線性互補問題(1)與

(αI+T)x=(αI-T)|x|-q

(2)

的不動點問題等價，求x∈Rn的方程(2)。

1)若x是方程(2)的解，則

w=α(|x|-x),z=|x|+x

是問題(1)的解。

2)若向量w、z是問題(1)的解，則x=(z-α-1w)/2是不動點方程(2)的解。

若對角線元素tij=ti-j，則矩陣Tn=(tij)n×n∈Rn×n稱為Toeplitz矩陣。一個m×m塊Toeplitz矩陣

且每塊T(k)(k=0,±1,…,±(m-1))均是n×n的Toeplitz矩陣，則稱為BTTB矩陣。若每塊都是對稱的，則稱為二級對稱BTTB矩陣。若元素滿足v-k=vn-k，1≤k≤n-1，則矩陣V∈Rn×n為循環(huán)矩陣，

一個m×m塊的塊循環(huán)矩陣，若每塊均為n×n的循環(huán)矩陣，則稱為BCCB矩陣。

2 塊預(yù)處理算子Cmn

其中0≤j≤m-1，0≤k≤n-1。

方法1

1)置k=0，取x0∈Rmn。

2)共軛梯度法計算

Pmn(αImn+Tmn)xk+1=Pmnbk，

(3)

其中Pmn=Cmn(Tmn)，bk=(αImn-Tmn)|xk|-q，k=0,1,2,…。

3)計算zk+1=|xk+1|+xk+1，wk+1=Tzk+1+q。

Pmn(αImn+Tmn)x=Pmn(αImn-Tmn)|x|-q。

(4)

將式(4)減式(3)得，

Pmn(αImn+Tmn)(x*-xk+1)=

Pmn(αImn-Tmn)(|x*|-|xk|)，

等價于

(x*-xk+1)=[Pmn(αImn+Tmn)]-1×

[Pmn(αImn-Tmn)]×

(|x*|-|xk|)=(αImn+Tmn)-1×

(αImn-Tmn)(|x*|-|xk|)，

(5)

則

‖x*-xk+1‖≤‖(αImn+Tmn)-1×

(αImn-Tmn)‖‖x*-xk‖。

3 數(shù)值實驗

方法1的實驗中，利用2個方面的數(shù)據(jù)驗證實驗效果：1)方法的迭代步數(shù)；2)CPU的運行時間t。記

e(zk)=min(|(zk+1，wk+1|)，

另x0=(0,0,…,0)T∈Rnm，b=(1,1,…,1)T，迭代終止條件為e(zk)<10-7。

例1系數(shù)矩陣對角線元素[7]為

為了驗證方法1的實驗效果，方法1式(3)中令Pmn=Imn，即未預(yù)處理，與使用塊預(yù)處理算子Cmn的實驗數(shù)據(jù)進行對比。實驗結(jié)果如表1、2所示。

表1 例1在α=1時方法1的實驗結(jié)果

表2 例1在α=1.36時方法1的實驗結(jié)果

為了進一步驗證方法1的實驗效果，方法1式(3)中令Pmn=Imn，即未預(yù)處理，與使用塊預(yù)處理算子Cmn的實驗數(shù)據(jù)進行對比。實驗結(jié)果如表3、4所示。

表3 例2在α=1時方法1的實驗結(jié)果

表4 例2在α=2時方法1的實驗結(jié)果

為了進一步驗證方法1的實驗效果，方法1式(3)中令Pmn=Imn，即未預(yù)處理，與使用塊預(yù)處理算子Cmn的實驗數(shù)據(jù)進行對比。實驗結(jié)果如表5、6所示。

表5 例3在α=1時方法1的實驗結(jié)果

表6 例3在α=2時方法1的實驗結(jié)果

例4生成函數(shù)f(x)=cos3x+1，系數(shù)矩陣對角線元素[10]為

為了進一步驗證方法1的實驗效果，方法1式(3)中令Pmn=Imn，即未預(yù)處理，與使用塊預(yù)處理算子Cmn的實驗數(shù)據(jù)進行對比。實驗結(jié)果如表7、8所示。

表7 例4在α=1時方法1的實驗結(jié)果

從表1～8可看出，無論是否使用塊預(yù)處理算子，迭代步數(shù)都相同，這也可從定理1的證明過程中看出預(yù)處理矩陣的使用不影響收斂速度。但使用預(yù)處理矩陣的計算時間比未使用預(yù)處理矩陣的短，這是因為在方法1中利用了BCCB塊預(yù)處理算子，減少了方法1式(3)即內(nèi)迭代相應(yīng)的運行時間，提高了算法整體的運行效率。

表8 例4在α=1.49時方法1的實驗結(jié)果

4 結(jié)束語

為快速求解基于正定的二級對稱BTTB矩陣的線性互補問題，利用模系矩陣迭代方法，結(jié)合預(yù)優(yōu)共軛梯度法，采用文獻[8]給出的塊處理算子構(gòu)造預(yù)處理矩陣，提出了求解帶BTTB矩陣的線性互補問題的塊預(yù)處理模系矩陣迭代方法，實驗結(jié)果表明該方法是可行有效的。