康玉潔,王亞子,王朝君
(周口師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南周口466001)
線性方程組的數(shù)值解法比較與分析
康玉潔,王亞子,王朝君
(周口師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南周口466001)
討論求解線性方程組的Gauss直接消元法、Jacobi迭代法和Gauss-Seidel迭代法,并對(duì)迭代法的收斂性、收斂速度和誤差估計(jì)進(jìn)行比較,對(duì)上述兩種方法的適用范圍及優(yōu)缺點(diǎn)進(jìn)行了總結(jié).
直接法;迭代法;Gauss消去法;Jacobi迭代法;Gauss-Seidel迭代
在實(shí)際生活中,大部分科學(xué)與工程問題最終都要?dú)w結(jié)為線性方程組的求解問題.目前線性方程組求解的數(shù)值方法主要有兩類,即直接解法和迭代法.直接法特點(diǎn)是從方程組系數(shù)矩陣A和右端常數(shù)項(xiàng)b出發(fā),經(jīng)有限次的代數(shù)計(jì)算,在不考慮舍入誤差的前提下,可獲得方程組的精確解.而線性方程組的迭代法也與其他迭代法相似,首先將線性組改寫成一種迭代格式,然后再給出初始解向量,經(jīng)過迭代產(chǎn)生一組解向量的近似值,在收斂時(shí),其序列的極限向量為線性方程組的解.
直接消去法即指在不考慮舍入誤差的前提下,通過有限次迭代運(yùn)算,可獲得方程組的精確解.解線性方程組最常用的方法是Gauss消去法,即逐步消去變?cè)南禂?shù),把方程組化為等價(jià)而系數(shù)矩陣為三角形的線性方程組,再用回代過程解此等價(jià)的方程組.設(shè)n元線性代數(shù)方程組:

若記A=(aij)n×n,x=(x1,x2,…,xn)T,b=(b1,b2,…,bn)T,則上述n元線性方程組可簡記為Ax=b.
采取順序Gauss消去法[1-4],其算法如下:記

消元過程:對(duì)于k=1,2,…,n-1執(zhí)行
①若a(k)kk=0,則算法失效,停止計(jì)算;否則轉(zhuǎn)②
②對(duì)于i=k+1,k+2,…,n計(jì)算

回代過程:

由上述算法可統(tǒng)計(jì)出順序Gauss消元法求解n元線性方程組的乘除法運(yùn)算總次數(shù)為n).消元的過程與回代的過程一起組成了解線性方程組Gauss消去法的整個(gè)過程.
由于三角方程組簡單易于求解,所以考慮當(dāng)系數(shù)矩陣A的各階順序主子式均不為零時(shí),對(duì)系數(shù)矩陣采用Doolittle分解[2-6],即A=LU,其中L為下三角單位矩陣,U為上三角矩陣,且分解結(jié)果唯一.
令

這時(shí),方程組(1)就可以化為兩個(gè)容易求解的三角形方程組

先由Ly=b解出向量y,再由Ux=y(tǒng)解出向量x,這就是原方程組(1)的解向量.有分解式(3)以及矩陣相乘法,可知

當(dāng)k=2,3,…,n時(shí),有

由上述各式得到矩陣A=(aij)n×n的Doolittle分解計(jì)算公式:對(duì)于k=1,2,…,n計(jì)算

在計(jì)算機(jī)上計(jì)算時(shí),可把ukj和lik的數(shù)值分別存放在原akj和aik的存放單元內(nèi).分解計(jì)算完后,原矩陣A就成為,求解下三角方程組Ly=b和上三角方程組Ux=y(tǒng)的計(jì)算公式為

由(4)、(5)兩式,就可求解方程組(1)的解向量.
例1 用Doolittle分解求解方程組

解 運(yùn)用(4)、(5)兩式,可編寫程序?qū)栴}采取Doolittle分解進(jìn)行求解,

