1、第3章 線性方程組的數(shù)值解法三、 矩陣的三角分解及其應(yīng)用四、 線性代數(shù)方程組的性態(tài)與誤差分析五、 迭代法一、 線性方程組二、 高斯(Gauss)消元法六、 共軛梯度法求解上一頁 下一頁 返回 一、線性方程組 實際問題中的線性方程組分類：按系數(shù)矩陣中零元素的個數(shù)：稠密線性方程組稀疏線性方程組按未知量的個數(shù)：高階線性方程組低階線性方程組(如1000)(80%)按系數(shù)矩陣的形狀對稱正定方程組三角形方程組三對角占優(yōu)方程組 實際中，存在大量的解線性方程組的問題。很多數(shù)值方法到最后也會涉及到線性方程組的求解問題：如樣條插值的M和m關(guān)系式，曲線擬合的法方程，方程組的Newton迭代等問題。解線性方程組的兩類

2、數(shù)值方法：直接法: 經(jīng)過有限次四則運算后可求得方程組精確解的方法(不計舍入誤差!)迭代法：從解的某個近似值出發(fā)，通過構(gòu)造一個無窮序列去逼近精確解的方法。（一般有限步內(nèi)得不到精確解）特點：準確，可靠，理論上得到的解是精確的。適合于中小規(guī)模的方程組。特點：速度快，但有誤差。適合于大規(guī)模的方程組。Cramer法則是一種直接法，但無法實用 直接法概述直接法是將原方程組化為一個或若干個三角形方程組的方法對于線性方程組根據(jù)Cramer(克萊姆)法則,若則都是三角形方程組上述方法稱為直接三角形分解法不論是Gauss消元法還是直接三角形分解法,最終都歸結(jié)為解三角形方程組一、三角形線性代數(shù)方程組的直接法若記下三

3、角形線性方程組上三角形線性方程組 Gauss消去法二、直接法 其解為其解為:按回代過程求解對方程組，作如下的變換，解不變：交換兩個方程的次序一個方程的兩邊同時乘以一個非0的數(shù)一個方程的兩邊同時乘以一個非0數(shù)，加到另一個方程二、Gauss消元法因此，對增廣矩陣(A，b )作如下的變換，解不變：交換矩陣的兩行某一行乘以一個非0的數(shù)某一行乘以一個非0數(shù)，加到另一行消元法就是對增廣矩陣作上述行的變換，變?yōu)槿切尉€性方程組，而后求根。用矩陣表示：=首先將A化為上三角陣，再回代求解。例3.1.1解化為同解方程組化為上三角方程組按回代過程求解上述求解三元線性代數(shù)方程組的方法即為Gauss消元法。(1) 定義

4、行乘數(shù)相當(dāng)于第i個方程-第一個方程行乘數(shù)新的第i方程同解！第一個方程不動！ 且(1) 定義行乘數(shù)相當(dāng)于第i個方程-第二個方程行乘數(shù)新的第i方程同解！第一、二個方程不動！ 系數(shù)矩陣與常數(shù)項： Gauss消元法的運算量乘法次數(shù)：除法次數(shù)：全部回代過程需作乘除法的總次數(shù)為于是Gauss消元法的乘除法運算總的次數(shù)為數(shù)級類似地，可得Gauss消元法的加減法運算總的次數(shù)為由于計算機完成一次乘除運算所耗時間要遠遠多于完成一次加減運算的時間，而且按統(tǒng)計規(guī)律，在一個算法中，加減運算和乘除運算次數(shù)大體相當(dāng)，故在衡量一個算法的運算量時只需統(tǒng)計乘除的運算次數(shù)。Gauss消元法乘除法約為2700次而如果用Cramer法

