《復旦大學Matlab課件—— (9).ppt》由會員分享，可在線閱讀，更多相關《復旦大學Matlab課件—— (9).ppt（22頁珍藏版）》請在匯文網(wǎng)上搜索。

1、Matlab Math線性方程組Cleve Morler著陳文斌（)1 Ax=b 有x=A-1b，但實際上并不顯式求A-1線性方程組例子：7x=21 x=21/7=3如果求逆 x=7-1 21=.142857 21=2.99997這就需要一次除和一次乘，且精度更低2Backslash運算符 AX=BX=AB 左除 XA=BX=B/A 右除33-by-3的例子4算法矩陣表示單位下三角陣上三角陣置換陣5置換陣PA交換A的行，AP交換A的列置換向量 p=4 1 3 2PA=A(p,:)Px=b x=PTb置換方程求解6三角陣的求解Ux=b x=zeros(n,1);for k=n:-1:1 x(k)

2、=b(k)/U(k,k);i=(1:k-1);b(i)=b(i)x(k)*U(i,k);end x=zeros(n,1);for k=n:-1:1 j=k+1:n x(k)=(b(k)U(k,j)*x(j)/U(k,k);end 求出x(n)，然后把x(n)消掉7LU分解C.F.Gauss:高斯消去法(GE)19551977:pivot選主元，舍入誤差的影響高斯消去法分為兩步：向前消去和向后替換8LU分解例子9線性方程組求解Ax=b LU=PA Ly=Pb Ux=y x=Ab10舍去和舍入六位有效數(shù)字：舍去浮點表示（chopped floating-pointing representatio

3、n）舍入浮點表示（rounded floating-pointing representation）11 為何要選主元？5位有效數(shù)字在消元時第二個方程乘上2.51036.001 2.5 103=1.50025 104 舍入為1.5002 104,再加上2.5,再一次舍入，成為 1.5004 104，則第三個方程為-0.001x2+(6)(0.99993)=6.001x3=1.5004 104/(1.5005 104)=0.99993 1.5005 104 x3=1.5004 104 x2=-1.5x1=-0.35真解：(0,-1,1)T12部分選主元如果每次乘數(shù)在量上是小于或等于1的，計算出來

4、的解就是精確的。要使乘數(shù)小于或等于1，可以通過部分選主元來實現(xiàn)，選要消去列的最大元素到對角元，交換A和b的行（x不要交換）。lutx.m:matlab內置lu分解的一個版本 bslashtx.m：matlab內置backslash運算 lugui.m：動態(tài)顯示LU分解13舍入誤差的影響設計算的解為x*,真解為x=A-1b,誤差為：e=x x*,殘量為 r=b Ax*六位精度算法 (選主元）得到精確解三位精度算法(選主元）：14分析在三位精度的算法中，x2是由兩個與舍入誤差相同階的數(shù)相除得到的，因此x2可以是任意的數(shù)。然后完全任意的x2代入第一個方程得到x1，所以第一方程的殘量是小的。我們也可以

5、期望第二個方程的殘量也是小的，因為矩陣是接近奇異的。部分主元的高斯消去法可以確保產(chǎn)生小的殘量部分主元的高斯消去法可以確保產(chǎn)生小的殘量“小的”舍入誤差的階與三個數(shù)量有關：原來矩陣悉數(shù)矩陣單元的大小，消去過程中系數(shù)矩陣單元的大小和計算解的元素的大小。如果其中一個是“大的”，殘量在絕對意義上不必是小的。即使殘量是小的，但不能保證解的誤差是小的。它們之間的關系可以部分由矩陣的條件數(shù)條件數(shù)來刻畫15范數(shù)如果Ax=b,如何測量x關于A和b的敏感性？x=(1:4)/5 n1=norm(x,1),n2=norm(x)ninf=norm(x,inf)16條件數(shù)條件數(shù)：17擾動分析條件數(shù)是相對誤差放大因子，也是刻

6、畫矩陣奇異性的測度。18例子norm119例子假設矩陣的第一個元素a11=0.1,所有其他元素和b是整數(shù)，且A的條件數(shù)是1012。假設用IEEE雙精度53位浮點算法，且我們假設可以精確求解計算機中的系統(tǒng)，那么由唯一的誤差是由0.1的二進制表示引起的，我們可以預期也就是說，簡單的矩陣系數(shù)的存儲就可能引起四位有效位的喪失。20矩陣范數(shù)J.H.Wilkinson 證明高斯消去法的計算結果x*精確滿足(A+E)x*=b21矩陣范數(shù) cond(A)或cond(A,2):使用svd(A)cond(A,1)計算k1(A)：使用inv(A)，工作量比cond(A,2)小 cond(A,inf):同cond(A,1)condest(A)估計k1(A)：用lu(A)用Higham和Tisseur的算法，適合大的稀疏矩陣 rcound(A)估計1/k1(A)：用lu(A)和一個在LINPACK和LAPACK工程中的一個老算法（歷史原因）22