由Ly=b和Ux=y(tǒng)可分別解得:y=(72,90.2,59.222 2)T,x=(11,12,13)T
此時(shí)x=(11,12,13)T為方程組的精確解.
所謂迭代法,就是根據(jù)所給的線性方程組設(shè)計(jì)一個(gè)迭代公式,由迭代公式構(gòu)造一個(gè)向量序列.若此向量序列收斂,則它的收斂極限就是相應(yīng)的線性方程組的解.其中Jacobi迭代法[3-4,7-10]是基本的迭代方法,下面就考慮其具體的迭代過程和收斂性.
2.1 Jacobi迭代
對(duì)非奇異線性代數(shù)方程組Ax=b,設(shè)系數(shù)矩陣A=(aij)n×n∈Rn×n滿足條件aii≠0(i=1,2,…,n),令

由已知條件,D-1存在,那么(1)可以寫成

其中B=D-1(L1+U1),g=D-1b.
若任意給定初始向量x0=(x(0)1x(0)2…x(0)n)T,并代入(7)右端,就可以計(jì)算出一個(gè)新的向量x1,即x1=Bx0+g;再把x1代入(7)右端,又得到一個(gè)向量x2;依次類推有xk=Bxk-1+g,(k=1,2,3,…),其分量形式是

2.2 Gauss-Seidel迭代
在Jacobi迭代公式(8)的使用過程中,要同時(shí)保留兩個(gè)近似解向量xk和xk+1.如果把迭代公式改成


其中M=(D-L1)-1U1,g=(D-L1)-1b.
這樣,在整個(gè)計(jì)算過程中,只需用n個(gè)單元存貯近似解的分量.而且通常認(rèn)為,新近似解可能比老近似解更接近精確解,因此,可望這樣迭代會(huì)更有效.
選取初始向量x0,用迭代公式(9)產(chǎn)生近似解序列xk{},這種方法叫Gauss-Seidel迭代法,式(9)為Gauss-Seidel迭代法的運(yùn)算公式.
2.3 Jacobi迭代和Gauss-Seidel迭代的收斂性
通過Jacobi迭代和Gauss-Seidel迭代的過程可以看出新得出的近似解xk+1與已知近似解之間存在著線性關(guān)系,而且僅與xk有關(guān),即可表示成如下形式

其中,對(duì)Jacobi迭代:M=D-1(L1+U1),g=D-1b;對(duì)Gauss-Seidel迭代:M=(D-L1)-1U1,g=(D-L1)-1b.
定理1[3-4]解方程組(1)的單步線性迭代(11)收斂的充要條件是M的譜半徑ρ(M)<1.
定理2[3-4]若迭代矩陣M的范數(shù)‖M‖=q<1,并假定范數(shù)‖I‖=1,則線性迭代法(11)的第k次近似和準(zhǔn)確解x*之差有估計(jì)式:

從近似解的誤差估計(jì)式(12),根據(jù)事先給定的精度ε,可以計(jì)算要得到滿足精度要求的近似解需要迭代的次數(shù)k:

例2 用Jacobi迭代求解例1,并且要求各個(gè)分量的誤差絕對(duì)值不超過10-4.

則

再用Jacobi迭代法求解線性方程組.由Jacobi迭代法的計(jì)算公式(8),有

取x(0)=(0,0,0)T,代人上式,得:

如此繼續(xù)下去,運(yùn)用程序,迭代結(jié)果見表1.
在Jacobi迭代法收斂的情況下,迭代9次,得近似解x(9)=(10.999 4,11.999 4,12.999 2)T,由例1可得方程組的精確解為x*=(11,12,13)T,從表1可以看出,當(dāng)?shù)螖?shù)增加時(shí),迭代結(jié)果逐步接近精確解.

表1 Jacobi迭代所得解
近似解要求各分量的誤差絕對(duì)值不超過10-4.若對(duì)上述方程組的近似解要求各分量的誤差絕對(duì)值不超過10-4,則由

例3 用Gauss-Seidel迭代求解例1.
解 先判斷Gauss-Seidel迭代法的收斂性

M的譜半徑為ρ(M)=0.073<1,所以Gaus-Seidel迭代法收斂.
再用Gaus-Seidel迭代法求解線性方程組.由Gaus-Seidel迭代法的計(jì)算公式,有

