于春肖， 楊艷芳， 井丁卉

(燕山大學(xué) 理學(xué)院 河北 秦皇島 066004)

0 引言

在科學(xué)和工程領(lǐng)域中,很多問題最終都可歸結(jié)為大型線性方程組的求解,文獻(xiàn)[1]提出的基于Galerkin原理的GMRES(m)算法是公認(rèn)的求解大型非對稱線性方程組較為成功的方法之一,同時(shí)也是快速多極邊界元法中的求解器。

隨著GMRES算法和相關(guān)知識(shí)被廣泛學(xué)習(xí),很多學(xué)者也對該算法進(jìn)行了優(yōu)化和改進(jìn)。一種方法是采用預(yù)處理技術(shù)[2-4]。劉曉明等[5]將GMRES方法應(yīng)用于油藏?cái)?shù)值模擬計(jì)算問題,利用矩陣分塊技術(shù),采用塊擬消去法(PE)對系數(shù)矩陣進(jìn)行了預(yù)處理。于春肖等[6]提出ILU預(yù)處理Householder-GMRES(m)算法,通過降低系數(shù)矩陣的條件數(shù)達(dá)到加快收斂的目的。關(guān)朋燕等[7]將多項(xiàng)式預(yù)處理與加權(quán)GMRES(m)相結(jié)合,構(gòu)造出一種新的方法,減少了計(jì)算量。于春肖等[8]通過改變重啟動(dòng)參數(shù)m的值提出VRP-HGMRES(m)算法,該算法有效地避免了因重啟動(dòng)參數(shù)的選取不當(dāng)而造成的問題。 文獻(xiàn)[9]將稀疏矩陣的小特征值對應(yīng)的特征向量加入Krylov子空間。于春肖等[10]提出一種截?cái)嘈虸GMRES(m)新算法,此算法是基于FMM的Krylov子空間投影類方法。郝雪景等[11]提出基于不完全正交的變參數(shù)IGMRES(m)算法。Walker[12]對GMRES算法進(jìn)行改進(jìn),提出一種基于Householder 變換的GMRES(簡稱H-GMRES)算法。

本文將利用截?cái)嗉夹g(shù)[13-14]給出一種不完全正交的變參數(shù)H-IGMRES(m)算法,使計(jì)算量大大減少,并結(jié)合數(shù)值算例分析該算法的高效性。

1 H-GMRES(m)算法

在實(shí)際計(jì)算中,H-GMRES(m)算法構(gòu)造Hessenberg矩陣元素和Krylov向量時(shí)都要用到較長的遞推式,從而導(dǎo)致消耗的計(jì)算時(shí)間很長,迭代過程數(shù)值也不穩(wěn)定。這種算法的具體步驟為:

1) 給定重啟動(dòng)參數(shù)m和誤差上界ε>0,取初始值x(0)∈Rn,計(jì)算r(0)=b-Ax(0),β=‖r(0)‖2。

2) 計(jì)算P1=I-2ωωΤ,其中：‖ω‖2=1；P1r(0)=βe1；v1=P1e1。

3) 對于j=1,2,…,m，令cj=PΤAvj(當(dāng)j=1時(shí),有P=P1成立),然后計(jì)算

vj+1=Pj+1Pj…P1ej+1，hj=Pj+1cj。

4) 如果‖r(0)-Azm‖2<ε，則x=x(m),迭代停止；否則x(0)=x(m),m=m+1且轉(zhuǎn)向1)再一次迭代,直到滿足給定的誤差上界為止。

2 H-IGMRES(m)算法與收斂性分析

2.1 H-IGMRES(m)算法

1) 初始化:給定重啟動(dòng)參數(shù)m和誤差上界ε>0,設(shè)置截?cái)嘀笜?biāo)q0,取初始值x(0)∈Rn,計(jì)算

r(0)=b-Ax(0),β=‖r(0)‖2。

2) 計(jì)算P1=I-2ωωΤ,其中‖ω‖2=1；P1r(0)=βe1；v1=P1e1。

3) 對于j=1,2,…,m有不完全正交化,令cj=PΤAvj(當(dāng)j=1時(shí),有P=P1成立),然后計(jì)算

vj+1=Pej+1，hj=Pj+1cj。

求解最小二乘問題

得出原方程組的近似解x(m)=x(0)+zm=x(0)+Vmym。

4) 如果‖r(0)-Azm‖2<ε，則x=x(m),迭代停止；否則x(0)=x(m),m=m+1且轉(zhuǎn)向1)再一次迭代,直到滿足給定的誤差上界為止。

2.2 收斂性分析

為研究這種截?cái)嘈虷-IGMRES(m)算法的收斂形態(tài),這里將H-IGMRES(m)算法和H-GMRES(m)算法結(jié)合起來,以建立H-IGMRES(m)算法的收斂性理論。

定義Vm+1=(v1,v2,…,vm+1),可以證明,如果H-IGMRES(m)算法在m步中斷且Vm+1列滿秩,則得到的數(shù)值解為準(zhǔn)確解x*。另外,當(dāng)q0=1時(shí),H-IGMRES(m)算法與H-GMRES(m)算法相同。為討論方便,選取Km=span{r(0),Ar(0),…,Am-1r(0)},將H-IGMRES(m)算法進(jìn)行m步迭代后的數(shù)值解表示為x(m)(HI),殘余向量表示為r(m)(HI); 將H-GMRES(m)算法進(jìn)行m步迭代后的數(shù)值解表示為x(m)(H),殘余向量表示為r(m)(H)。

