第二章解線性方程組的直接方法2017秋

1高斯消去法主元素法直接三角分解法平方根法與改進平方根法誤差分析第二章線性方程組的直接方法2討論n元線性方程組的直接解法.3其中方程組的矩陣形式為若矩陣A非奇異,即det(A)≠0,則方程組有唯一解.4直接解法是指,若不考慮計算過程中的舍入誤差，經過有限次算術運算就能求出線性方程組的精確解的方法.但由于實際計算中舍入誤差的存在,用直接解法一般也只能求出方程組的近似解.Cramer法則是一種不實用的直接法，下面介紹幾種實用的直接法.解線性方程組的方法:直接方法和迭代法.迭代法是從解的某個近似值出發(fā),通過構造一個無窮序列去逼近精確解的方法.一般地,有限步內得不到精確解.5Gauss消元法是一種規(guī)則化的消元法，其基本思想是通過逐次消元計算，把一般線性方程組的求解轉化為等價的上三角形方程組的求解?！?Gauss消去法6消去后兩個方程中的x1得再消去最后一個方程的x2得消元結束.經過回代得解:例考慮線性方程組順序Gauss消去法7消元過程:先逐次消去變量x1,x2,將方程組化為同解的上三角形方程組.上過方法可推廣到一般情況回代過程:按方程相反的順序求解上三角形方程組.8第一步.設依次用乘矩陣的第1行加到第i行上,得到矩陣:則線性方程組的增廣矩陣為記9其中第二步.設依次用乘矩陣的第2行加到第i行,得到矩陣:10其中11這就完成了消元過程.如此繼續(xù)消元下去,第n1步結束后得到矩陣:12對此方程組進行回代，就可求出方程組的解.對應的方程組變成：13能用順序Gauss消去法求解的條件是在消元過程中得到的主元必須全不為0，即順序Gauss消去法通常也簡稱為Gauss消去法.主元素都不為零矩陣A的各階順序主子式都不為零.

順序Gauss消去法中的稱為主元素.14順序Gauss消去法求解n元線性方程組的乘除運算量

第1次消元乘除運算量:

消元過程乘除運算量求li1:

(n－1)求aij(2):(n－1)2求bi(2):(n－1)共(n2－1)次15

回代過程乘除運算量：求xn:

1求xn－1:2求x1:n….16n=30時,順序Gauss消去法只需9890次乘除法運算.順序Gauss消去法求解n元線性方程組的乘除運算量17高斯消去法優(yōu)缺點：簡單易行要求主元均不為零,因而適用范圍小數值穩(wěn)定性差18例：單精度解方程組

精確解8個8個用順序Gauss消去法計算：8個小主元可能導致計算失敗.§2主元素法19若將方程組改寫成:用順序Gauss消去法,消元得回代得解:x2=1,x1=1與準確解非常接近.可見,第一種算法是不穩(wěn)定的,第二種算法是穩(wěn)定的.20此例說明,在消元過程中,應避免選取絕對值較小的數作主元.如例中的第二種解法,通過交換方程次序,選取絕對值較大的元素作為主元.基于這種想法導出了主元法.為了提高計算的數值穩(wěn)定性,在消元過程中采用選擇主元的方法.常采用的是列主元消去法和全主元消去法.21給定線性方程組Ax=b,記A(1)=A,b(1)=b,列主元Gauss消去法的具體過程如下:

首先在增廣矩陣B(1)=(A(1),b(1))的第一列元素中,取

然后進行第一步消元得增廣矩陣B(2)=(A(2),b(2)).列主元消去法22

然后進行第二步消元得增廣矩陣B(3)=(A(3),b(3)).按此方法繼續(xù)進行下去,經過n1步選主元和消元運算,得到增廣矩陣B(n)=(A(n),b(n)).則方程組A(n)x=b(n)是與原方程組等價的上三角形方程組,可進行回代求解.只要|A|0,列主元Gauss消去法就可順利進行.

再在矩陣B(2)=(A(2),b(2))的第二列元素中,取23

全主元素法每一步選絕對值最大的元素為主元素，保證Stepk:①選?、贗fik

k

then交換第k行與第ik

行;Ifjk

k

then交換第k列與第jk

列;③消元注：列交換改變了xi

