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

《數(shù)值分析 第三章 基于MATLAB的科學計算—線性方程組1》由會員分享，可在線閱讀，更多相關(guān)《數(shù)值分析 第三章 基于MATLAB的科學計算—線性方程組1（22頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、科學計算—理論、方法 及其基于ＭＡＴＬＡＢ的程序?qū)崿F(xiàn)與分析 三、 解線性方程組（線性矩陣方程） 解線性方程組是科學計算中最常見的問題。所說的“最常見”有兩方面的含義： １） 問題的本身是求解線性方程組； ２） 許多問題的求解需要或歸結(jié)為線性方程組的求解。 　　關(guān)于線性方程組 　　　　　　（１） 其求解方法有兩類： １） 直接法：高斯消去法（Gaussian Elimination）； ２） 間接法：各種迭代法（Iteration）。 １、高斯消去法 １） 引例 考慮如下（梯形）線性方程組： 高斯消去法的求解思路：把一般的線性方程組（１）化成（上或下）梯形的形式。

2、 ２）高斯消去法——示例 考慮如下線性方程組： １） 第一個方程的兩端乘加到第二個方程的兩端，第一個方程的兩端乘 －１加到第三個方程的兩端，得 ２）　第二個方程的兩端乘加到第三個方程的兩端，得 ３） 從上述方程組的第三個方程依此求解，得 ３）高斯消去法的不足及其改進——高斯（全、列）主元素消去法 在上例中，由于建模、計算等原因，系數(shù)2.001而產(chǎn)生0.0005的誤差，實際求解的方程組為 注：數(shù)值穩(wěn)定的算法 高斯列主元素消去法就是在消元的每一步選?。校┲髟亍涣兄薪^對值最大的元取做主元素，高斯列主元素消去法是數(shù)值穩(wěn)定的方法。 列主元素消去法的基本思想

3、：在每輪消元之前，選列主元素（絕對值最大的元素）,使乘數(shù). 列主元素消去法的步驟：設(shè)已經(jīng)完成第1步到第步的按列選主元、交換兩行、消元計算，得到矩陣 . 第步計算如下：對于， （1） 選列主元素，即確定使； （2） 如果，則方程組解不唯一，或者接近奇異矩陣，停止運算； （3） 如果，則交換第行與第行元素； （4） 消元計算： （5） 回代計算： 完全主元素消去法即是每次選主元時，依次按行、列選取絕對值最大的元素作為主元素，然后交換兩行、兩列，再進行消元計算. 完全主元素消去法的步驟：設(shè)已經(jīng)完成第1步到第步的選主元、交換行和列、消元計算，得到矩陣 . 第步計算選主

4、元素的范圍為，即確定使. 第步計算如下：對于， （1） 選主元素，即確定使； （2） 如果，則方程組解不唯一，或者接近奇異矩陣，停止運算； （3） 如果，則交換第行與第行元素；如果，則交換第列與第列元素； （4） 消元計算： （5） 回代求解. 【注】 完全主元消去法是解低階稠密矩陣方程組的有效方法，但完全主元消去法解方程組，在選主元素時要化費較多的計算機時間，行主元消去法與列主元消去法運算量大體相同，實際計算時，用列主元消去法即可滿足一定的精度要求. 對同一數(shù)值問題，用不同的計算方法，所得結(jié)果的精度大不一樣.對于一個算法來說，如果計算過程中舍入誤差能得到控制，對計算結(jié)果

5、影響較小，則稱此算法是數(shù)值穩(wěn)定的；否則，如果計算過程中舍入誤差增長迅速，計算結(jié)果受舍入誤差影響較大，則稱此算法為數(shù)值不穩(wěn)定的.因此，我們解數(shù)值問題時，應選擇和使用數(shù)值穩(wěn)定的算法，否則如果使用數(shù)值不穩(wěn)定的算法，就可能導致計算失敗. 例 用高斯列主元消去法解方程組 解 . 所以，方程組的解為. ４）高斯列主元素消去法的MATLAB實現(xiàn)：，意為． 例 LinearEquiation02.m open LinearEquiation02 LinearEquiation02 一個典型的例子: Hilbert矩陣: 注: 非奇異矩陣的條件數(shù): ５）ＬＵ分解　