5、則的乘除法運算次數(shù)約為用行列式定義回代求解 從Gauss消元法的過程自然想到：如果在消元過程中把對角線上方未知量的系數(shù)也化為零，使方程組化成每個方程只含一個未知量的方程組，則無須回代就可得到解。這種消元法稱為GaussJordan消元法。但經(jīng)計算，它需要的乘除法次數(shù)為 顯然，其運算量超過Gauss消元法。所以用GaussJordan消元法求解線性代數(shù)方程組是不可取的。不過，用它來求逆矩陣卻有方便之處。三、選主元Gauss消元法例3.1.3用Gauss消元法解線性方程組(用3位十進制浮點數(shù)計算）解本方程組的精度較高的解為用Gauss消元法求解(用3位十進制浮點數(shù)計算)回代后得到主元9999與精確

6、解相比,該結(jié)果相當(dāng)糟糕！究其原因,在求行乘數(shù)時用了很小的數(shù)0.0001作除數(shù)! 上述例題表明：Gauss消元法是不穩(wěn)定的算法，其中小主元是不穩(wěn)定的根源。因而在Gauss消元法中要避免小主元的出現(xiàn)。這就需要采用“選主元”的技術(shù)。所謂選主元，簡單地說就是選取絕對值最大的元素作主元。 全選主元Gauss消元法在此范圍內(nèi)選主元需進行元素比較 列選主元Gauss消元法在此范圍內(nèi)選主元不必選主元的情況： 當(dāng)系數(shù)矩陣A對稱正定或嚴格對角占優(yōu)或不可約對角占優(yōu)時，可不必選主元?，F(xiàn)用列選主元消元法重解例3.1.3：將方程組的1,2行交換，即得回代后得到這是一個相當(dāng)不錯的結(jié)果！0.9999例3.1.4解線性方程組(

7、用8位十進制尾數(shù)的浮點數(shù)計算)解這個方程組和例3.1.3一樣，若用Gauss消元法計算會有小數(shù)作除數(shù)的現(xiàn)象，若采用換行的技巧（即列選主元消元法） ，則可避免。絕對值最大不需換行經(jīng)過回代后可得事實上,方程組的準確解為 列選主元Gauss消元法的計算步驟? 列選主元Gauss消元法的流程圖開始輸出無解信息消元換行停機回代求解本算法主要由四個模塊組成 一些特殊情況, 主元就在對角線上,不需選主元. 元素滿足如下條件的矩陣 即對角線上每一元素的絕對值均大于同行其他各元素絕對值之和,這樣的矩陣稱為按行嚴格對角占優(yōu)矩陣,簡稱嚴格對角占優(yōu)矩陣.例：性質(zhì)：嚴格對角占優(yōu)矩陣必定非奇異.上一頁 下一頁 返回 補充

8、一、 LU分解法 /* LU Factorization. */就 n=3的情況分析順序消去法的消元過程.上一頁 下一頁 返回 矩陣的三角(LU)分解法三、直接法 可見, 消元過程相當(dāng)于下述矩陣乘法運算：由分塊乘法可得：直接計算可得由(*)式可得上一頁 下一頁 返回 Step 1:記M(1) =，則Step n 1: 其中 M (k)= 以上分析推廣到n階線性方程組可得上一頁 下一頁 返回 記為L單位下三角陣/* unitary lower-triangular matrix */記 U =A 的 LU 分解/* LU factorization */也稱 A 的Doolittle分解上一頁

9、下一頁 返回 道立特分解法 /* Doolittle Factorization */： LU 分解的緊湊格式 /* compact form */根據(jù)矩陣乘法法則,先比較等式兩邊第1行和第1列元素有：上一頁 下一頁 返回 設(shè)已定出U 的第1行到第k-1行的元素 L 的第1列到第k-1列的元素上一頁 下一頁 返回 (1),(2)兩式就是 LU分解的一般計算公式, 其結(jié)果與高斯消去法所得結(jié)果完全一樣. 避免了中間過程的計算,故稱為A的直接分解公式.上一頁 下一頁 返回 按上圖逐框求出矩陣A的LU分解,這種計算方法稱為緊湊格式法。上一頁 下一頁 返回 定理 若矩陣A非奇異, 則A能分解為LU 的充

10、分必要條件是A的所有順序主子式 均不為0.定理 若非奇異矩陣A有LU 分解,則此分解是唯一的.上一頁 下一頁 返回 矩陣的Doolittle分解法的Matlab程序function 1,u=lu_Doolittle1(A)%A為可逆方陣，l為返回下三角矩陣，u為上三角陣n=length(A);u=zeros(n);l=eye(n);u(1,:)=A(1,:);l(2:n,1)=A(2:n,1)/u(1,1);for k=2:n u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n); l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1

11、:k-1,k)/u(k,k);end例：利用系數(shù)矩陣的LU分解, 求解方程組解：LU分解的緊湊格式為上一頁 下一頁 返回 在Matlab命令窗口執(zhí)行A=1 1 1 1;1 2 1 3;4 3 2 1;6 11 3 7;l,u=lu_Doolittle1(A)則求解原方程組等價于求解下面兩個方程組Ly=b及Ux=y上一頁 下一頁 返回 在Matlab命令窗口執(zhí)行A=1 1 1 1;1 2 1 3;4 3 2 1;6 11 3 7;b=6 12 14 45;y,x=LU_s(A,b)A=LU若U為單位上三角矩陣，則稱為Crout分解.Crout分解相應(yīng)的求解公式可類似得到。上一頁 下一頁 返回 二

