數(shù)值分析 第三章 基于MATLAB的科學(xué)計(jì)算—線性方程組1

科學(xué)計(jì)算—理論、方法 及其基于ＭＡＴＬＡＢ的程序?qū)崿F(xiàn)與分析 三、 解線性方程組（線性矩陣方程） 解線性方程組是科學(xué)計(jì)算中最常見(jiàn)的問(wèn)題。所說(shuō)的“最常見(jiàn)”有兩方面的含義： １） 問(wèn)題的本身是求解線性方程組； ２） 許多問(wèn)題的求解需要或歸結(jié)為線性方程組的求解。 　　關(guān)于線性方程組 　　　　　?。ǎ保? 其求解方法有兩類： １） 直接法：高斯消去法（Gaussian Elimination）； ２） 間接法：各種迭代法（Iteration）。 １、高斯消去法 １） 引例 考慮如下（梯形）線性方程組： 高斯消去法的求解思路：把一般的線性方程組（１）化成（上或下）梯形的形式。 ２）高斯消去法——示例 考慮如下線性方程組： １） 第一個(gè)方程的兩端乘加到第二個(gè)方程的兩端，第一個(gè)方程的兩端乘 －１加到第三個(gè)方程的兩端，得 ２）　第二個(gè)方程的兩端乘加到第三個(gè)方程的兩端，得 ３） 從上述方程組的第三個(gè)方程依此求解，得 ３）高斯消去法的不足及其改進(jìn)——高斯（全、列）主元素消去法 在上例中，由于建模、計(jì)算等原因，系數(shù)2.001而產(chǎn)生0.0005的誤差，實(shí)際求解的方程組為 注：數(shù)值穩(wěn)定的算法 高斯列主元素消去法就是在消元的每一步選?。校┲髟亍涣兄薪^對(duì)值最大的元取做主元素，高斯列主元素消去法是數(shù)值穩(wěn)定的方法。 列主元素消去法的基本思想：在每輪消元之前，選列主元素（絕對(duì)值最大的元素）,使乘數(shù). 列主元素消去法的步驟：設(shè)已經(jīng)完成第1步到第步的按列選主元、交換兩行、消元計(jì)算，得到矩陣 . 第步計(jì)算如下：對(duì)于， （1） 選列主元素，即確定使； （2） 如果，則方程組解不唯一，或者接近奇異矩陣，停止運(yùn)算； （3） 如果，則交換第行與第行元素； （4） 消元計(jì)算： （5） 回代計(jì)算： 完全主元素消去法即是每次選主元時(shí)，依次按行、列選取絕對(duì)值最大的元素作為主元素，然后交換兩行、兩列，再進(jìn)行消元計(jì)算. 完全主元素消去法的步驟：設(shè)已經(jīng)完成第1步到第步的選主元、交換行和列、消元計(jì)算，得到矩陣 . 第步計(jì)算選主元素的范圍為，即確定使. 第步計(jì)算如下：對(duì)于， （1） 選主元素，即確定使； （2） 如果，則方程組解不唯一，或者接近奇異矩陣，停止運(yùn)算； （3） 如果，則交換第行與第行元素；如果，則交換第列與第列元素； （4） 消元計(jì)算： （5） 回代求解. 【注】 完全主元消去法是解低階稠密矩陣方程組的有效方法，但完全主元消去法解方程組，在選主元素時(shí)要化費(fèi)較多的計(jì)算機(jī)時(shí)間，行主元消去法與列主元消去法運(yùn)算量大體相同，實(shí)際計(jì)算時(shí)，用列主元消去法即可滿足一定的精度要求. 對(duì)同一數(shù)值問(wèn)題，用不同的計(jì)算方法，所得結(jié)果的精度大不一樣.對(duì)于一個(gè)算法來(lái)說(shuō)，如果計(jì)算過(guò)程中舍入誤差能得到控制，對(duì)計(jì)算結(jié)果影響較小，則稱此算法是數(shù)值穩(wěn)定的；否則，如果計(jì)算過(guò)程中舍入誤差增長(zhǎng)迅速，計(jì)算結(jié)果受舍入誤差影響較大，則稱此算法為數(shù)值不穩(wěn)定的.因此，我們解數(shù)值問(wèn)題時(shí)，應(yīng)選擇和使用數(shù)值穩(wěn)定的算法，否則如果使用數(shù)值不穩(wěn)定的算法，就可能導(dǎo)致計(jì)算失敗. 例 用高斯列主元消去法解方程組 解 . 所以，方程組的解為. ４）高斯列主元素消去法的MATLAB實(shí)現(xiàn)：，意為． 例 LinearEquiation02.m open LinearEquiation02 LinearEquiation02 一個(gè)典型的例子: Hilbert矩陣: 注: 非奇異矩陣的條件數(shù): ５）ＬＵ分解?。ǎ蹋铡actorization）（高斯消去法、Doolittle分解） 　　高斯消去法的消元過(guò)程，從代數(shù)運(yùn)算的角度看就是用一個(gè)下三角矩陣左乘方程組的系數(shù)矩陣Ａ，且乘積的結(jié)果為上三角矩陣，即 (２) 可通過(guò)直接用A元素計(jì)算矩陣A的三角分解矩陣L和U.這種直接計(jì)算A的三角分解的方法有實(shí)用上的好處.下面利用矩陣乘法規(guī)則來(lái)確定三角矩陣L和U. . 第一步：利用A的第一行、第一列元素確定U的第一行、L的第一列元素.由矩陣乘法， ， ， 得到 ,. （3.7） 設(shè)已經(jīng)計(jì)算出U的第1至r -1行元素，L的第1至r -1列元素，現(xiàn)在要計(jì)算U的第r行元素及L的第r列元素. 第r步：利用A的第r行、第r列剩下的元素確定U的第r行、L的第r列元素.由矩陣乘法，有 ， 得U的第r行元素為 . （3.8） 由 ， 得 . （3.9） 例5 用LU分解法求解方程組 . 解 對(duì)系數(shù)矩陣A進(jìn)行LU分解， . 由，有. ，. 因此 . 解方程組，得. 解方程組，得. 6） LU 分解的MATLAB實(shí)現(xiàn)：或 例 A=rand(5); [L,U，P]=lu(A) A=rand(5); [L,U,P]=lu(A) L=P\L 當(dāng)是主對(duì)角占優(yōu)的三對(duì)角矩陣時(shí)，基于Doolittle分解可得到解這類方程組的追趕法。 ２、Cholesky分解 (Cholesky Factorization) 　　對(duì)稱正定矩陣的Cholesky分解和以為系數(shù)矩陣地的線性方程組的改進(jìn)的平方根法： 設(shè)階方程組，是對(duì)稱正定矩陣(Positive Definite Matrix)，則有三角分解 . 再將分解為 , 則. （1） 對(duì)稱正定矩陣有唯一的分解 這是由于，，且對(duì)稱陣，則有 再利用三角分解的唯一性，得.因此，對(duì)稱正定矩陣有唯一的分解. （2） 是正定對(duì)角陣（即） 由于對(duì)稱正定的充要條件是對(duì)稱正定，其中是階可逆方陣.取，就推知是正定對(duì)角陣. 因此的對(duì)角元素，記，其中 , 則. （3） 喬萊斯基（Cholesky）分解 將記為，則稱為Cholesky分解.利用Cholesky直接分解公式，推導(dǎo)出的解方程組方法，稱為Cholesky方法或平方根法. （4） 解方程組的平方根法（Cholesky方法） 由Cholesky分解，有 . （3.10） 利用矩陣乘法，逐步確定的第行元素. 由（當(dāng)時(shí)，），有分解公式：對(duì)于 （3.11） 將對(duì)稱正定矩陣作Cholesky分解后，則解方程組就轉(zhuǎn)化為解兩個(gè)三角方程組. 例7 用Cholesky方法解方程組 . 解 對(duì)系數(shù)矩陣作Cholesky分解得到 . 解，得. 解，得. cholesky分解的MATLAB的實(shí)現(xiàn)：L=chol(A)。 3、追趕法 在許多實(shí)際問(wèn)題中，如，常微分方程兩點(diǎn)邊值問(wèn)題、三次樣條插值方法等，往往遇到線性方程組的求解，其中 . （3.13） 稱具有公式（3.13）形式的系數(shù)矩陣為三對(duì)角陣,稱相應(yīng)的線性方程組為三對(duì)角方程組（Tridiagonal Linear Systems).具有這種形式的方程組在實(shí)際問(wèn)題中是經(jīng)常遇到的，而且往往是對(duì)角占優(yōu)（Diagonally Dominant）的.滿足條件： ， ， . 這類方程組的解存在唯一（非奇異），可以直接利用高斯消去法或直接分解法，而其解答可以用極其簡(jiǎn)單的遞推公式表示出來(lái)，即下面介紹的追趕法.追趕法通常是數(shù)值穩(wěn)定的. 對(duì)作LU分解（Doolitle分解），可以發(fā)現(xiàn)L、U具有非常簡(jiǎn)單的形式. . 由矩陣乘積，得 . 比較等式兩端，得到 （3.14） 因?yàn)樯鲜龇纸猓瑒t方程組的求解轉(zhuǎn)化為解兩個(gè)簡(jiǎn)單的三角方程組和，從而得到求解方程組的算法公式. 先解，即 . （3.15） 再解，即 . （3.16） 這種把三對(duì)角方程組的解用遞推公式（3.14）、（3.15）、（3.16）表示出來(lái)的方法形象化地叫做追趕法，其中（3.14）、（3.15）是關(guān)于下標(biāo)由小到大的遞推公式稱為追的過(guò)程，而（16）卻是下標(biāo)由大到小的遞推公式稱為趕的過(guò)程，一追一趕構(gòu)成了求解的追趕法. 例9 用追趕法解三對(duì)角方程組 . 解 系數(shù)矩陣分解得到 . 解，得. 解，得. 調(diào)用函數(shù)LU_Factorization.m解例9.輸入 A=[4 -1 0;-1 4 -1;0 -1 4];b=[1;3;2];[x,L,U,index]=LU_Factorization(A,b) 得到方程組的解及相應(yīng)的LU分解矩陣： x = 0.5179 L= 1.0000 0 0 U= 4.0000 -1.0000 0 1.0714 -0.2500 1.0000 0 0 3.7500 -1.0000 0.7679 0 -0.2667 1.0000 0 0 3.7333 為了對(duì)線性方程組的直接法作出誤差分析，為了討論方程組迭代法的收斂性，需要對(duì)向量和矩陣的大小進(jìn)行度量，進(jìn)而引入了 　　范數(shù)━ 用于度量“量”的大小的概念 引言 實(shí)數(shù)的絕對(duì)值：是數(shù)軸上的點(diǎn)到原點(diǎn)的距離； 復(fù)數(shù)的模：是平面上的點(diǎn)到原點(diǎn)的距離； 還有其他刻畫復(fù)數(shù)大小的方法（準(zhǔn)則）：如 １）； ２） 　向量的內(nèi)積、范數(shù)及維空間距離的度量 令是一數(shù)域，是上的向量空間，如果函數(shù)有如下性質(zhì)： １、共軛對(duì)稱性：，； ２、非負(fù)性：，，； ３、線性性：，， ； 則稱是上的一個(gè)向量?jī)?nèi)積（inner product），向量空間上的向量?jī)?nèi)積通常用符號(hào)表示，定義了內(nèi)積的向量空間稱為內(nèi)積空間（inner product space）。記做表示。 例１．１　，，，容易驗(yàn)證函數(shù) 　　　　　　　　　（１．１） 定義了上的一個(gè)內(nèi)積。 令是一數(shù)域，是上的向量空間，如果函數(shù)有如下性質(zhì)： １、非負(fù)性：，，； ２、齊次性：，，； ３、三角不等式：，； 則稱是上的一個(gè)向量范數(shù)（norm），向量空間上的范數(shù)通常用符號(hào)表示。定義了范數(shù)的向量空間稱為賦范空間（normed space）。記做表示。 例１．２　，，，容易驗(yàn)證函數(shù) 　　 　　　　?。ǎ保玻? 定義了上的一個(gè)范數(shù)，這樣定義的范數(shù)稱為由內(nèi)積（１．１）誘導(dǎo)的范數(shù)。 例１．３　上常用的向量范數(shù)： ， １、１—范數(shù)：； ２、２—范數(shù)：； ３、—范數(shù)：； 令是一數(shù)域，是上的向量空間，如果實(shí)值函數(shù)有如下性質(zhì)： １、對(duì)稱性：，； ２、非負(fù)性：，， ３、三角不等式：，； 則稱是上的一個(gè)距離（函數(shù)）（distance function）或度量（metric），定義了度量的向量空間稱為度量空間（metric space），記做表示。 例１．4　上常用的（由范數(shù)誘導(dǎo)的）度量：， １、１—范數(shù)誘導(dǎo)的度量：； ２、２—范數(shù)誘導(dǎo)的度量：； ３、—范數(shù)誘導(dǎo)的度量：； 矩陣的范數(shù) 矩陣是線性映射（當(dāng)時(shí)為線性變換）的一種表現(xiàn)形式。因此，除了可以把矩陣看做向量而定義其范數(shù)外，更為基本、更為重要的是表征其線性映射的算子范數(shù)（operator norm），以的情況為例： 　　　　　　　　　?。ǎ保常? 其中（１．３）右端的范數(shù)是賦范空間中向量的范數(shù)，由矩陣算子范數(shù)的定義（１．３）容易證明（對(duì)映像大小的估計(jì)）不等式： ，　， ?。ǎ保矗? 稱滿足不等式（１．４）的矩陣范數(shù)是與對(duì)應(yīng)的向量范數(shù)相容的。 例１．５　常用的矩陣范數(shù)： １、１—范數(shù)（列范數(shù)）： ； ２、２—范數(shù)（譜范數(shù)）： ； ３、—范數(shù)（行范數(shù)）： ； 上述三種范數(shù)是如下定義的矩陣—范數(shù)的特例： ４、由向量的—范數(shù)：，，定義： 　　　　　　　?。ǎ保担? ５、Ｆ—范數(shù)（Frobenius）: ； 例 設(shè)，求范數(shù). 解 . . . 由于，所以，因此. 向量和矩陣范數(shù)的MATLAB實(shí)現(xiàn) 在MATLAB中，norm( )函數(shù)求向量與矩陣的范數(shù)，其命令格式為norm(X,p). 當(dāng)X為向量或矩陣時(shí)，norm(X,p)表示向量或矩陣X的p范數(shù)，例如，norm(X,1)表示X的1范數(shù)，norm(X,2)表示X的2范數(shù)，norm(X,inf)表示X的范數(shù).對(duì)于矩陣X，norm(X,fro)表示X的F范數(shù). 矩陣的條件數(shù)與直接法的誤差分析 解線性方程組的直接法產(chǎn)生誤差的主要原因：1）不同的算法及舍入誤差的影響；2）方程組本身固有的問(wèn)題（病態(tài)或良態(tài)）本節(jié)我們將分析方程組的狀態(tài)并估計(jì)算法的誤差，即原始數(shù)據(jù)擾動(dòng)對(duì)解的影響. 考慮階線性方程組，其中為非奇異矩陣. 由于（或）的數(shù)值是測(cè)量得到的，或者是計(jì)算的結(jié)果，在第一種情況下（或）常帶有某些觀測(cè)誤差，在后一種情況（或）又包含有舍入誤差.因此我們處理的實(shí)際矩陣是（或），下面我們來(lái)研究數(shù)據(jù)（或）的微小誤差對(duì)解的影響.首先考慮一個(gè)例子. 例 設(shè)方程組，即，它的精確解為. 現(xiàn)在考慮系數(shù)矩陣和右端項(xiàng)的微小變化對(duì)方程組解的影響，即考察方程組 ， 其解變?yōu)?這種現(xiàn)象的出現(xiàn)是完全由方程組的性態(tài)決定的. 定義：如果矩陣或常數(shù)項(xiàng)的微小變化，引起方程組解的巨大變化，則稱此方程組為病態(tài)方程組，矩陣稱為病態(tài)矩陣（相對(duì)于方程組而言），否則稱方程組為良態(tài)方程組，矩陣稱為良態(tài)矩陣.需要一種能刻劃矩陣和方程組“病態(tài)”程度的標(biāo)準(zhǔn). 線性方程組的誤差分析 設(shè)線性方程組為 ， 其中，，且非奇異.為精確解，為解的誤差，記. 設(shè)為的誤差，為的誤差.下面分別討論與，的關(guān)系. b有誤差而A無(wú)誤差情形 將帶有誤差的右端項(xiàng)和帶誤差的解向量代入方程組，則 . 由于，而得到，從而. 另一方面，由(3.24)式取范數(shù)，有，即.可得 定理 設(shè)是非奇異矩陣，，且，則有誤差估計(jì)式 , 其中稱為方陣A的條件數(shù). 【注】 （1） 解的相對(duì)誤差是右端項(xiàng)的相對(duì)誤差的cond(A)倍； （2） 如果條件數(shù)越大，則解的相對(duì)誤差就可能越大； （3） 條件數(shù)成了刻劃矩陣的病態(tài)程度和方程組解對(duì)或擾動(dòng)的敏感程度. 【定義7】 稱條件數(shù)很大的矩陣為“病態(tài)”矩陣；稱病態(tài)矩陣對(duì)應(yīng)的方程組為病態(tài)方程組.反之，則稱矩陣為良態(tài)矩陣，對(duì)應(yīng)的方程組為良態(tài)方程組. A及b都有誤差的情形 定理7 設(shè)方程組，為非奇異矩陣，及都有誤差，若的誤差非常小，使，則有誤差估計(jì)式 . (3.26) 例 已知方程組中 ， 若時(shí)，估計(jì)解的相對(duì)誤差. 解 由于由式（3.25）有誤差估計(jì) ， 比右端項(xiàng)的相對(duì)誤差擴(kuò)大了2015倍. 例 下列希爾伯特（Hilbert）矩陣是一族著名的病態(tài)矩陣. . 用MATLAB函數(shù)計(jì)算條件數(shù).輸入 for n=3:8 cond(hilb(n)) end 得到3至8階希爾伯特矩陣的條件數(shù)分別為： 524.0568 1.5514e+004 4.7661e+005 1.4951e+007 4.7537e+008 1.5258e+010 由此可見(jiàn)，隨著的增加，的病態(tài)可能越嚴(yán)重.常出現(xiàn)在數(shù)據(jù)擬合和函數(shù)逼近中. 對(duì)于病態(tài)方程組，通常的方法無(wú)法得到它的準(zhǔn)確解，需要采用一些特殊的處理方法. 病態(tài)方程組的處理 對(duì)于病態(tài)方程組，可采用高精度的算術(shù)運(yùn)算，如雙精度或擴(kuò)充精度，以改善或減輕方程組的病態(tài)程度.如果用無(wú)限精度運(yùn)算即不存在舍入誤差的話，即使條件數(shù)很大，也沒(méi)有病態(tài)可言.我們也可對(duì)病態(tài)方程組作預(yù)處理，使改善方程組系數(shù)矩陣的條件數(shù). 例 設(shè)方程組，即 ， 試驗(yàn)證其為病態(tài)方程組，且對(duì)其作預(yù)處理，使. 解 （1）用MATLAB函數(shù)計(jì)算系數(shù)矩陣的條件數(shù)，得，顯然方程組為病態(tài)方程組. （2） 令，使，其中 ， 則有.顯然，經(jīng)過(guò)預(yù)處理后的方程組是良態(tài)的. 奇異值分解法解病態(tài)方程組. 奇異值分解（Singular-Value Decomposition）簡(jiǎn)稱SVD，對(duì)于階矩陣，必存在正交矩陣，和對(duì)角陣，使得有奇異值分解 . （3.27） 奇異值分解非常有用，對(duì)于矩陣階矩陣，存在階正交矩陣，階對(duì)角陣，滿足。和中分別是的奇異向量，而是的奇異值。的正交單位特征向量組成，特征值組成，的正交單位特征向量組成，特征值（與相同）組成。因此，奇異值分解和特征值問(wèn)題緊密聯(lián)系（注：奇異值分解對(duì)也適用）。 在MATLAB中，函數(shù)svd( )表示矩陣的奇異值分解，其命令格式為 [U,S,V]=svd(A) 其中，，為正交矩陣，為對(duì)角陣.例如，求4階Hilbert矩陣的奇異值分解，輸入 [U,S,V]=svd(hilb(4)) 得到 U = -0.7926 0.5821 -0.1792 -0.0292 -0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146 S = 1.5002 0 0 0 0 0.1691 0 0 0 0 0.0067 0 0 0 0 0.0001 V = -0.7926 0.5821 -0.1792 -0.0292 -0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146 矩陣對(duì)角線上的元素稱為奇異值，由大到小排列.