數(shù)值計(jì)算方法第三章講義課件

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

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

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

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

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

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

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

8、一、 LU分解法 /* LU Factorization. */就 n=3的情況分析順序消去法的消元過程.上一頁 下一頁 返回 矩陣的三角(LU)分解法三、直接法 可見, 消元過程相當(dāng)于下述矩陣乘法運(yùn)算：由分塊乘法可得：直接計(jì)算可得由(*)式可得上一頁 下一頁 返回 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分解的一般計(jì)算公式, 其結(jié)果與高斯消去法所得結(jié)果完全一樣. 避免了中間過程的計(jì)算,故稱為A的直接分解公式.上一頁 下一頁 返回 按上圖逐框求出矩陣A的LU分解,這種計(jì)算方法稱為緊湊格式法。上一頁 下一頁 返回 定理 若矩陣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)則求解原方程組等價(jià)于求解下面兩個(gè)方程組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、、平方根法 對(duì)稱 正定矩陣的分解法將對(duì)稱 正定陣 A 做 LU 分解記為定理 設(shè)n階對(duì)稱正定矩陣A，則存在唯一的單位下三角陣L及對(duì)角陣D 使得 。 稱為矩陣A 的喬累斯基分解上一頁 下一頁 返回 定理 設(shè)矩陣A對(duì)稱正定，則存在唯一的對(duì)角元為正的下三角陣 L，使得 。 稱為對(duì)稱正定矩陣A 的喬累斯基分解 利用喬累斯基（Cholesky）分解式來求解Ax=b的方法也稱Cholesky方法或平方根法三、三對(duì)角方程組追趕法若A滿足Gauss消去法可行的條件，則可用LU分解法求解其中：上一頁 下一頁 返回 解方程組Ax=d分為兩步，即求解Ly=d和Ux=y，計(jì)算公式如下：上述方法為求解三對(duì)角方程組的追趕

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

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

15、|+|y| (三角不等式).則稱f(x)=|x|是Rn(或Cn)上的一個(gè)向量范數(shù)(或模).定義3.2(向量的范數(shù)) 如果向量xRn(或Cn)的某個(gè)實(shí)值函數(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, 使得對(duì)一切xRn有 證明 只要就|x|s=|x|證明上式成立即可, 即證明存在常數(shù)c1,c20, 使得 考慮n元函數(shù)記S=x|x|=1, xRn, 則S是一個(gè)有界閉集. 由于f

16、(x)為S上的連續(xù)函數(shù), 所以f(x)于S上達(dá)到最大值和最小值, 即存在x, xS使得設(shè)xRn且x0, 則x/|x|S, 從而有即 可以證明常用范數(shù)滿足下述關(guān)系以及 注意 定理3.2不能推廣到無窮維空間. 由定理3.2可得結(jié)論: 如果在一種范數(shù)意義下向量序列收斂時(shí), 則在任何一種范數(shù)意義下該向量序列均收斂. 定理3.3 , 其中|為向量的任一種范數(shù). 證明 顯然, 而對(duì)于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的某個(gè)實(shí)值函數(shù)N(A)=|A|，滿足條件: (1) |A| 0 (|A|=0當(dāng)且僅當(dāng)A=0) (正定性)； (2) |cA|=|c| |A|, c為實(shí)數(shù) (齊次性)； (3) 對(duì)任意A,B有|A+B| |A|+|B| (三角不等式) (4) 對(duì)任意A,B有|AB| |A| |B| (相容性)；則稱 N(A) 是Rnn上的一個(gè)矩陣范數(shù)(或模). 上面我們定義的F(A)=|A|F就是Rnn上的一個(gè)矩陣范數(shù). 上述條件(即不等式)就稱為矩陣范數(shù)與向量范數(shù)的

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

19、5.7)是顯然的. 現(xiàn)只驗(yàn)證定義4中條件(4). 由(5.7), 有當(dāng)x0, 有故顯然這種矩陣的范數(shù)|A|v依賴于向量范數(shù)|x|v的具體含義. 也就是說, 當(dāng)給出一種具體的向量范數(shù)|x|v時(shí), 相應(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看出, 計(jì)算一個(gè)矩陣的|A|, |x|1還是比較容易的, 而矩陣的2-范數(shù)在|x|2計(jì)算上不方便, 但是矩陣的2-范數(shù)具有許多好的性質(zhì), 它在理論上是非常有用的. 我們指出, 對(duì)于復(fù)矩陣(即ACnn)定理3.5中1,2顯然也成立, 對(duì)于3應(yīng)改為例7 設(shè)則求相應(yīng)的各種范數(shù).解因?yàn)榱畹肁TA的特征值故可見定理3.7(特征值上界) 設(shè)ARnn, 則(A)|A|,即A的譜半徑不超過A的任何一種算子范數(shù)(對(duì)|A|F亦成立). 證明 設(shè)是A的任一特征值, x為相應(yīng)的特征向量,則Ax=x, 由(5.7)得注意到|x|0, 即得|A|. 定理3.6 , 其中|為矩陣的任一種范數(shù).定理3.8定理

21、3.9 如果|A|1時(shí),則(6.3)是“病態(tài)”的(即A是“病態(tài)”矩陣, 或者說A是壞條件的, 相對(duì)于方程組), 當(dāng)A的條件數(shù)相對(duì)的小, 則(6.3)是“良態(tài)”的(或者說A是好條件的). 注意, 方程組病態(tài)性質(zhì)是方程組本身的特性. A的條件數(shù)越大, 方程組的病態(tài)程度越嚴(yán)重, 也就越難用一般的計(jì)算方法求得比較準(zhǔn)確的解.注: cond (A) 的具體大小與 | | 的取法有關(guān)，但相對(duì)大小一致。 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為對(duì)稱矩陣時(shí)其中1, n為A的絕對(duì)值最大和絕對(duì)值最小的特征值. 條件數(shù)的性質(zhì): 1. 對(duì)任何非奇異矩陣, 都有cond(A)v1. 事實(shí)上 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.精確解為例計(jì)算cond 2(A) 。A1 = 解：考察 A 的特征根39206 1 測(cè)試病態(tài)程度：給一個(gè)擾動(dòng)，其相對(duì)誤差為此時(shí)精確解為2.0102 200%上一頁 下一頁 返回 例：Hilbert 陣cond (H2) =27cond (H3) 748cond (H6) =2.9 106cond (Hn) as n 注：一般判斷矩陣是否病態(tài)，并不計(jì)算A1，而由經(jīng)驗(yàn)得出。 行列式的值很大或很?。ㄈ缒承┬?、列近似相關(guān)）； 元素間的數(shù)量級(jí)相差大，且無規(guī)則； 主元消去過程中出現(xiàn)小主元； 特征值相差大數(shù)量級(jí)。上一頁 下一頁 返回 直接法: 經(jīng)過有限

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

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