第5章 線性方程組的數(shù)值解法

在工程技術(shù)、自然科學和社會科學中，經(jīng)常遇到的許多問題最終都可歸結(jié)為解線性方程組。◎電學中的網(wǎng)絡問題、◎最小二乘法求實驗數(shù)據(jù)的曲線擬合問題、◎三次樣條函數(shù)的插值問題、◎經(jīng)濟運行中的投入產(chǎn)出問題、◎大地測量、機械與建筑結(jié)構(gòu)的設計計算問題等等，歸結(jié)為求解線性方程組或非線性方程組的數(shù)學問題。因此，掌握線性方程組的解法對于解決實際問題是極其重要的。第5章線性方程組的數(shù)值解法

1簡記為Ax=b，其中：

一般地，設b≠0,當系數(shù)矩陣A非奇異(即detA≠0)時，如下方程組有惟一解。2線性方程組的數(shù)值解法一般有兩類：直接法：就是經(jīng)過有限步算術(shù)運算，可求得方程組精確解的方法（若計算過程中沒有舍入誤差）。直接法中具有代表性的算法是高斯(Gauss)消去法。（低階稠密矩陣）迭代法:

就是用某種極限過程去逐步逼近線性方程組的精確解的方法。也就是從解的某個近似值出發(fā)，通過構(gòu)造一個無窮序列去逼近精確解的方法。(一般有限步內(nèi)得不到精確解)。（高階稀疏矩陣）解線性方程組有Gram法則，但Gram法則不能用于計算方程組的解，如n＝100，1033次/秒的計算機要算10120年.3三角形方程組的解乘除法：（n＋1）n/2次運算(1)上三角方程的一般形式

5.1高斯消去法

解：

(2)下三角方程的一般形式解：4消去后兩個方程中的x1得再消去最后一個方程的x2得消元結(jié)束,經(jīng)過回代得解:Gauss消去法5上述求解的消元過程可用矩陣表示為:(A,b)=這是Gauss消去法的計算形式,新的增廣矩陣對應的線性方程組就是上三角形方程組,可進行回代求解?，F(xiàn)在介紹求解一般線性方程組的順序Gauss消去法:記則,線性方程組的增廣矩陣為6第一步.設,依次用乘矩陣的第1行加到第i行,得到矩陣:其中7第二步.設,依次用乘矩陣的第2行加到第i行,得到矩陣:其中8如此繼續(xù)消元下去,第n-1步結(jié)束后得到矩陣:這就完成了消元過程。對應的方程組變成：對此方程組進行回代，就可求出方程組的解。9順序Gauss消去法求解n元線性方程組的乘除運算量是：

n2-1+(n-1)2-1+…+22-1

+1+2+…+n

n=20時,順序Gauss消去法只需3060次乘除法運算.順序Gauss消去法通常也簡稱為Gauss消去法.順序Gauss消去法中的稱為主元素.主元素都不為零

Gauss消去法才能執(zhí)行下去.

主元素都不為零矩陣A的各階順序主子式都不為零.

10步1

輸入系數(shù)矩陣A，右端項b，置k:=1；步2.

消元：對k=1,2,…,n-1，計算步3.

回代對k=n-1,…,1，計算算法5.1順序Guass消去法

11

列主元Gauss消元法解:[方法1]小主元a11回代求解x2=1.0000,x1=0.0000.計算結(jié)果相當糟糕.

原因:求行乘數(shù)時用了”小主元”0.0001作除數(shù).例1

在F(10,4,L,U)上用Gauss消去法求解線性方程組方程組的準確解為x*=(10000/9999,9998/9999)T.12解:[方法2]若上題在求解時將1,2行交換,即回代后得到

x=(1.0000,1.0000)T與準確解x*=(9998/9999,10000/9999)T相差不大.主元a11

方法2所用的是在Gauss消去法中加入選主元過程,利用換行,避免”小主元”作除數(shù),該方法稱為Gauss列主元消去法。13由于計算機和方程組本身的限制，在Gauss消元過程中會出現(xiàn)零主元和小主元，而使Gauss消去法無法進行或出現(xiàn)較大舍入誤差，為避免此類現(xiàn)象，可以通過行交換來選取絕對值較大的主元素.常采用的是列主元消去法和全主元消去法.給定線性方程組Ax=b,記A(1)=A,b(1)=b,列主元Gauss消去法的具體過程如下:首先在增廣矩陣B(1)=(A(1),b(1))的第一列元素中,取然后進行第一步消元得增廣矩陣B(2)=(A(2),b(2)).14Gauss列主元消去法因為|mi,k|≤1,列主元消去法能避免”小主元”作除數(shù),誤差分析與數(shù)值試驗表明：Gauss列主元消去法“基本穩(wěn)定”.再在矩陣B(2)=(A(2),b(2))的第二列元素中,取按此方法繼續(xù)進行下去,經(jīng)過n-1步選主元和消元運算,得到增廣矩陣B(n)=(A(n),b(n)).則方程組A(n)x=b(n)是與原方程組等價的上三角形方程組,可進行回代求解.然后進行第二步消元得增廣矩陣B(3)=(A(3),b(3)).易證,只要|A|0,列主元Gauss消去法就可順利進行.15算法3.2（列主元Gauss消去法）：步1、輸入系數(shù)矩陣A，右端項b，置k=1;步2、對于(1)按列選主元：選取l