6、（ＬＵ　Factorization）（高斯消去法、Doolittle分解） 　　高斯消去法的消元過程，從代數(shù)運算的角度看就是用一個下三角矩陣左乘方程組的系數(shù)矩陣Ａ，且乘積的結(jié)果為上三角矩陣，即 (２) 可通過直接用A元素計算矩陣A的三角分解矩陣L和U.這種直接計算A的三角分解的方法有實用上的好處.下面利用矩陣乘法規(guī)則來確定三角矩陣L和U. . 第一步：利用A的第一行、第一列元素確定U的第一行、L的第一列元素.由矩陣乘法， ， ， 得到 ,. （3.7） 設(shè)已經(jīng)計算出U的第1至r -1行元素，L的第1至r -1列元素，現(xiàn)在要計算U的第r

7、行元素及L的第r列元素. 第r步：利用A的第r行、第r列剩下的元素確定U的第r行、L的第r列元素.由矩陣乘法，有 ， 得U的第r行元素為 . （3.8） 由 ， 得 . （3.9） 例5 用LU分解法求解方程組 . 解 對系數(shù)矩陣A進行LU分解， . 由，有. ，. 因此 . 解方程組，得. 解方程組，得. 6） LU 分解的MATLAB實現(xiàn)：或 例 A=rand(5); [L,U，P]=lu(A) A=rand(5); [L,U,P]=lu(A) L=P\L 當是主對角占優(yōu)的三對角矩陣時，基于Dooli

8、ttle分解可得到解這類方程組的追趕法。 ２、Cholesky分解 (Cholesky Factorization) 　　對稱正定矩陣的Cholesky分解和以為系數(shù)矩陣地的線性方程組的改進的平方根法： 設(shè)階方程組，是對稱正定矩陣(Positive Definite Matrix)，則有三角分解 . 再將分解為 , 則. （1） 對稱正定矩陣有唯一的分解 這是由于，，且對稱陣，則有 再利用三角分解的唯一性，得.因此，對稱正定矩陣有唯一的分解. （2） 是正定對角陣（即） 由于對稱正定的充要條件是對稱正定，其中是階可逆方陣.取，就推知是正定對角陣. 因此的對角元素，

9、記，其中 , 則. （3） 喬萊斯基（Cholesky）分解 將記為，則稱為Cholesky分解.利用Cholesky直接分解公式，推導出的解方程組方法，稱為Cholesky方法或平方根法. （4） 解方程組的平方根法（Cholesky方法） 由Cholesky分解，有 . （3.10） 利用矩陣乘法，逐步確定的第行元素. 由（當時，），有分解公式：對于 （3.11） 將對稱正定矩陣作Cholesky分解后，則解方程組就轉(zhuǎn)化為解兩個三角方程組. 例7 用Cholesky方法解方程組 . 解 對系數(shù)矩陣作Cholesky分解得到 . 解

10、，得. 解，得. cholesky分解的MATLAB的實現(xiàn)：L=chol(A)。 3、追趕法 在許多實際問題中，如，常微分方程兩點邊值問題、三次樣條插值方法等，往往遇到線性方程組的求解，其中 . （3.13） 稱具有公式（3.13）形式的系數(shù)矩陣為三對角陣,稱相應的線性方程組為三對角方程組（Tridiagonal Linear Systems).具有這種形式的方程組在實際問題中是經(jīng)常遇到的，而且往往是對角占優(yōu)（Diagonally Dominant）的.滿足條件： ， ， . 這類方程組的解存在唯一（非奇異），可以直接利用高斯消去法或直接分解法，而其

11、解答可以用極其簡單的遞推公式表示出來，即下面介紹的追趕法.追趕法通常是數(shù)值穩(wěn)定的. 對作LU分解（Doolitle分解），可以發(fā)現(xiàn)L、U具有非常簡單的形式. . 由矩陣乘積，得 . 比較等式兩端，得到 （3.14） 因為上述分解，則方程組的求解轉(zhuǎn)化為解兩個簡單的三角方程組和，從而得到求解方程組的算法公式. 先解，即 . （3.15） 再解，即 . （3.16） 這種把三對角方程組的解用遞推公式（3.14）、（3.15）、（3.16）表示出來的方法形象化地叫做追趕法，其中（3.14）、（3.15

12、）是關(guān)于下標由小到大的遞推公式稱為追的過程，而（16）卻是下標由大到小的遞推公式稱為趕的過程，一追一趕構(gòu)成了求解的追趕法. 例9 用追趕法解三對角方程組 . 解 系數(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) 得到方程組的解及相應的LU分解矩陣： x = 0.5179 L= 1.0000 0 0 U= 4.0000 -1.0000

