鞍點問題中的位移分裂預(yù)條件技術(shù)

劉世紅，黃卓紅，蘇 翃

(1．四川工程職業(yè)技術(shù)學(xué)院基礎(chǔ)教學(xué)部，四川德陽618000; 2．重慶理工大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，重慶400054)

鞍點問題中的位移分裂預(yù)條件技術(shù)

劉世紅1，黃卓紅2，蘇 翃2

(1．四川工程職業(yè)技術(shù)學(xué)院基礎(chǔ)教學(xué)部，四川德陽618000; 2．重慶理工大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，重慶400054)

提出一類位移分裂預(yù)條件技術(shù)(SSP)，用于求解大型稀疏非正定鞍點方程組，其中該方程組的系數(shù)矩陣具有非對稱正定的(1，1)子塊，同時，對于任意迭代參數(shù)α＞0，證明這一類位移分裂迭代法是無條件收斂的，最后通過數(shù)值算例進一步驗證這類預(yù)條件技術(shù)的有效性和穩(wěn)定性．

鞍點問題;位移分裂迭代法;預(yù)處理技術(shù);Krylov子空間方法;無條件收斂

考慮如下具有2×2塊結(jié)構(gòu)的大型稀疏非對稱鞍點問題Ax=b，即

其中，B∈Rn×n是非對稱正定的()是

對稱正定的)，E∈Rn×m(n≥m)具有行滿秩(即Rank(E)=m)，u，f∈Rn以及v，g∈Rm．本文使用(·)T和(·)*分別表示某矩陣或向量的轉(zhuǎn)置和共軛轉(zhuǎn)置．

鞍點問題(1)是一類特殊的增廣線性系統(tǒng)．文獻［1］對這類問題進行了詳細(xì)的綜述，認(rèn)為這類問題廣泛應(yīng)用于科學(xué)計算和工程技術(shù)中，如帶約束條件的優(yōu)化問題、內(nèi)點問題、最小二乘問題、計算流體力學(xué)、不可壓縮性流體問題、結(jié)構(gòu)分析以及無網(wǎng)格方法等．

當(dāng)問題(1)中的系數(shù)矩陣A是大型稀疏不定矩陣時，采用直接求解方法不僅要計算(1，1)塊子矩陣B的逆還需要計算B的Schur補ETB－1E的逆，求這些矩陣逆往往需要花費昂貴的求解代價，例如求解時間以及存儲空間等．因此，迭代求解方法受到了研究人員的廣泛關(guān)注，特別是使用高效迭代算法以及新型預(yù)條件技術(shù)對Krylov子空間方法進行預(yù)處理，已經(jīng)成為求解大型稀疏鞍點問題的重要手段．

近幾十年以來，為了更有效地求解大型稀疏線性系統(tǒng)，許多高效迭代算法和新型預(yù)條件技術(shù)已經(jīng)被提出來，例如，埃爾米特和反埃爾米特分裂迭代方法(HSS)［2－3］、塊對角及塊三角預(yù)處理技術(shù)［4］、Uzawa類迭代算法［5］、超松弛迭代法(SOR)［6］、約束預(yù)處理技術(shù)［7］、位移分裂迭代法(SSP)以及廣義位移分裂迭代法(GSSP)［8－13］．

為了更有效地求解非埃爾米特正定鞍點系統(tǒng)，文獻［8］首先提出位移分裂迭代法思想，并且基于如下的矩陣分裂形式

他們提出了如下形式的位移分裂預(yù)條件子

這里正常數(shù)α為迭代參數(shù)，而In+m表示n+m階單位矩陣．

根據(jù)上述位移分裂思想，文獻［9－10］同時提出了一類具有如下形式的廣義位移分裂方法

并給出了如下形式的廣義位移分裂預(yù)條件子

其中正常數(shù)α和β為迭代參數(shù)．

特別地，當(dāng)?shù)鷧?shù)α=β時，文獻［11］提出了位移分裂迭代法，并給出了如下形式的位移分裂預(yù)條件子

文獻［12］應(yīng)用廣義位移分裂迭代法求解了大型稀疏鞍點問題(1)．根據(jù)文獻［8－13］中的理論分析，可以很容易看出，無論迭代參數(shù)α和β取什么值，當(dāng)位移迭代算法和廣義位移迭代算法被應(yīng)用于求解大規(guī)模稀疏鞍點問題時，它們都是無條件收斂的．