使

(2)如果，交換(A(k),b(k))

的第k行與第l

行元素(3)

消元計算：3、回代計算16Gauss-Jordan列主元消去法

高斯消去法有消元和回代兩個過程，消去的是對角線下方的元素。對消元過程稍加改變便可使方程組化為對角陣。這時求解就不需要回代。這種將主元素化為1,并用主元將其所在列的冗余元素全都消為0,即消去對角線上方與下方元素的方法稱為Gauss-Jordan消去法。這時等號右端即為方程組的解。17例2用Gauss-Jordan消去法求方程組的解解方程組相應的增廣矩陣為：

列選主元故得：

18定理5.1設A為非奇異矩陣,方程組Ax=I的增廣矩陣為C=A

,I]

，如果對C應用高斯-約當消去法化為I

,T，則T

=A-1。的逆矩陣.

例3用Gauss-Jordan消去法求

解C=A

I

=

195.2矩陣的三角分解

順序Gauss消去法的矩陣表示若令記則有A(2)=對于給定方陣A，如果能求出下三角方陣L和上三角方陣U，使A=LU，就說對方陣A進行了三角分解，三角分解也稱LU分解。20記令若則有A(3)=如此進行下去,第n-1步得到:

A(n)=其中也就是:

A(n)=Ln-1A(n-1)=Ln-1Ln-2A(n-2)=…=Ln-1Ln-2…L2L1A(1)21所以有:A=A(1)=L1-1L2-1…Ln-1-1A(n)=LU其中L=L1-1L2-1…Ln-1-1.而且有式A=LU稱為矩陣A的三角分解.A(n)=Ln-1Ln-2…L2L1A(1)令A(n)=U，即22

設A為n階方陣，如果解Ax=b用Gauss消去法能夠完成(即akk(k)不為零),則矩陣A可分解為單位下三角矩陣L和上三角矩陣U的乘積,且這種分解是唯一的。

證明只證唯一性,設有兩種分解A=LU=L1U1,則有=I所以得定理5.2.1定理5.2.2約化的主元素

的充要條件是矩陣A的順序主子式于是Ax=b

LUx=b

令

Ux=y

得23LU分解的直接推導直接利用矩陣乘法來計算

LU分解的元素比較等式兩邊的第一行得：u1j=a1j比較等式兩邊的第一列得：比較等式兩邊的第二行得：比較等式兩邊的第二列得：(j=1,…,n)(i=2,…,n)(j=2,…,n)(i=3,…,n)U

的第一行L

的第一列U

的第二行L

的第二列24LU分解算法（續(xù)）第k

步：此時U

的前k-1行和

L

的前k-1列已經(jīng)求出比較等式兩邊的第k行得：比較等式兩邊的第k列得：直到第n

步，便可求出矩陣L

和U

的所有元素。(j=k,…,n

)(i=k+1,…,n

)25LU分解算法算法5.3：(LU

分解

)For

k=1,2,...,nEndForj=k,…,ni=k+1,…,n為了節(jié)省存儲空間，通常用A

的絕對下三角部分來存放L(對角線元素無需存儲)，用

A

的上三角部分來存放U。運算量：(n3-n)/3Doolittle三角分解26由可得這就是求解方程組Ax=b的Doolittle三角分解方法。27例利用三角分解方法解線性方程組

解因為所以先解

,得28,得解線性方程組Ax=b的Doolittle三角分解法的計算量約為n3/3,與Gauss消去法基本相同.其優(yōu)點在于求一系列同系數(shù)的線性方程組Ax=bk,(k=1,2,…,m)時,可大大節(jié)省運算量.此外，還有另一種LU分解，稱為Crout三角分解.再解29

A正定：A的各階順序主子式均大于零。對稱正定矩陣的三角分解如果A為對稱正定矩陣，則存在一個實的非奇異下三角矩陣L，使A=LLT