的順序，須記錄交換次序，解完后再換回來。24例用主元素法求解線性方程組計算過程保留三位小數,方程的精確解為x1*=1,x2*=2,x3*=3.25解1.按列主元素法，求解過程如下消元回代得26解2.按全主元素法，求解過程如下回代得27全主元素法的精度優(yōu)于列主元素法,這是由于全主元素是在全體系數中選主元,故它對控制舍入誤差十分有效.但全主元素法在計算過程中,需同時作行與列的互換,因而程序比較復雜,計算時間較長.

列主元素法的精度雖然稍低于全主元素法,但其計算簡單,工作量大為減少,且計算經驗與理論實踐均表明,它與全主元素法同樣具有良好的數值穩(wěn)定性.列主元素法是求解中小型稠密線性方程組的最好方法之一.

例3的計算結果表明28Gauss消元法的矩陣表示§3直接三角分解法兩者等價29其中兩者等價n=3時Gauss消元法的矩陣表示30其中兩者等價3132條件:矩陣A經Gauss消元法后得到的上三角矩陣.33例求矩陣的三角分解.解34上述n=3的情形可以推廣到一般情形,例如n=43536條件:矩陣A經Gauss消元法后得到的上三角矩陣.37例38條件:矩陣A經Gauss消元法后得到的上三角矩陣.Gauss消元法的矩陣表示Doolittle分解39矩陣的三角分解定理

設A為n階方陣,若A

的前n階順序主子式Ai

(i=1,2,…,n)均不為0,則矩陣A存在唯一的Doolittle分解.存在性：唯一性：反證法40

Doolittle分解中LU

元素的求解次序Doolittle分解中LU

的另一求法41U的第一行L的第一列

Doolittle分解中LU

元素的求解右端用矩陣乘法展開,比較兩邊的第1行和第1列得42假設已求得U的前(r

1)行和L的前(r

1)列,r>1.下面求U的第r行和L的第r列.右端用矩陣乘法,比較兩邊的第r行的后(n－r+1)個元素arj

=L的第r行行向量與U的第j列列向量的內積(j

r)U的第r行43U的第r行44右端用矩陣乘法,比較兩邊的第r列的后(n－r)個元素air=L的第i行行向量與U的第r列列向量的內積(i

>r)L的第r列45L的第r列46對r=2,3,…,n,矩陣三角分解的緊湊格式47例用緊湊格式求矩陣的三角分解.解48先求Ly=b,得y;再求Ux=y,得x.直接三角分解法或Doolittle分解法.Doolittle分解法49Ly=b乘除運算量為50Ux=y乘除運算量為51例用直接三角分解法解解52先求Ly=b得y再求Ux=y

得x53由于在求出uij,lij和yi后,aij和bi就無需保留了,故上機計算時,可把L,U和y存在A,b所占單元,回代時x取代y,整個計算過程中不需要增加新的存儲單元.在求一系列系數矩陣相同而右端項不同的線性方程組Ax=b(k),(k=1,2,…,m)時(如求逆矩陣),用三角分解法更為簡便.每解一個方程組Ax=b(k)

僅需要增加n2

次乘除法運算.解線性方程組Ax=b的Doolittle三角分解法的計算量約為n3/3,與Gauss消去法相同.54Crout分解定理

設A為n階方陣,若A

的前n階順序主子式Ai

(i=1,2,…,n)均不為0,則矩陣A可以唯一分解為其中L為下三角陣,U為單位上三角陣.55特殊的稀疏矩陣解三對角方程組的追趕法解三對角方程組56追趕法是求三對角線性方程組的三角分解法.

追趕法本質上是Gauss消元法.二階常微分方程邊值問題的差分離散方程組,熱傳導方程以及船體數學中建立的三次樣條函數的三轉角或三彎矩方程組均為三對角方程組.57例

4階三對角矩陣的三角分解58例

4階三對角矩陣的三角分解單位下二對角陣上二對角陣59定理設三對角矩陣A滿足下列條件則它可以分解成A=LU三對角矩陣的三角分解對角占優(yōu)陣單位下二對角陣上二對角陣工程中得到的三對角陣多數滿足此條件60三對角矩陣三角分解中LU的求解次序61對k=3,…,n