針對鞍點問題(1)的特殊結(jié)構(gòu)，為了獲得更有效的近似解，本文采用文獻［13］中提出的理論驗證方法，對文獻［12］中的理論推導(dǎo)進行了簡化，提出了較為簡單的驗證方法;同時，令迭代參數(shù)α=β，提出了一般的位移分裂迭代方法，通過理論分析，證明了位移分裂迭代法是無條件收斂的;最后，給出了數(shù)值算例來驗證廣義位移分裂預(yù)條件子的有效性和實用性．

1 位移分裂迭代法的收斂性分析

文獻［12］應(yīng)用廣義位移分裂迭代法求解了大型稀疏鞍點問題(1)，提出了如下形式的預(yù)條件子

并且給出了如下形式的廣義位移分裂迭代矩陣

其中，B是非對稱正定的，E具有行滿秩正常數(shù)，而α和β為正常數(shù)．通過詳細(xì)的理論分析，他們證明了這類廣義位移分裂迭代方法無條件地收斂到鞍點問題(1)的唯一解．

本文根據(jù)文獻［12］中提出的廣義位移分裂預(yù)條件子(5)，提出了如下形式的預(yù)條件子

并且給出了如下形式的位移分裂迭代矩陣

其中，B是非對稱正定的，E具有行滿秩正常數(shù)，而α為正常數(shù)．

采用文獻［13］中提出的理論驗證方法，使用一種比文獻［14］提出的更加簡單的方法來證明這類廣義位移分裂迭代法無條件地收斂到鞍點問題(1)的唯一解，并且驗證了本文提出的位移分裂迭代方法也是無條件收斂的．

記ρ(T)和 λ分別為迭代矩陣 T(α，β)(或T(α))的譜半徑和任一特征值，設(shè)(u*，v*)*為對應(yīng)于特征值λ的特征向量．根據(jù)文獻［14］中提出的分裂迭代法的收斂性理論可知，(廣義)位移分裂迭代法收斂的充要條件是ρ(T)＜1．考慮如下廣義特征值問題

方程(9)等價為

引理1．1 設(shè)B是非對稱正定的，E具有行滿秩，α和β是任意給出的正常數(shù)．另外假設(shè)λ是廣義位移分裂迭代矩陣T(α，β)的任一特征值，則λ≠±1．

證明 假設(shè)(u*，v*)*是相應(yīng)于特征值λ的特征向量．采用類似于文獻［11］中的證明方式，如果λ=1，根據(jù)(10)式，可以得到如下結(jié)論

因為系數(shù)矩陣A是非奇異的［1］，因此，u=0且v=0．然而(u*，v*)*是相應(yīng)于特征值λ的特征向量，顯然u和v不可能同時為零，從而λ≠1．

同理可證λ≠－1．從而引理得證．

定理1．1 設(shè)B是非對稱正定的，E具有行滿秩，α和β是任意給出的正常數(shù)．記ρ(T)是迭代矩陣T(α，β)的譜半徑，則ρ(τ(α，β))＜1，即廣義位移分裂迭代法無條件地收斂到鞍點問題(1)的唯一解．

證明 由引理1．1可得 λ≠ ±1．為了證明ρ(τ(α，β))＜1，這里需要考慮如下2種情形．

(i)若0≠u∈N(ET)，根據(jù)方程(10)中的第一個方程，可以得到

設(shè)u*Bu=a+bi，因為B是正定的，顯然a＞0，因此，下列不等式成立

(ii)若u?N(ET)，不失一般性，假設(shè)‖u‖2= 1，使用類似于文獻［13］中的理論驗證方法，在方程(10)的第一個及第二個方程中分別乘以向量u*及v*，并且對第二個方程進行共軛轉(zhuǎn)置，從而獲得如下結(jié)論

由方程(11)中的第二個方程可得

將方程(12)代入(11)式中的第一個方程，則有

進一步，可以計算出Ψ的實部

因此，可以獲得位移分裂迭代矩陣的特征值滿足如下不等式