，且當限定的對角元素為正時，這種分解是唯一的。三、平方根法A是正定矩陣，即對于任意n維非零列向量x

，恒有xTAx

>0，A對稱：AT=A

定理5.2.3

若A為對稱正定矩陣，則A可唯一分解為A=LDLT。其中L為單位下三角，D為元素全為正數(shù)的對角陣。30U=uij=u11uij/uii111u22unn記為

A對稱且分解是唯一的即記D1/2=則仍是下三角陣設A為對稱正定矩陣,則有唯一分解A=LU,且ujj>0.分解A=LLT稱為對稱正定矩陣的Cholesky分解.如何求L？31一般有(第j列)稱為平方根法，因為帶了開方運算。對于k=1,2,…,n計算（一列一列算），k=1時32Cholesky分解法特點：

2.對于對稱正定陣A，從可知對任意ki

有。即L

的元素不會增大，誤差可控，不需選主元。缺點：存在開方運算。優(yōu)點：1.可以減少存儲單元。為一般矩陣LU分解計算量的一半。For

j=1,2,...,n算法5.4：(Cholesky

分解，按列)EndFori=k+1,…,n運算量：n3/6

+n2/2+n/3

33改進Cholesky分解法改進的cholesky分解法,為了避免開方運算，分解A=LDLT現(xiàn)確定計算L的第j列的元素。考察A的元素aij，ajj34算法5.5：(改進的Cholesky

分解法

)步1

輸入對稱正定矩陣A，右端項b；步2.

改進的Cholesky分解：對j=2,3,…,n，計算：步3.解下三角形方程組LDy=b步4.解上三角形方程組LTx=y改進后的LDLT分解乘除法運算量約為n3/4+O(n2),但存儲量增加了。35追趕法Ax=f其中：36用矩陣的Doolitte分解，有得追：趕：計算量為5n-4次乘除法追趕法是數(shù)值穩(wěn)定的,追趕法具有計算程序簡單,存貯量少,計算量小的優(yōu)點.37例

用追趕法求解三對角線性方程組Ax=f，其中解按公式分解有385.3向量和矩陣的范數(shù)

1．向量的范數(shù)定義1：設x

R

n，x表示定義在Rn上的一個實值函數(shù)，稱之為x的范數(shù)，它滿足下列性質(zhì)：（3）三角不等式：即對任意兩個向量x、y

R

n，恒有(1)非負性：即對一切x

R

n，x

0,x>0(2)齊次性：即對任何實數(shù)a

R,有（4）||

-x||=|

-1

|

×

||

x||=||

x||;（5）對任意的x、y?Rn，恒有

|

||

x||-||

y||

|≤

||

x-y||.有性質(zhì)：39設x=(x1,x2,…,xn)T

三個常用的范數(shù)：

在Rn上定義的任一向量范數(shù)

都與范數(shù)

等價，成立。結(jié)論：Rn上定義的任何兩個范數(shù)都是等價的。

對常用范數(shù)，容易驗證下列不等式：

即存在正數(shù)M與m(M>m)對一切xRn，不等式40定義2：設給定Rn中的向量序列{}，即若對任何i(i=1,2,…,n)都有則向量

或者說向量序列{}依坐標收斂于向量，記為，稱為向量序列{

}的極限，定理：向量序列{x(k)}依坐標收斂于x*的充要條件是向量序列依范數(shù)收斂與依坐標收斂是等價的。依范數(shù)收斂412矩陣的范數(shù)

定義3設‖‖是以n階方陣為變量的實值函數(shù),且滿足條件:

(1)非負性:‖A‖0,且‖A‖=0當且僅當A=0(2)齊次性:‖kA‖=|k|‖A‖,kR

(3)三角不等式:‖A+B‖‖A‖+‖B‖(4)相容性:‖AB‖‖A‖‖B‖則稱‖A‖為矩陣A的范數(shù).

記A=(aij),常用的矩陣范數(shù)有:

矩陣的1-范數(shù):‖A‖1

,也稱矩陣的列范數(shù).

矩陣的2-范數(shù):‖A‖2

,也稱為譜范數(shù).（4'）對任意向量xRn，和任意矩陣A，有

矩陣的-范數(shù):‖A‖

,也稱為矩陣的行范數(shù).42例設矩陣求矩陣A的范數(shù)‖A‖p,p=1,2,.

解

‖A‖1=

‖A‖=5,4，43在有限維線性空間上任何兩種矩陣范數(shù)是等價的。即

m

‖A‖‖A‖

M

‖A‖,ARnn對常用范數(shù)，容易驗證下列不等式都成立矩陣的范數(shù)與矩陣的特征值之間也有密切的聯(lián)系.設是矩陣A的特征值,