三對角矩陣三角分解中LU的求解右端用矩陣乘法展開,比較兩邊元素得乘除運算量2n－262解三對角方程組的追趕法先求Ly=d,得y;再求Ux=y,得x.追：消元過程趕：回代過程63Ly=d乘除運算量為n－164Ux=y乘除運算量為2n－1追趕法總的乘除運算量5n－465追趕法的實質就是Gauss消元法,只是由于系數中出現了大量的零,在計算過程中將它們撇開,從而使計算公式大大簡化,也大大減少了計算量.為節(jié)省計算機存儲單元,計算得到的lk,uk分別存放在ak,bk的存儲單元內,而yk,xk

存放在dk

的存儲單元內.當系數矩陣為滿足定理條件的嚴格對角占優(yōu)陣時,追趕法具有良好的數值穩(wěn)定性.66§4平方根法與改進的平方根法求解對稱正定方程組,即其中A為對稱正定陣.利用對稱性,平方根法的計算量是Gauss消去法的一半.67對稱正定陣的Cholesky分解定理

若A

是對稱正定陣,則存在唯一的非奇異下三角陣L,使得且L的對角元素皆為正數,即lii>0(i=1,2,…,n).證明

A可進行Doolittle分解.A為對稱正定陣,它的n個順序主子式均大于0,設其中為單位下三角陣,U為上三角矩陣.用反證法68對U進一步分解.

A對稱

A正定對角矩陣D的對角元均為正.令則其中為非奇異下三角陣,且對角元均為正.

A=LLT.69對稱

A為3階對稱正定陣,A=LLT,怎樣求L?對稱70例對下列矩陣進行Cholesky分解71例對下列矩陣進行Cholesky分解Matlab解法：A=[4-26;-2175;6522]R=chol(A)L=R’72

A為n階對稱正定陣,A=LLT,怎樣求L?L中元素的求解次序依次求L的第一列,第二列,…,第n列.73………...對稱

A為n階對稱正定陣,A=LLT,怎樣求L?li為L的第i

個行向量74對稱

A為n階對稱正定陣,A=LLT,怎樣求L?75求L的第一列求L的第二列

A為n階對稱正定陣,A=LLT,怎樣求L?76

A為n階對稱正定陣,A=LLT,怎樣求L?設已經求得L的前k－1列,現求L的第k列(k=3,4,…,n)77例求正定陣的Cholesky分解.解78例求正定陣的Cholesky分解.解79先求Ly=b,得y,再求LTx=y,得x.解正定線性方程組的平方根法或Cholesky分解法.平方根法或Cholesky分解法設A為對稱正定陣80定理

若A

是對稱正定陣,則存在唯一的單位下三角陣L和對角陣D,使得且D的對角元素皆為正數.對稱正定陣的LDLT分解證明

A對稱

A正定對角矩陣D的對角元均為正.81對稱正定陣的LDLT分解本質上是對A作Doolittle分解,即LU分解.LDLT分解中的D=LU分解中的U的對角部分LDLT分解中的L=LU分解中的L82對稱正定矩陣A的LU分解,計算量可以節(jié)省一半求U的第1行求L的第1列對稱正定陣的LDLT分解中L,D的計算先對對稱正定陣A作LU分解83求U的第k行(k=2,3,…,n)求L的第k列(k=2,3,…,n)對稱正定陣的LDLT分解中L,D的計算節(jié)省了計算量求D84例求矩陣的LDLT分解.解8586解正定線性方程組的改進平方根法或LDLT分解法.先求Ly=b,得y,再求LTx=D1y,得x.改進平方根法或LDLT分解法設A為對稱正定陣87平方根法與改進的平方根法的優(yōu)點計算無須選主元,由于正定性,計算過程是數值穩(wěn)定的計算量是Gauss消元法的一半由于對稱性,實際計算可存儲一半是求解中小型稠密正定線性方程組的好算法88§5誤差分析用直接法解線性方程組，初始數據會有誤差，計算過程同樣會產生誤差，這就需要對這些誤差作一些分析.主要數學工具：向量范數，矩陣范數，條件數89向量范數向量范數是用來度量向量長度的,它可以看成是二、三維解析幾何中向量長度概念的推廣.定義(向量范數)對任一向量xRn，按照一定規(guī)則確定一個實數與它對應，該實數記為||x||，若||x||滿足下面三個性質：1)||x||0；||x||=0當且僅當x=0(零向量)