這里使用Re(Ψ)和Im(Ψ)分別表示Ψ的實部及虛部，即

因此，廣義位移分裂迭代法無條件地收斂到鞍點問題(1)的唯一解．

從而，定理得證．

接下來，考慮如下形式的位移分裂迭代法．給定初始向量x0，計算

直到迭代序列xk(k=0，1，2，…)收斂到給定精度，則迭代終止．

在方程(9)和(10)中令α=β，則可以獲得如下形式的廣義特征值問題

方程(14)等價為

定理1．2 設(shè)B是非對稱正定的，E具有行滿秩，α是任意給出的正常數(shù)．記ρ(T)是迭代矩陣T(α)的譜半徑，則ρ(τ(α))＜1，即位移分裂迭代法無條件地收斂到鞍點問題(1)的唯一解．

證明 根據(jù)引理1．1，則 λ≠ ±1．為了證明ρ(τ(α))＜1，這里僅需要考慮如下2種情形．

(i)若0≠u∈N(ET)，根據(jù)方程(15)中的第一個方程有

設(shè)u*Bu=a+bi，因為B是正定的，那么a＞0，因此，下列不等式成立

(ii)若u?N(ET)，不失一般性，假設(shè)‖u‖2= 1，使用類似于文獻［13］中的理論驗證方法，在方程(15)的第一個及第二個方程中分別乘以向量u*及v*，從而獲得如下結(jié)論

對方程(16)中的第二個方程進行變形，可以獲得如下形式的等式

將方程(17)代入方程(16)中的第一個等式可得

進一步，可以計算出Ψ的實部

因此，可以獲得位移分裂迭代矩陣的特征值滿足如下不等式

這里使用Re(Ψ)和Im(Ψ)分別表示Ψ的實部及虛部．因為λ≠±1，所以ρ(τ(α))＜1．

因此，位移分裂迭代法無條件地收斂到鞍點問題(1)的唯一解．從而，定理得證．

2 數(shù)值實驗

下面將給出一個數(shù)值例子來進一步研究我們的理論分析并驗證本文提出的位移分裂預(yù)條件子的實用性及可行性．通過使用位移分裂預(yù)條件子SSP(α)及預(yù)條件子HSS預(yù)處理重啟的GMRES［15］迭代方法，將兩類預(yù)條件子SSP(α)及HSS進行了數(shù)值比較，這里HSS(α)預(yù)條件子定義為

其中

這里A是定義在問題(1)中的系數(shù)矩陣．

使用“IT”和“CPU”分別表示迭代所需迭代步數(shù)和計算時間(以秒記)．在計算中，選用零向量作為初始迭代向量，即u0=(0，0，…，0)，對于具體問題如果迭代次數(shù)超過1 000次或者第k步相對誤差滿足

則終止迭代過程．

例2．1 考慮如下Oseen方程

在這個數(shù)值實驗中，使用H．C．Elman等［16］開發(fā)的“IFISS”軟件包中的Q2－Q1有限元方法離散Oseen方程中的二維“l(fā)id－driven cavity”問題．在實際計算過程中，這里黏性系數(shù)ν=0．1，通常選擇16× 16、32×32以及64×64的四邊形一致網(wǎng)格來離散問題區(qū)域Ω．

在實際計算中，采用16×16、32×32、64×64及128×128的均勻網(wǎng)格來離散問題域Ω．表1(α= 0．01和υ=0．001)和表2(α=0．002和υ=0．01)中，對重啟數(shù)取為20的廣義極小殘量方法進行預(yù)處理，將本文提出的SSP方法和文獻［3］中的HSS方法進行了數(shù)值比較．從這2個表格中的數(shù)據(jù)可以看出:無論從迭代數(shù)還是求解時間方面，SSP方法要略優(yōu)于文獻［3］中的HSS方法．但要從理論上驗證SSP方法優(yōu)于HSS方法以及確定SSP預(yù)條件子的最優(yōu)迭代參數(shù)仍需要進一步研究．

表1 迭代步數(shù)及計算時間Table 1 Number of iterations and solution time in seconds

表2 迭代步數(shù)及計算時間Table 2 Number of iterations and solution time in seconds

［1］BENZI M，GOLUB G H，LIESEN J．Numerical solution of saddle point problems［J］．Acta Numer，2005，14:1－137．