x是對應的特征向量,則有

Ax=x利用向量和矩陣范數(shù)的相容性,則得

||‖x‖=‖x‖于是||‖A‖設n階矩陣A的n個特征值為1,2,…,n,則稱為矩陣A的譜半徑.對矩陣的任何一種相容范數(shù)都有

(A)‖A‖=‖Ax‖‖A‖‖x‖44

結(jié)論：定義5.3.4設有Rnn中的矩陣序列{A(k)},其中A(k)=(aij(k))nn,若對任何i

和j(i,j=1,2,…,n)，都有則矩陣A*=(aij*)nn,稱為矩陣序列{A(k)}的極限，記為依范數(shù)收斂455.4線性方程組的性態(tài)和解的誤差分析方程組的精確解為x=(2,0)T.設有線性方程組Ax=b,A為非奇異方陣，x為精確解。由于A和b元素是測量或觀測得到的，常帶有某些觀測誤差，或A,b的元素是計算的結(jié)果，因而包含有舍入誤差。

實際處理矩陣是A+δA，b+δb，其中δA是矩陣A的微小誤差矩陣或擾動；δb是向量b的微小誤差或擾動。需要研究方程組數(shù)據(jù)A,b的微小誤差（擾動）對解x的影響，即考慮方程組的解y和x的差的估計問題。例5.4.146考察常數(shù)項有微小變動則方程組為即.常數(shù)項相對誤差

,由于常數(shù)項微小誤差，而引起解的相對誤差較大，為常數(shù)項相對誤差的104倍，也就是說，此方程組的解對方程組的數(shù)據(jù)A,b非常敏感。這種由于原始數(shù)據(jù)微小變化而導致解嚴重失真的線性方程組稱為病態(tài)方程組,相應的系數(shù)矩陣稱為病態(tài)矩陣.而引起解的相對誤差為47設一般線性方程組Ax=b1)系數(shù)矩陣是精確的,常數(shù)項有誤差δb,此時記解為x+δx

,則

A(x+δx)=b+δb于是

Aδx=δb所以又由于因此即‖A-1‖‖δb‖‖δx

‖=‖A-1δb‖‖A‖‖x‖‖b‖=‖Ax‖‖A‖‖A-1‖‖δb‖‖x‖‖δx

‖‖b‖48

2)再設b是精確的,A有誤差δA,此時記解為x+δx

,則

(A+δA)(x+δx)=b則有所以

δx

=-A-1δA(x+δx)于是也就是記Cond(A)=‖A‖‖A-1‖,稱為矩陣A的條件數(shù).

Aδx

+δA(x+δx)=0

‖δx‖‖A-1‖‖δA‖‖x+δx‖49條件數(shù)的性質(zhì)：

1）cond(A)≥12）cond(kA)=cond(A)(k為非零常數(shù))3）若，則經(jīng)常使用的條件數(shù)有：Condp(A)=‖A‖p‖A-1‖pp=1,2,。記Cond(A)=‖A‖‖A-1‖,稱為矩陣A的條件數(shù).是良態(tài)方程組，或A是良態(tài)的。例5.4.2

對于例5.4.1的方程組，計算A的條件數(shù)。.50由于計算條件數(shù)運算量較大,實際計算中若遇到下述情況之一,方程組就有可能是病態(tài)的:行列式很大或很小（如某些行、列近似相關）；元素間相差大數(shù)量級，且無規(guī)則；主元消去過程中出現(xiàn)小主元；特征值相差大數(shù)量級。Hilbert矩陣病態(tài)方程組的計算1)用雙精度或更高精度計算2)采用數(shù)值穩(wěn)定性較好的算法，如全主元法3)用預處理方法。515.5

解線性方程組的迭代法迭代法的基本思想是構(gòu)造一串收斂到解的序列，即建立一種從已有近似解計算新的近似解的規(guī)則。思路將Ax=b等價改寫為x=Bx+f,建立迭代格式：x(k+1)=Bx(k)+f,從初值x(0)出發(fā)，得到迭代序列{x(k)}(k=0,1,...).如果其迭代序列的極限存在,即則稱迭代法收斂,

x*即是線性方程組Ax=b的解.否則，稱迭代法發(fā)散.矩陣B稱為迭代矩陣，x(k)為第k次迭代近似解，稱e(k)=x(k)-x*為第k次迭代誤差向量.研究內(nèi)容：如何建立迭代格式？

迭代序列收斂速度？

迭代序列的收斂條件？

近似解的誤差估計？525.5.1Jacobi迭代法若系數(shù)矩陣非奇異，且