正定性2)對任意實數，||x||=||||x||齊次性3)對任意向量x,yRn，||x+y||||x||+||y||三角不等式則稱該實數||x||為向量x的范數.90Rn中常用的三種范數其中x1,x2,…,xn分別是x的n個分量.1-范數2-范數-范數用向量范數的定義來驗證,即驗證它們滿足向量范數定義中的三個性質.向量的模91定義(范數等價)

在Rn中有兩個范數||||和||||,若存在實數M,m>0,使得對任意的n維向量x,都有則稱這兩個范數等價.Rn中范數的重要性質:范數等價定理92范數等價定理:

Rn中任意兩個范數等價.Rn中范數的重要性質:范數等價定理例1-范數,2-范數和-范數是兩兩等價的.93當不需要指明使用哪一種向量范數時，就用記號||.||泛指任何一種向量范數.有了向量的范數就可以用它來衡量向量的大小和表示向量的誤差.設x

為Ax=b

的精確解，x*為其近似解絕對誤差相對誤差94矩陣范數是用于定義矩陣“大小”的量，類似于向量范數，可以定義n階方陣A的范數.定義(矩陣范數)

設A為n階方陣，按照一定規(guī)則有一實數與之對應，記為||A||,若||A||滿足：1)||A||0,||A||=0當且僅當A=0時；2)對任意實數，||A||=||||A||；3)對任意兩個n階方陣A,B,都有||A+B||||A||+||B||;4)||AB||||A||||B||(相容性條件)則稱||A||為矩陣A的范數.矩陣范數95常用的三種矩陣范數1-范數或列范數2-范數或譜范數-范數或行范數用矩陣范數的定義來驗證,即驗證它們滿足矩陣范數定義中的四個性質.96定理設A為n階方陣，||||是Rn中的向量范數,則是一種矩陣范數,稱其為由向量范數||||誘導出的矩陣范數.由向量范數誘導出的矩陣范數矩陣范數與向量范數的相容性性質：設矩陣范數是由向量范數誘導出的，則對任意n階方陣A,以及任意的n維向量x,有979899常用的三種矩陣范數均是由向量范數誘導出的.

對于給定的向量范數1-范數,2-范數及-范數,可以證明由它們誘導出的矩陣范數分別為1-范數或列范數2-范數或譜范數-范數或行范數100101實際計算常用1-范數與-范數因其計算比較簡單,理論證明常用2-范數因為它有一些好性質.由矩陣范數與向量范數的相容性性質知對任意n階方陣A,以及任意的n維向量x,有誤差分析中常用102矩陣的誤差可用矩陣范數表示.

設A*是A的近似矩陣絕對誤差相對誤差矩陣范數的等價定理也成立.103Matlab函數：normNORM(X)isthelargestsingularvalueofX,max(svd(X)).NORM(X,2)isthesameasNORM(X).NORM(X,1)isthe1-normofX,thelargestcolumnsum,=max(sum(abs((X)))).NORM(X,inf)istheinfinitynormofX,thelargestrowsum,=max(sum(abs((X')))).NORM(X,'fro')istheFrobeniusnorm,sqrt(sum(diag(X'*X))).NORM(X,P)isavailableformatrixXonlyifPis1,2,infor'fro'.

104NormofAMatrixMatlabcommand:norm(A,1)norm(A,2)=norm(A)norm(A,inf)norm(A,’fro’)105方程組的狀態(tài)與條件數例方程組I方程組II右端項有0.00001的差別,最大相對誤差為0.5105,但解分量的相對誤差為50%.平面上兩條接近于平行的直線的交點，當其中一條直線稍有變化時,新的交點可與原交點相差甚遠.106當一個方程組,由于系數矩陣或右端項有微小擾動,而引起解發(fā)生巨大變化時,則稱該方程組是“病態(tài)”的.分以下兩種情形加以討論只有右端項有擾動只有系數矩陣有擾動怎樣刻畫線性方程組“病態(tài)”程度.用條件數來描繪.107只有右端項有擾動原方程組擾動后方程組右端項有擾動b,引起解的變化為x,其相對誤差為它究竟有多大?108故又這表明:當右端項有擾動b時,解的相對誤差不超過右端項的相對誤差的||A||||A1||倍.只有右端項有擾動從而109只有系數矩陣有擾動原方程組擾動后方程組系數矩陣有擾動A,引起解的變化為x,其相對誤差為它究竟有多大?110故這表明:當系數矩陣有擾動A,解的擾動不超過系數矩陣相對誤差的||A|