12、、平方根法 對稱 正定矩陣的分解法將對稱 正定陣 A 做 LU 分解記為定理 設(shè)n階對稱正定矩陣A，則存在唯一的單位下三角陣L及對角陣D 使得 。 稱為矩陣A 的喬累斯基分解上一頁 下一頁 返回 定理 設(shè)矩陣A對稱正定，則存在唯一的對角元為正的下三角陣 L，使得 。 稱為對稱正定矩陣A 的喬累斯基分解 利用喬累斯基（Cholesky）分解式來求解Ax=b的方法也稱Cholesky方法或平方根法三、三對角方程組追趕法若A滿足Gauss消去法可行的條件，則可用LU分解法求解其中：上一頁 下一頁 返回 解方程組Ax=d分為兩步，即求解Ly=d和Ux=y，計算公式如下：上述方法為求解三對角方程組的追趕

13、法，也稱Thomas算法.上一頁 下一頁 返回 追趕法公式簡單，計算量和存儲量都小，整個求解過程只需要5n-4次乘除運算。 在分析方程組的解的誤差及下章中迭代法的收斂性時，常產(chǎn)生一個問題,即如何判斷向量 x 的“大小”，對矩陣也有類似的問題. 本節(jié)介紹n維向量范數(shù)和nn矩陣的范數(shù). 向量范數(shù)是三維歐氏空間中向量長度概念的推廣，在數(shù)值分析中起著重要作用.四、線性代數(shù)方程組的性態(tài)與誤差分析 首先將向量長度概念推廣到Rn(或Cn)中.稱為向量x, y的數(shù)量積(內(nèi)積). 將非負實數(shù)稱為向量x的歐氏范數(shù). 定義3.1 設(shè)x=(x1,x2,xn)T, y= (y1,y2,yn)TRn , 將實數(shù) 下面定理

14、可在線性代數(shù)書中找到. 定理 設(shè)x, yRn (或Cn), 則 1. (x, x)=0, 當(dāng)且僅當(dāng)x=0時成立; 我們還可以用其他辦法來度量Rn中向量的“大小”. 例如, 對于x=(x1,x2)TRn, 可以用一個x的函數(shù)N(x)=max|xi|來度量 x 的“大小”, 而且這種度量 x 的“大小”的方法計算起來比歐氏范數(shù)方便. 在許多應(yīng)用中, 對度量x的“大小”的函數(shù)N(x)都要求是正定的、齊次的且滿足三角不等式. 下面我們給出向量范數(shù)的一般定義. (1) |x|0(|x|=0當(dāng)且僅當(dāng)x=0) (正定性), (2) |x|=| |x|, 對任何R(或C)(齊次性), (3) |x+y| |x