取x(0)=(0,0,0)T,代人上式,得:如此繼續(xù)下去,運(yùn)用程序,迭代結(jié)果見表2.


表2 Gaus-Seidel迭代所得解
2.4 高階稀疏矩陣應(yīng)用迭代法求解
對(duì)于高階線性方程組,尤其是高階稀疏矩陣應(yīng)用迭代法求解一般有下面幾個(gè)步驟:
(1)構(gòu)造迭代序列;
(2)構(gòu)造的迭代序列的收斂性;
(3)如果收斂,收斂的速度.在解題時(shí)應(yīng)給出量的刻畫,用以比較各種迭代方法收斂的快慢;
(4)在有限次計(jì)算過程中,要考慮和估計(jì)近似解的誤差,以及處理迭代過程產(chǎn)生的中斷等問題,這就需要分析舍入誤差[5,10-12].
根據(jù)題中所給系數(shù)矩陣,首先判斷該迭代方法是否收斂;如果收斂,再根據(jù)求解所要求達(dá)到的精度,可以設(shè)計(jì)算法計(jì)算出要迭代的次數(shù),采用Jacobi迭代、Gauss-Seidel迭代或其他迭代方法,經(jīng)過有限步迭代,可求解出較為精確的近似解.
由以上對(duì)Jacobi迭代法和Gauss-Seidel迭代法的論述知,當(dāng)向量序列收斂時(shí),其極限就是方程組的解.當(dāng)方程組的系數(shù)矩陣為高階稀疏矩陣時(shí),運(yùn)用Jacobi迭代法即可以維持系數(shù)矩陣的稀疏特性,又容易計(jì)算和編寫計(jì)算機(jī)程序,并且在很多問題中,具有較快的收斂速度,適用于大型線性方程組,尤其是稀疏線性方程的求解問題.Gauss-Seidel迭代法對(duì)Jacobi迭代法進(jìn)行改進(jìn),在Jacobi迭代公式的使用過程中,要同時(shí)保留兩個(gè)近似解向量xk和xk+1.對(duì)于Gauss-Seidel迭代法每算出新近似解,在算下一個(gè)近似解時(shí),用新近似解代替老近似解進(jìn)行計(jì)算.這樣,在整個(gè)計(jì)算過程中,只需用n個(gè)單元存儲(chǔ)近似解的分量.而且通常認(rèn)為,新近似解可能比老近似解更靠近精確解,因此,可望這樣迭代會(huì)更有效.
在計(jì)算機(jī)存儲(chǔ)量允許的情況下,直接法是求解中小型線性方程組的有效方法,其中Gauss直接消元法作為最基本的一種直接法,其運(yùn)算量為.由此可見,對(duì)于中小規(guī)模線性方程組(即方程組的階數(shù)不高,比如小于1 000),它具有運(yùn)算量小,效率高等優(yōu)點(diǎn),但該方法用于高階方程組時(shí)計(jì)算量太大而不實(shí)用.它一般用于線性方程組的系數(shù)矩陣稠密(即系數(shù)矩陣中,大部分?jǐn)?shù)據(jù)均非零),且又沒有任何特殊結(jié)構(gòu)的線性方程組.
現(xiàn)今,隨著運(yùn)算技術(shù)的飛速發(fā)展,電子計(jì)算機(jī)的貯存量逐漸加大,所需計(jì)算時(shí)間急速減少.在計(jì)算機(jī)上,對(duì)于高階的線性方程組,可運(yùn)用直接法(如Gauss消去法)來進(jìn)行求解,但是直接法往往需要對(duì)系數(shù)矩陣A進(jìn)行Doolittle分解,因而大多數(shù)情況下丟失了系數(shù)矩陣A的稀疏性.但是在現(xiàn)實(shí)生活中,特別是求解偏微分方程的數(shù)值解時(shí),常遇到的恰恰就是大型稀疏線性方程的求解問題.
對(duì)于高階線性方程組,如果運(yùn)用直接消元法進(jìn)行求解,運(yùn)算量過大,故采用迭代法.迭代法具有存儲(chǔ)空間小,程序簡單的優(yōu)點(diǎn),所以對(duì)于大型的線性方程組,迭代法更有優(yōu)越性.線性方程組的迭代解法與直接消元法不同,其精確解不能經(jīng)過有限次的計(jì)算而獲得,而是通過迭代逐步逼近它.因此,凡是迭代解法都有一個(gè)收斂性問題,不收斂的方法是不能使用的.但迭代解法具有稀疏的系數(shù)矩陣,所需存儲(chǔ)空間少,而且程序設(shè)計(jì)簡單,適用于自主計(jì)算,用較少的運(yùn)算量,就可以得到較合理的結(jié)果.因此,對(duì)于線性方程組,特別是系數(shù)矩陣是大型稀疏的線性方程組來說,一種非常重要方法就是運(yùn)用迭代法求解.
Jacobi迭代法和Gauss-Seidel迭代法是兩種最基本、最簡單的方法,其中Jacobi迭代法具有很好的并行性,而Gauss-Seidel迭代法是典型的串行算法.對(duì)于同一個(gè)線性方程組,這兩種方法可能同時(shí)收斂,也可能同時(shí)發(fā)散,也可能其一收斂,而另一發(fā)散.但當(dāng)兩者皆收斂時(shí),一般來說Gauss-Seidel迭代法比Jacobi迭代法收斂的快.
[1]劉長安.數(shù)值分析教程[M].西安:西北工業(yè)大學(xué)出版社,2005:98-101.
[2]顏慶津.數(shù)值分析[M].北京:北京航空航天大學(xué)出版社,2006:73-75.
[3]徐樹方,高立,張平文.數(shù)值線性代數(shù)[M].北京:北京大學(xué)出版社,2004:104-105.
[4]A D Gunawardena,S K Jain,L Snyder.Modified iterative methods for consistentlin earsy stems[M].1991:123-143.
[5]曹志浩.數(shù)值線性代數(shù)[M].上海:復(fù)旦大學(xué)出版社,1996:125-127.
[6]丁麗娟,程杞元.數(shù)值計(jì)算方法(2版)[M].北京:北京理工大學(xué)出版社,2005:63-65.
[7]馮有前,王國正,李炳杰,等.數(shù)值分析[M].北京:清華大學(xué)出版社,北京交通大學(xué)出版社,2006:47-49.
[8]張耀仁.C++程序設(shè)計(jì)[M].北京:中國鐵道出版社,2006:285-294.
[9]薛秋芳.解線性方程組的幾種迭代法的收斂性分析[D].西安:陜西師范大學(xué),2005.
[10]S P Han.A lobally convergent method for nonlinear programming[J].Journal of Optimization Theory and Applications,1977(22):297-309.
[11]J V Burke,S P Han.A robust sequential quadratic programming method[J].Mathematical Programming,1989(43):277-303.
[12]M Delfour,M Fortin,G Payre.Finrte Difference Solution of a Non-linear Schrodinger[J],Equation.J.Comput.Phys.,1981(3):277-288.
Comparison and analysis of numerical methods for solving linear equations
KANG Yujie,WANG Yazi,WANG Chaojun
(School of Mathematics and Statistics,Zhoukou Normal University,Zhoukou 466001,China)
This paper researches the method for the linear algebra equations:the direct Gauss elimination method,the Jacobi iterative method and the Gauss-Seidel iterative method,compare the astringency,the convergent rate and the error estimation,then analyze the scope of application and the advantage and disadvantage of the two methods.
direct method;iterative method;Gauss elimination method;Jacobi iterative method;Gauss-Seidel iterative method
O241.6
A
1671-9476(2017)02-0036-07
10.13450/j.cnkij.zknu.2017.02.009
2016-08-09;
2016-12-22
河南省科技廳科技發(fā)展計(jì)劃項(xiàng)目(No.162102310619);河南省高等學(xué)校重點(diǎn)科研項(xiàng)目(No.17A110038);周口師范學(xué)院教育教學(xué)改革研究項(xiàng)目(No.J2016006)
康玉潔(1986-),女,河南周口人,助教,碩士,主要從事不確定理論研究.
周口師范學(xué)院學(xué)報(bào)2017年2期