13、 0 1.0714 -0.2500 1.0000 0 0 3.7500 -1.0000 0.7679 0 -0.2667 1.0000 0 0 3.7333 為了對線性方程組的直接法作出誤差分析，為了討論方程組迭代法的收斂性，需要對向量和矩陣的大小進行度量，進而引入了 　　范數(shù)━ 用于度量“量”的大小的概念 引言 實數(shù)的絕對值：是數(shù)軸上的點到原點的距離； 復數(shù)的模：是平面上的點到原點的距離； 還有其他刻畫復數(shù)大小的方法（準則）：如

14、１）； ２） 　向量的內(nèi)積、范數(shù)及維空間距離的度量 令是一數(shù)域，是上的向量空間，如果函數(shù)有如下性質(zhì)： １、共軛對稱性：，； ２、非負性：，，； ３、線性性：，， ； 則稱是上的一個向量內(nèi)積（inner product），向量空間上的向量內(nèi)積通常用符號表示，定義了內(nèi)積的向量空間稱為內(nèi)積空間（inner product space）。記做表示。 例１．１　，，，容易驗證函數(shù) 　　　　　　　　?。ǎ保保? 定義了上的一個內(nèi)積。 令是一數(shù)域，是上的向量空間，如果函數(shù)有如下性質(zhì)： １、非負性：，，； ２、齊次性：，，； ３、三角不等式：，； 則稱是上的一個向量范數(shù)（nor

15、m），向量空間上的范數(shù)通常用符號表示。定義了范數(shù)的向量空間稱為賦范空間（normed space）。記做表示。 例１．２　，，，容易驗證函數(shù) 　　 　　　　?。ǎ保玻? 定義了上的一個范數(shù)，這樣定義的范數(shù)稱為由內(nèi)積（１．１）誘導的范數(shù)。 例１．３　上常用的向量范數(shù)： ， １、１—范數(shù)：； ２、２—范數(shù)：； ３、—范數(shù)：； 令是一數(shù)域，是上的向量空間，如果實值函數(shù)有如下性質(zhì)： １、對稱性：，； ２、非負性：，， ３、三角不等式：，； 則稱是上的一個距離（函數(shù)）（distance function）或度量（metric），定義了度量的向量空間稱為度量空

16、間（metric space），記做表示。 例１．4　上常用的（由范數(shù)誘導的）度量：， １、１—范數(shù)誘導的度量：； ２、２—范數(shù)誘導的度量：； ３、—范數(shù)誘導的度量：； 矩陣的范數(shù) 矩陣是線性映射（當時為線性變換）的一種表現(xiàn)形式。因此，除了可以把矩陣看做向量而定義其范數(shù)外，更為基本、更為重要的是表征其線性映射的算子范數(shù)（operator norm），以的情況為例： 　　　　　　　　　?。ǎ保常? 其中（１．３）右端的范數(shù)是賦范空間中向量的范數(shù)，由矩陣算子范數(shù)的定義（１．３）容易證明（對映像大小的估計）不等式： ，　， ?。ǎ保矗? 稱滿足不等式（１．４）的矩陣范

17、數(shù)是與對應的向量范數(shù)相容的。 例１．５　常用的矩陣范數(shù)： １、１—范數(shù)（列范數(shù)）： ； ２、２—范數(shù)（譜范數(shù)）： ； ３、—范數(shù)（行范數(shù)）： ； 上述三種范數(shù)是如下定義的矩陣—范數(shù)的特例： ４、由向量的—范數(shù)：，，定義： 　　　　　　　　（１．５） ５、Ｆ—范數(shù)（Frobenius）: ； 例 設(shè)，求范數(shù). 解 . . . 由于，所以，因此. 向量和矩陣范數(shù)的MATLAB實現(xiàn) 在MATLAB中，norm( )函數(shù)求向量與矩陣的范數(shù)，其命令格式為norm(X,p). 當X為向量或矩陣時，norm(X,p)表示向量或矩陣X的p范數(shù)，例如，norm

18、(X,1)表示X的1范數(shù)，norm(X,2)表示X的2范數(shù)，norm(X,inf)表示X的范數(shù).對于矩陣X，norm(X,fro)表示X的F范數(shù). 矩陣的條件數(shù)與直接法的誤差分析 解線性方程組的直接法產(chǎn)生誤差的主要原因：1）不同的算法及舍入誤差的影響；2）方程組本身固有的問題（病態(tài)或良態(tài)）本節(jié)我們將分析方程組的狀態(tài)并估計算法的誤差，即原始數(shù)據(jù)擾動對解的影響. 考慮階線性方程組，其中為非奇異矩陣. 由于（或）的數(shù)值是測量得到的，或者是計算的結(jié)果，在第一種情況下（或）常帶有某些觀測誤差，在后一種情況（或）又包含有舍入誤差.因此我們處理的實際矩陣是（或），下面我們來研究數(shù)據(jù)（或）的微小誤差對