15、|+|y| (三角不等式).則稱f(x)=|x|是Rn(或Cn)上的一個向量范數(shù)(或模).定義3.2(向量的范數(shù)) 如果向量xRn(或Cn)的某個實值函數(shù)f(x)=|x|，滿足條件: 由(3)可推出不等式. (4) | |x|-|y| | |x-y|. 下面給出幾種常用的向量范數(shù), 設(shè)x=(x1,x2,xn)T.容易證明前三種范數(shù)是的p-范數(shù)特殊情況, 其中向量的-范數(shù)(最大范數(shù))向量的1-范數(shù)向量的p-范數(shù)(0p0, 使得對一切xRn有 證明 只要就|x|s=|x|證明上式成立即可, 即證明存在常數(shù)c1,c20, 使得 考慮n元函數(shù)記S=x|x|=1, xRn, 則S是一個有界閉集. 由于f

16、(x)為S上的連續(xù)函數(shù), 所以f(x)于S上達到最大值和最小值, 即存在x, xS使得設(shè)xRn且x0, 則x/|x|S, 從而有即 可以證明常用范數(shù)滿足下述關(guān)系以及 注意 定理3.2不能推廣到無窮維空間. 由定理3.2可得結(jié)論: 如果在一種范數(shù)意義下向量序列收斂時, 則在任何一種范數(shù)意義下該向量序列均收斂. 定理3.3 , 其中|為向量的任一種范數(shù). 證明 顯然, 而對于Rn上任一種范數(shù)|, 由定理15, 存在常數(shù)c1,c20, 使于是又有 下面我們將向量范數(shù)概念推廣到矩陣上去. 可以將Rnn中的矩陣A=(aij)nn當(dāng)作n2維向量, 則由向量的2-范數(shù)可以得到Rnn中矩陣的一種范數(shù)稱為A的F

17、robenius范數(shù). |A|F顯然滿足正定性、齊次性及三角不等式. 下面我們給出矩陣范數(shù)的一般定義. 定義3.4(矩陣的范數(shù)) 如果矩陣ARnn的某個實值函數(shù)N(A)=|A|，滿足條件: (1) |A| 0 (|A|=0當(dāng)且僅當(dāng)A=0) (正定性)； (2) |cA|=|c| |A|, c為實數(shù) (齊次性)； (3) 對任意A,B有|A+B| |A|+|B| (三角不等式) (4) 對任意A,B有|AB| |A| |B| (相容性)；則稱 N(A) 是Rnn上的一個矩陣范數(shù)(或模). 上面我們定義的F(A)=|A|F就是Rnn上的一個矩陣范數(shù). 上述條件(即不等式)就稱為矩陣范數(shù)與向量范數(shù)的

18、相容性條件. 由于在大多數(shù)與估計有關(guān)的問題中, 矩陣和向量會同時參與討論, 所以希望引進一種矩陣的范數(shù), 它是和向量范數(shù)相聯(lián)系而且和向量范數(shù)相容的, 即對任何向量xRn及矩陣ARnn都成立|Ax|A| |x|. 為此我們再引進一種矩陣的范數(shù). 定義3.5(矩陣的算子范數(shù)) 設(shè)xRn, ARnn, 給出一種向量范數(shù)|x|v(如v=1,2或), 相應(yīng)地定義一個矩陣的非負函數(shù)可驗證|A|v滿足定義4(見下面定理), 所以|A|v是Rnn上矩陣的一個范數(shù), 稱為A的算子范數(shù). 定理3.4 設(shè)|x|v是Rn上的一個向量范數(shù), 則|A|v是Rnn上矩陣的范數(shù), 且滿足相容條件 證明 由(5.6)相容條件(

19、5.7)是顯然的. 現(xiàn)只驗證定義4中條件(4). 由(5.7), 有當(dāng)x0, 有故顯然這種矩陣的范數(shù)|A|v依賴于向量范數(shù)|x|v的具體含義. 也就是說, 當(dāng)給出一種具體的向量范數(shù)|x|v時, 相應(yīng)地就得到了一種矩陣范數(shù)|A|v. 定理3.5 設(shè)xRn, A=(aij)Rnn, 則有常用的算子范數(shù)（稱為A的行范數(shù)）,（稱為A的列范數(shù)）,（稱為A的2-范數(shù)）.其中max(ATA)表示ATA的最大特征值, 由于矩陣2-范數(shù)與ATA 的特征值有關(guān)，所以又稱為A的譜范數(shù). 定義3.6 設(shè)ARnn的特征值為i (i=1,2,n) , 稱為A的譜半徑.例子 設(shè), 求A的譜半徑.解得A的特征值故得A的譜半徑

20、為 由定理3.5看出, 計算一個矩陣的|A|, |x|1還是比較容易的, 而矩陣的2-范數(shù)在|x|2計算上不方便, 但是矩陣的2-范數(shù)具有許多好的性質(zhì), 它在理論上是非常有用的. 我們指出, 對于復(fù)矩陣(即ACnn)定理3.5中1,2顯然也成立, 對于3應(yīng)改為例7 設(shè)則求相應(yīng)的各種范數(shù).解因為令得ATA的特征值故可見定理3.7(特征值上界) 設(shè)ARnn, 則(A)|A|,即A的譜半徑不超過A的任何一種算子范數(shù)(對|A|F亦成立). 證明 設(shè)是A的任一特征值, x為相應(yīng)的特征向量,則Ax=x, 由(5.7)得注意到|x|0, 即得|A|. 定理3.6 , 其中|為矩陣的任一種范數(shù).定理3.8定理

21、3.9 如果|A|1時,則(6.3)是“病態(tài)”的(即A是“病態(tài)”矩陣, 或者說A是壞條件的, 相對于方程組), 當(dāng)A的條件數(shù)相對的小, 則(6.3)是“良態(tài)”的(或者說A是好條件的). 注意, 方程組病態(tài)性質(zhì)是方程組本身的特性. A的條件數(shù)越大, 方程組的病態(tài)程度越嚴重, 也就越難用一般的計算方法求得比較準確的解.注: cond (A) 的具體大小與 | | 的取法有關(guān)，但相對大小一致。 cond (A) 取決于A，與解題方法無關(guān)。常用條件數(shù)有：cond 1(A)cond (A)cond2 (A)條件數(shù)的性質(zhì)： A可逆，則 cond p (A) 1； A可逆， R 則 cond ( A) =

22、cond (A) ； A正交，則 cond 2 (A) =1； A可逆，R正交，則 cond 2 (RA) = cond 2 (AR) = cond (A)2 。上一頁 下一頁 返回 通常使用的條件數(shù), 有 (1) cond(A) =|A-1|A| ; (2) A的譜條件數(shù);當(dāng)A為對稱矩陣時其中1, n為A的絕對值最大和絕對值最小的特征值. 條件數(shù)的性質(zhì): 1. 對任何非奇異矩陣, 都有cond(A)v1. 事實上 cond(cA)v=cond(A)v. 2. 設(shè)A為非奇異矩陣且c0(常數(shù)), 則 3. 如果A為正交矩陣, 則cond(A)2=1; 如果A為非奇異矩陣, R為正交矩陣, 則 c

23、ond(RA)2=cond(AR)2=cond(A)2.精確解為例計算cond 2(A) 。A1 = 解：考察 A 的特征根39206 1 測試病態(tài)程度：給一個擾動，其相對誤差為此時精確解為2.0102 200%上一頁 下一頁 返回 例：Hilbert 陣cond (H2) =27cond (H3) 748cond (H6) =2.9 106cond (Hn) as n 注：一般判斷矩陣是否病態(tài)，并不計算A1，而由經(jīng)驗得出。 行列式的值很大或很小（如某些行、列近似相關(guān)）； 元素間的數(shù)量級相差大，且無規(guī)則； 主元消去過程中出現(xiàn)小主元； 特征值相差大數(shù)量級。上一頁 下一頁 返回 直接法: 經(jīng)過有限

24、次運算后可求得方程組精確解的方法(不計舍入誤差!) 直接法比較適用于中小型方程組。直接法得到的解理論上是準確的，但它們的計算量都是n3數(shù)量級，存儲量為n2數(shù)量級，這在n比較小的時候還比較合適（n400）。對于現(xiàn)在的很多實際問題，往往要我們求解n很大的高階方程組，而且這些矩陣往往是稀疏的。 對于這類矩陣，用直接法時在運算中很難保持稀疏性，程序設(shè)計復(fù)雜，且會耗費大量的時間和存儲單元。因此我們有必要引入一類新的方法：迭代法。四、迭代法 迭代法能保持矩陣的稀疏性，具有計算簡單，存儲量小，編制程序容易等優(yōu)點，并在許多情況下收斂較快。故迭代法能有效地解一些高階方程組。 迭代法的基本思想是構(gòu)造一串收斂到解的

25、序列，即建立一種從已有近似解計算新的近似解的規(guī)則。由不同的計算規(guī)則得到不同的迭代法，本節(jié)介紹的幾種迭代法都各有其限定的應(yīng)用范圍，這是因為對于一個給定的方程組，某些迭代法收斂得很快，而有些迭代法可能不收斂，或收斂得很慢，以至于無實用價值。 迭代法: 從解的某個近似值出發(fā),通過構(gòu)造一個無窮序列去逼近精確解的方法。(一般有限步內(nèi)得不到精確解) 迭代法的一般形式及其收斂性四、迭代法幾種古典迭代算法（一） 雅可比(Jacobi)迭代法五、迭代法建立迭代格式：記 Jacobi迭代的矩陣形式易知，Jacobi 迭代有 Jacobi迭代法的計算過程Jacobi迭代法分量形式的Matlab程序function

26、x,k=jacobif(A,b,x0,ep,Nmax)n=length(A);k=0;if nargin5 Nmax=500;endif nargin5 ep=1e-5;endif narginep&kNmax,k=k+1;x0=x; for i=1:n y(i)=b(i); for j=I y(i)=y(i)-A(i,j)*x0(j); end; end if abs(A(i,i)1e-10|k=Nmax warning(A(i,i)太小); return end y(i)=y(i)/A(i,i); end x=y;end（二） 逐次代換法（GaussSeidel迭代） GaussSeidel迭代的矩陣形式 Gauss-Seidel