［2］BAI Z Z．Optimal parameters in the HSS－like methods for saddle－point problems［J］．Numer Linear Algebra Appl，2009，16: 447－479．

［3］BENZI M，GOLUB G H．A preconditioner for generalized saddle point problems［J］．SIAM J Matrix Anal Appl，2004，26:20－41．

［4］BAI Z Z，CHAN R H，REN Z R．On order－reducible sinc discretizations and block－diagonal preconditioning methods for linear third－order differential equations［J］．Numer Linear Algebra Appl，2014，21:108－135．

［5］CHEN P Y，HUANG J G，SHENG H S．Some Uzawa methods for steady incompressible Navier－Stokes equations discretized by mixed element methods［J］．J Comput Appl Math，2015，273:313－325．

［6］雷剛．預(yù)處理(I+S)后改進的SOR迭代法收斂性討論［J］．四川師范大學(xué)學(xué)報(自然科學(xué)版)，2011，34(4):528－531．

［7］BERGAMASCHI L．On eigenvalue distribution of constraint－preconditioned symmetric saddle point matrices［J］．Numer Linear Algebra Appl，2012，19:754－772．

［8］BAI Z Z，YIN J F，SU Y F．Shift－splitting preconditioner for non－Hermitian positive definite matrices［J］．J Comput Math，2006，24:539－552．

［9］曹陽，陶懷仁．鞍點問題的廣義位移分裂預(yù)條件子［J］．計算數(shù)學(xué)，2014，36:16－26．

［10］CHEN C R，MA C F．A generalized shift－splitting preconditioner for saddle point problems［J］．Appl Math Lett，2015，43:49－55．

［11］CAO Y，DU J，NIU Q．Shift－splitting preconditioners for saddle point problems［J］．J Comput Appl Math，2014，272:239－250．

［12］CAO Y，LI S，YAO L Q．A class of generalized shift－splitting preconditioners for nonsymmetric saddle point problems［J］．Appl Math Lett，2015，49:20－27．

［13］SALKUYEH D K，MASOUDI M，HEZARI D．On the generalized shift－splitting preconditioner for saddle point problems［J］．Appl Math Lett，2015，48:55－61．

［14］YOUNG D M．Iterative Solution of Large Linear Systems［M］．New York:Academic Press，1971．

［15］SAAD Y，SCHULTZ M H．GMRES:a generalized minimal residual algorithm for solving nonsymmetric linear systems［J］．SIAM J Sci Statist Comput，1986，7:856－869．

［16］ELMAN H C，RAMAGE A，SILVESTER D J．IFISS:a Matlab toolbox for modelling incompressible flow［J］．ACM Trans Math Software，2007，33:Article14．

［17］于善玲，張耀明．二維Helmholtz方程的邊界元法［J］．重慶理工大學(xué)學(xué)報(自然科學(xué)版)，2015，29(11):139－143．

On the Shift-splitting Preconditioning Technique for the Saddle Point Problems

LIU Shihong1，HUANG Zhuohong2，SU Hong2

(1．Department of General Studies，Sichuan Engineering Technical College，Deyang 618000，Sichuan; 2．School of Mathematics and Statistics，Chongqing University of Technology，Chongqing 400054)

In this paper，the shift-splitting iteration method is used to solve the large sparse indefinite saddle point systems with nonsymmetric positive definite(1，1)-block．It is analysed that the shift-splitting iteration method converges unconditionally to a unique solution of saddle point system for any iteration parameter α＞0．A numerical experiment shows the correctness and the feasibility of the shift-splitting iteration method．

saddle point problems;shift-splitting iteration method;preconditioning technique;Krylov subspace methods;unconditional convergence

O151．21

A

1001－8395(2016)04－0549－05

10．3969/j．issn．1001－8395．2016．04．016

(編輯 李德華)

2015－11－25

重慶市基礎(chǔ)與前沿研究一般項目(CSTC2015JCYJAA1432)和重慶市教委科技項目(KJ1500941)

劉世紅(1960—)，男，講師，主要從事代數(shù)學(xué)、數(shù)值代數(shù)及其應(yīng)用的研究，E－mail:lsh@scetc．edu．cn

2010 MSC:65F10