定理1假設(shè)Vm+1列滿秩，則H-IGMRES(m)算法和H-GMRES(m)算法在Km上的殘值范數(shù)滿足

‖r(m)(HI)‖≤S(Vm+1)‖r(m)(H)‖。

(1)

其中Vm+1由不完全正交化過程生成,S(Vm+1)=‖Vm+1‖‖Vm+1+‖,Vm+1+是Vm+1的廣義逆矩陣。如果H-IGMRES(m)算法在m步中斷,即hm+1,m=0,則x(m)(HI)=x*。

因此

(2)

因此,式(2)變成

‖r(m)(HI)‖=‖Vm+1(Vm+1+r(0)-Vm+1+AVmym)‖。

因此可得

所以式(1)成立。

3 數(shù)值算例

圖1 精確解與數(shù)值解的比較Figure 1 Comparison of exact solution and numerical solution

Ay=b,A∈R(n-1)2×(n-1)2,y,b∈R(n-1)2。

如果選取p=0.01,q=0.5,n=40,初始向量為y(0)=(0,…,0)Τ,誤差上界為t=1e-6,當(dāng)m=10時(shí),取截?cái)嘀笜?biāo)為9,用H-IGMRES(m)算法求解該線性方程組,將得到的數(shù)值解與精確解進(jìn)行比較,如圖1所示。從圖1可看出,數(shù)值解和精確解完全吻合,證明此算法的可行性。

下面選取n=100,分別用H-IGMRES(m)算法和H-GMRES(m)算法求解這個(gè)線性方程組。當(dāng)重啟動(dòng)參數(shù)不同時(shí),兩種算法所需的迭代次數(shù)和得到的殘值范數(shù)如圖2所示,兩種算法所需的運(yùn)行時(shí)間見表1。從圖2和表1可看出,重啟動(dòng)參數(shù)不同時(shí),H-IGMRES(m)算法與H-GMRES(m)算法相比,前者迭代次數(shù)和運(yùn)行時(shí)間較少,說明新算法提高了計(jì)算效率,并且具有更高的計(jì)算精度。

圖2 兩種算法計(jì)算效率與計(jì)算精度的比較Figure 2 Comparison of computational efficiency and accuracy between two algorithms

表1 當(dāng)重啟動(dòng)參數(shù)不同時(shí),兩種算法運(yùn)行時(shí)間比較(m=10)Table 1 Comparison of running time under different meshes for the two algorithms(m=10)

例2考慮如下二維Possion方程,其中Ω{(x,y)|0≤x≤1,0≤y≤1}，u(x,y)為溫度分布函數(shù)。

令f(x,y)=2π2sin πxsin πy,將求解區(qū)域沿x方向和y方向進(jìn)行等分,由五點(diǎn)差分法可得Au=b,A∈Rn2×n2,u,b∈Rn2×1。

當(dāng)n=35時(shí),選取重啟動(dòng)參數(shù)m=20及截?cái)嘀笜?biāo)q0=9,分別用高斯消元法和H-IGMRES(m)算法求解線性方程組,其結(jié)果分別如圖3(a)和3(b)所示。

圖3 溫度分布Figure 3 Temperature distribution

由圖3可以看出,兩種算法的計(jì)算結(jié)果完全一致,這表明H-IGMRES(m)算法是正確的。截?cái)嘀笜?biāo)可以在1～19中取值,不同的截?cái)嘀笜?biāo)對計(jì)算速率、計(jì)算精度、運(yùn)行時(shí)間的影響如圖4所示。

圖4 截?cái)嘀笜?biāo)對計(jì)算速率、計(jì)算精度、運(yùn)行時(shí)間的影響Figure 4 Effect of truncation index on computional rate, computional accuracy, operation time

由圖4可看出,迭代次數(shù)和殘值范數(shù)隨截?cái)嘀笜?biāo)的增大呈下降趨勢,運(yùn)行時(shí)間也在突降之后逐漸趨于平穩(wěn),說明截?cái)嗉夹g(shù)在保證計(jì)算精度的前提下提高了計(jì)算效率。當(dāng)m=10時(shí),隨著計(jì)算規(guī)模n的不斷增大,當(dāng)殘值范數(shù)滿足給定的誤差上界時(shí),H-IGMRES(m)算法總的運(yùn)行時(shí)間要比H-GMRES(m)算法少很多,進(jìn)一步說明了H-IGMRES(m)算法的優(yōu)越性,見表2。

表2 不同網(wǎng)格下兩種算法運(yùn)行時(shí)間比較(m=10)Table 2 Comparison of running time under different meshes for the two algorithms(m=10)

4 結(jié)論

1) 提出了H-IGMRES(m)算法,對算法的收斂性進(jìn)行了理論分析。

2) 數(shù)值算例表明H-IGMRES(m)算法的計(jì)算結(jié)果和精確解吻合,說明該算法的有效性。

3) 分析了截?cái)嘀笜?biāo)的選取對算法的計(jì)算精度和效率的影響,證明了H-IGMRES(m)算法的計(jì)算效率和精度都高于H-GMRES(m)算法。