(i=1,2,…,n)，將方程組改成53然后寫成迭代格式（5.5.5）（5.5.6）式也可以簡單地寫為Jacobi迭代格式。54寫成矩陣形式：A=－L－UDBJacobi迭代陣(5.5.8)Jacobi迭代格式簡單，每迭代一次只需計算一次矩陣與向量乘法，計算機中只需要兩組工作單元用來保存x(k)及x(k+1),用來控制迭代終止。55算法5.7（Jacobi迭代法）1.輸入矩陣A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次數(shù)N，置k=0;

2.對i=1,2,…,n,計算3.若‖x(k+1)-x(k)‖<ε,則停算，輸出x(k+1)為近似解；4.若k<N,置k=k+1,轉(zhuǎn)步2；否則，停算，輸出迭代失敗信息.

56

例5.5.1

用雅可比迭代法解方程組解：雅可比迭代格式為取經(jīng)14次迭代,近似解為誤差為57為了加快收斂速度，同時為了節(jié)省計算機的內(nèi)存，作如下的改進：

5.5.2Gauss-Seidel迭代法逐一寫出來即為每算出一個分量的近似值，立即用到下一個分量的計算中去，即用迭代格式：58…………寫成矩陣形式：BGauss-Seidel

迭代陣(5.5.11)(5.5.13)只保存一組向量59算法5.8（Gauss-Seidel迭代法）：1.輸入矩陣A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次數(shù)N，置k=0;

2.計算對i=2,…,n-1,計算3.若‖x(k+1)-x(k)‖<ε,則停算，輸出x(k+1)為近似解；4.若k<N,置k=k+1,轉(zhuǎn)步2；否則，停算，輸出迭代失敗信息.

60

例用Gauss-Seidel迭代法解方程組取經(jīng)6次迭代,近似解為誤差為解：Gauss-Seidel迭代格式61松弛（5.5.18）是松馳因子(0<<2),當0<<1時叫低松弛，>1時叫超松弛，=1時，就是Gauss-Seidel迭代法。5.3

超松弛法逐次超松弛迭代法（SuccessiveOverRelaxationMethod.簡稱為SOR法）是Gauss-Seidel法的一種加速方法。這種方法將前一步的迭代結(jié)果xi(k)與Gauss-Seidel法的迭代值,

適當進行線性組合，以構(gòu)成一個收斂速度較快的近似解序列。62對Gauss－Seidel迭代格式有故SOR的迭代格式用矩陣語言描述為：Bωfω63算法5.9（SOR迭代法）1.輸入矩陣A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次數(shù)N，參數(shù)ω；置k=0;

2.計算對i=2,…,n-1,計算3.若‖x(k+1)-x(k)‖<ε,則停算，輸出x(k+1)為近似解；4.若k<N,置k=k+1,轉(zhuǎn)步2；否則，停算，輸出迭代失敗信息.

64例5.5.3用SOR法法解方程組初值x(0)=0，寫出計算格式。解：SOR迭代格式為

選取ω=1.3,第11次迭代結(jié)果為：

且滿足:65記誤差向量e(k)=x(k)-x*,則迭代法收斂就是e(k)0.

x(k+1)=Bx(k)+f

k=0,1,2,…

x*=Bx*+f

所以

e(k+1)=Be(k),

k=0,1,2,…遞推可得

e(k)=Bke(0),

k=0,1,2,…可見,當k時,

e(k)0Bk

O.對任意初始向量x(0),迭代法收斂

(B)<1.定理5.6.15.6迭代法的收斂性分析定理5.6.2若‖B‖=q<1,則對任意x(0),迭代法收斂,而且

由定理5.3.2知，

，再由定理5.6.1知，迭代法收斂。證明：先證收斂性。66證由于

x(k+1)=Bx(k)+f

,x(k)=Bx(k-1)+f,且

x*=Bx*+f；

所以

x(k+1)-x(k)=B(x

(k)-x(k-1)),x(k+1)–x*=B(x

(k)–x*)于是有

‖x(k+1)-x(k)‖q‖x

(k)-x(k-1)‖,‖x(k+1)–x*‖q‖x

(k)–x*‖

‖x(k+1)-x(k)‖=‖(x

(k+1)–x*)-(x(k)–x*)‖

‖x

(k)–x*‖-‖x(k+1)–x*‖(1-q)‖x(k)–x*‖所以定理5.6.2只是收斂的充分條件,并不必要,如則‖M‖1=1.2,‖M‖=1.3,‖M‖2=1.09,但(M)=0.8<1,所以迭代法是收斂的.可見,q越小收斂越快,且當‖x(