19、解的影響.首先考慮一個例子. 例 設(shè)方程組，即，它的精確解為. 現(xiàn)在考慮系數(shù)矩陣和右端項的微小變化對方程組解的影響，即考察方程組 ， 其解變?yōu)?這種現(xiàn)象的出現(xiàn)是完全由方程組的性態(tài)決定的. 定義：如果矩陣或常數(shù)項的微小變化，引起方程組解的巨大變化，則稱此方程組為病態(tài)方程組，矩陣稱為病態(tài)矩陣（相對于方程組而言），否則稱方程組為良態(tài)方程組，矩陣稱為良態(tài)矩陣.需要一種能刻劃矩陣和方程組“病態(tài)”程度的標準. 線性方程組的誤差分析 設(shè)線性方程組為 ， 其中，，且非奇異.為精確解，為解的誤差，記. 設(shè)為的誤差，為的誤差.下面分別討論與，的關(guān)系. b有誤差而A無誤差情形 將帶有誤差的右

20、端項和帶誤差的解向量代入方程組，則 . 由于，而得到，從而. 另一方面，由(3.24)式取范數(shù)，有，即.可得 定理 設(shè)是非奇異矩陣，，且，則有誤差估計式 , 其中稱為方陣A的條件數(shù). 【注】 （1） 解的相對誤差是右端項的相對誤差的cond(A)倍； （2） 如果條件數(shù)越大，則解的相對誤差就可能越大； （3） 條件數(shù)成了刻劃矩陣的病態(tài)程度和方程組解對或擾動的敏感程度. 【定義7】 稱條件數(shù)很大的矩陣為“病態(tài)”矩陣；稱病態(tài)矩陣對應的方程組為病態(tài)方程組.反之，則稱矩陣為良態(tài)矩陣，對應的方程組為良態(tài)方程組. A及b都有誤差的情形 定理7 設(shè)方程組，為非奇異矩陣，及都有誤差，

21、若的誤差非常小，使，則有誤差估計式 . (3.26) 例 已知方程組中 ， 若時，估計解的相對誤差. 解 由于由式（3.25）有誤差估計 ， 比右端項的相對誤差擴大了2015倍. 例 下列希爾伯特（Hilbert）矩陣是一族著名的病態(tài)矩陣. . 用MATLAB函數(shù)計算條件數(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 由此可

22、見，隨著的增加，的病態(tài)可能越嚴重.常出現(xiàn)在數(shù)據(jù)擬合和函數(shù)逼近中. 對于病態(tài)方程組，通常的方法無法得到它的準確解，需要采用一些特殊的處理方法. 病態(tài)方程組的處理 對于病態(tài)方程組，可采用高精度的算術(shù)運算，如雙精度或擴充精度，以改善或減輕方程組的病態(tài)程度.如果用無限精度運算即不存在舍入誤差的話，即使條件數(shù)很大，也沒有病態(tài)可言.我們也可對病態(tài)方程組作預處理，使改善方程組系數(shù)矩陣的條件數(shù). 例 設(shè)方程組，即 ， 試驗證其為病態(tài)方程組，且對其作預處理，使. 解 （1）用MATLAB函數(shù)計算系數(shù)矩陣的條件數(shù)，得，顯然方程組為病態(tài)方程組. （2） 令，使，其中 ， 則有.顯然，經(jīng)過預處理

23、后的方程組是良態(tài)的. 奇異值分解法解病態(tài)方程組. 奇異值分解（Singular-Value Decomposition）簡稱SVD，對于階矩陣，必存在正交矩陣，和對角陣，使得有奇異值分解 . （3.27） 奇異值分解非常有用，對于矩陣階矩陣，存在階正交矩陣，階對角陣，滿足。和中分別是的奇異向量，而是的奇異值。的正交單位特征向量組成，特征值組成，的正交單位特征向量組成，特征值（與相同）組成。因此，奇異值分解和特征值問題緊密聯(lián)系（注：奇異值分解對也適用）。 在MATLAB中，函數(shù)svd( )表示矩陣的奇異值分解，其命令格式為 [U,S,

24、V]=svd(A) 其中，，為正交矩陣，為對角陣.例如，求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

25、 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 矩陣對角線上的元素稱為奇異值，由大到小排列.