解線性方程組的直接方法

1、第2章 解線性方程組的直接法,考慮一個簡單的例。保持人體健康需要多種微量元素，主要通過食物攝取，并且不同的年齡，職業(yè)和身體條件要求的量也不同。,考慮 n元線性方程組,(2.1),方程組(2.1)的矩陣形式為,Ax=b,其中,若矩陣A非奇異,即det(A)≠0,則方程組(2.1)有唯一解。,所謂直接解法是指，若不考慮計算過程中的舍入誤差，,經(jīng)過有限次算術(shù)運算就能求出線性方程組精確解的方法。,但由于實際計算中舍入誤差的存在，用直接解法一

2、般也只能求出方程組的近似解。 在線性代數(shù)中我們學過用Cramer法則求解線性方程組。,,行列式,而每個n階行列式按定義D=,需要計算(n-1)n!次乘法,則Cramer法則至少需要(n2-1)n!,次乘法,當n=20時,有(202-1)20!?9.7×1020次乘法運算。,如果用每秒鐘計算1百萬次乘除運算的計算機，約需要：,÷60,÷24,÷60,÷365,≈3000萬年,Cra

3、mer法則求n元線性方程組Ax=b的解,需要計算n+1個n階,9.7×1020÷106,可見Cramer法則不是一種實用的直接法，下面介紹幾種實,用的直接法。,§1 Gauss消去法,Gauss消元法是一種規(guī)則化的加減消元法，其基本思,想是通過逐次消元計算，把一般線性方程組的求解轉(zhuǎn)化為,等價的上三角形方程組的求解。,,,,對角,上三角,下三角,§1.1 順序Gauss消去法,為了清楚起見,先看

4、一個簡單的例子.,考慮線性方程組,,,單位上三角,,單位下三角,上三角方程組 Ax=b:,,,,,可按xn,xn-1,…x1，順序進行回代求解,消去后兩個方程中的x1得,再消去最后一個方程的x2得,消元結(jié)束,經(jīng)過回代得解:,,上述求解的消元過程可用矩陣表示為:,(A, b)=,這是Gauss消去法的矩陣計算形式,新的增廣矩陣對應的,線性方程組就是上三角形方程組,可進行回代求解。,現(xiàn)在介紹求解線性方程組 Ax=b 的順序Gauss

5、消去法:,記,則,線性方程組 Ax=b 的增廣矩陣為,,第一步.設(shè) ,依次用,乘矩陣的第1行加到第i行(i=2,3,…n),得到矩陣:,,其中,第二步.設(shè) ,依次用,乘矩陣(A(2),b(2))的第2行加到第i行,得到矩陣:,其中,,如此繼續(xù)消元下去,第n-1步結(jié)束后得到矩陣:,這就完成了消元過程。,對應的方程組變成：,對此方程組進行回代，就可求出方程組的解。,,順序Gauss消去法求解n元線性方程組的乘除運算量

6、是：,,n=20時,順序Gauss消去法只需3060次乘除法運算.,順序Gauss消去法通常也簡稱為Gauss消去法.,順序Gauss消去法中的,稱為主元素.,主元素都不為零?矩陣A的各階順序主子式都不為零.,,例 用高斯消去法求解Ax=b, 其中,,解,,,,回代求解：x3=3,x2=2,x1=2,§1.2 主元Gauss消去法,(用十進制四位浮點計算):,此問題的精確解為x1*=1.0001 ,x2*=0.9999。

7、,解 用順序Gauss消去法, 消元得,回代得解:x2=1.00 ,x1=0.00,若將方程組改寫成:,例1 解線性方程組,,再用Gauss消去法, 消元得,回代得解:x2=1.00 ,x1=1.00,為了提高計算的數(shù)值穩(wěn)定性,在消元過程中采用選擇主元的方法.常采用的是列主元消去法和全主元消去法.,給定線性方程組Ax=b, 記A(1)=A,b(1)=b,列主元Gauss消去法的具體過程如下:,首先在增廣矩陣B(1)=(A(1),b

8、(1))的第一列元素中,取,然后進行第一步消元得增廣矩陣B(2)=(A(2),b(2)).,再在矩陣B(2)=(A(2),b(2))的第二列元素中,取,,然后進行第二步消元得增廣矩陣B(3)=(A(3),b(3)).,按此方法繼續(xù)進行下去, 經(jīng)過n-1步選主元和消元運算,得到增廣矩陣B(n)=(A(n), b(n)).則方程組A(n)x = b(n)是與原方程組等價的上三角形方程組,可進行回代求解.,易證,只要|A|?0,列主元Gaus

9、s消去法就可順利進行.,采用十進制四位浮點計算,分別用順序Gauss消去法和列主元Gauss消去法求解線性方程組:,例2,,方程組具有四位有效數(shù)字的精確解為 x1*=17.46, x2*=-45.76, x3*=5.546,解 1. 用順序Gauss消去法求解,消元過程為,回代得: x3=5.546, x2=100.0, x1=-104.0,,2. 用列主元Gauss消去法求解,消元過程

10、為,,回代得: x3=5.545, x2=-45.77, x1=17.46,可見，列主元Gauss消去法是在每一步消元前,在主元所在的一列選取絕對值最大的元素作為主元素.,全主元Gauss消去法是在每一步消元前,在所有元素中選取絕對值最大的元素作為主元素.但由于運算量大增,實際應用中并不經(jīng)常使用.,,,§2 直接三角分解法,§2.1 Gauss消去法的矩陣表示,對矩陣,若,令,記,,則有,A(2)=,記,令,若,

11、則有,,A(3)=,如此進行下去, 第n-1步得到:,A(n)=,其中,,也就是:,A(n)=Ln-1A(n-1),其中,=Ln-1Ln-2A(n-2),=…= Ln-1Ln-2…L2L1A(1),,所以有:,A=A(1)= L1-1L2-1…Ln-1-1A(n)=LU,而且有,其中L=L1-1L2-1…Ln-1-1, U=A(n) .,A=LU, L為單位下三角矩陣; U是上三角矩陣.,,分解式 A=LU稱為矩陣A的直接三角分解.,

12、§2.2 Doolittle分解法,設(shè)n階方陣A的各階順序主子式不為零,則存在唯一單位下三角矩陣L和上三角矩陣U使A=LU .,證明 只證唯一性, 設(shè)有兩種分解 A=LU,則有,=E,所以得,設(shè) A有分解 A=LU,于是Ax=b ? LUx=b, 令 Ux=y, 則得,定理2.1,,這是與Ax=b 等價的線性方程組,下面介紹矩陣三角分解的Doolittle分解方法,,利用矩陣乘法規(guī)則依次

13、求 U的第k行，L的第k列(k=1,2…),,對k=2,3,…,n,計算,設(shè),akj=lk1u1j+lk2u2j+…+lkk-1uk-1j+ukj,aik=li1u1k+li2u2k+…+likukk,,,,對k=2,3,…,n,計算,要按順序先求U的第1行，L的第1列，U的第2行，L的第２列，…,U的第k行，L的第k列，…?？蓪С鲇嬎愎剑?在完成分解 A=LU后，依次求解方程組 Ly=b和 Ux=y，既,可得,這就是求解方程組Ax

14、=b的Doolittle三角分解方法。,,利用三角分解方法解線性方程組,解 令,利用矩陣乘法規(guī)則，先求 U的第1行元素，可得：,例3,,2=u112=u12-3=u13,,,,再求L的第2列。利用矩陣乘法規(guī)則建立方程得到：2=2l21,3=2l31,則有,,,再求U的第2行。利用矩陣乘法規(guī)則建立方程得到：,,,,-1=2+u22,,3=-3+u23,,則有: u22=-3, u23=6,,然后，再求U的第3行，2=-9/2+1

15、0+u33,最終得,,,,得,先解,再解,,得,解線性方程組Ax=b的Doolittle三角分解法的計算量約為1/3n3,與Gauss消去法基本相同.其優(yōu)點在于求一系列同系數(shù)的線性方程組Ax=bk ,(k=1,2,…,m)時,可大大節(jié)省運算量.,,求解：Ly=b和Ux=y,例如,求如下 3階矩陣A的逆矩陣：,這等價于求X使AX=E。設(shè)X=(x1,x2,x3)，可分別取3維單位向量 ： e1=(1,0,0)T, e2=(0,1,

16、0)T, e3=(0,0,1)T,由矩陣分塊乘法：Ax1=e1,Ax2=e2,Ax3=e3,利用三角分解法得:,依次求解 Lyi=ei,Uxi=yi,i=1,2,3.則 X=(x1,x2,x3)為所求.,所以,,,為了提高數(shù)值穩(wěn)定性, 可考慮列主元三角分解法, 設(shè)已完成A=LU的k-1步分解計算, 矩陣分解成,設(shè),令 rk?ri,相當于取,為第k步分解的主元素.,但要注意方程組的常數(shù)項也要相應變換.,,例如,用列主元三角分解解例3中方程

17、組.則有,,設(shè)A為對稱正定矩陣, 則有唯一分解A=LU, 且ukk>0.,則有 A=LDM,又因為 (LDM)T=MTDLT=LDM 所以 M=LT,=LDLT,則有,§2.3 平 方 根 法,,分解A=GGT稱為對稱正定矩陣的Cholesky分解.,,若記G=(gij), 則有: 對k=1,2,…,n,按列求矩陣 G的元素。實際計算時,可采用緊湊格式,,Ax=b?GGTx=b?Gy=

18、b , GTx=y, 稱為平方根法,Ax=b?GGTx=b?Gy=b , GTx=y, 解三角方程可得,解 令,例4 解線性方程組,,,平方根法是求對稱正定系數(shù)線性方程組的三角分解法,對稱正定矩陣的Cholesky分解的計算量和存貯量均約為一般矩陣的LU分解的一半. 且Cholesky分解具有數(shù)值穩(wěn)定性.,,利用矩陣乘法規(guī)則可計算得,追趕法是求三對角線性方程組的三角分解法.即方程組,三對角矩陣A的各階順序主子式都不為零的一個充

19、分條件是:,|a1|>|c1|>0 ; |an|>|dn|>0 ; |ai|?|ci|+|di| , cidi? 0 ,i=2,3,…,n-1.,在此條件下, A=LDM=TM , 稱之為矩陣A的Crout分解.,對三對角矩陣A進行Crout分解,有,§2.4 追 趕 法,,其中,Ax=b?TMx=b? Ty=b , Mx=y , 解三角方程可得,稱之為解三對角方程組的追趕法. 計算量為 0(

20、n).,,例5 解線性方程組,,解.令,利用矩陣乘法規(guī)則依次計算,得到,當滿足條件 |a1|>|c1|>0 ; |an|>|dn|>0 ; |ai|?|ci|+|di| , cidi? 0 ,i=2,3,…,n-1.時, 追趕法是數(shù)值穩(wěn)定的, 追趕法具有計算程序

21、簡單, 存貯量少,計算量小的優(yōu)點.,,,§3 向量和矩陣的范數(shù),§3.1 向量的范數(shù),定義2.1 設(shè)‖?‖是以向量 x?Rn 為自變量的實值函數(shù), 且滿足條件:,(1)非負性: 對任何向量x?Rn ,‖x‖?0 ,且‖x‖=0當 且僅當x=0,(2)齊次性: 對任何向量x ?Rn 和實數(shù)? , ‖?x‖=|? |‖x‖,(3)三角不等式: 對任何向量x

22、 ,y?Rn ‖x+y‖?‖x‖+‖y‖,則稱‖x‖為向量x的范數(shù)，‖?‖為空間Rn上的范數(shù)。,,記x=(x1,x2,…,xn)T, 常用的三種向量范數(shù)有:,向量的1-范數(shù):‖x‖1=|x1|+|x2|+…+|xn|,向量的2-范數(shù):‖x‖2=,向量的?-范數(shù):‖x‖?=,例6 設(shè)向量x=(2,-4,3,1)T, 求向量范數(shù)‖x‖p ,p=1,2, ?.,解 由定義‖x‖1=10 , ‖x‖2=,,‖x‖?

23、=4 .,雖然不同范數(shù)的值可能不同，但它們間存在等價關(guān)系.,定理2.2 (范數(shù)的等價性),對于 Rn 上的任何兩種范數(shù)‖?‖?和‖?‖? ,存在正常數(shù)m,M,使得 m ‖x‖?? ‖x‖? ?M ‖x‖? ,? x?Rn,,常用的三種向量范數(shù)滿足如下等價關(guān)系 ‖x‖?? ‖x‖1 ? n‖x‖? ,? x?Rn,定義2.2 設(shè)向量序列,k=1,2,…,向量,如果,則稱

24、向量序列{x(k)}收斂于向量x*, 記作,易見, 按|| x ||1范數(shù) ：,,§3.2 矩陣的范數(shù),定義2.3 設(shè)‖?‖是以n階方陣為自變量的實值函數(shù),且滿足條件:,(1)非負性: ‖A‖?0 ,且‖A‖=0當且僅當A=0,(2)齊次性: ‖?A‖=|? |‖A‖, ??R,(3)三角不等式:‖A+B‖?‖A‖+‖B‖,(4)乘積不等式:‖AB‖?‖A‖‖B‖,則稱‖A‖為矩陣A的范數(shù).,記A=(aij) , 常用的矩

25、陣范數(shù)有:,矩陣的1-范數(shù):‖A‖1,,也稱矩陣的列范數(shù).,矩陣的2-范數(shù):‖A‖2,,也稱為譜范數(shù).,,矩陣的?-范數(shù):‖A‖?,,也稱為行范數(shù).,矩陣的F-范數(shù):‖A‖F(xiàn),例7 設(shè)矩陣,求矩陣A的范數(shù)‖A‖p ,p=1,2, ?,F.,解 ‖A‖1=4 , ‖A‖?=5 , ‖A‖F(xiàn),,設(shè)‖?‖是一種向量范數(shù), 則定義,稱之為由向量范數(shù)派生的矩陣算子范數(shù).,矩陣的算子范數(shù)滿足,‖Ax‖?‖A‖‖x‖, ?x?Rn,把滿足上式的

26、矩陣范數(shù)稱為與向量范數(shù)相容的矩陣范數(shù).,對于p=1,2,?,矩陣范數(shù)‖A‖p是由向量范數(shù)‖x‖p派生的矩陣算子范數(shù),所以‖A‖p是與‖x‖p相容的矩陣范數(shù).但‖A‖F(xiàn)不是一種算子范數(shù),卻與‖x‖2是相容的.,設(shè)‖?‖是一種算子范數(shù), 則,,矩陣的范數(shù)與矩陣的特征值之間也有密切的聯(lián)系.,設(shè)?是矩陣A的特征值, x是對應的特征向量,則有,Ax= ?x,利用向量和矩陣范數(shù)的相容性, 則得,|?|‖x‖=‖?x‖=‖Ax‖?‖A‖‖x‖,于是

27、 |?|?‖A‖,設(shè)n階矩陣A的n個特征值為?1, ?2, …, ?n, 則稱,為矩陣A的譜半徑.,對矩陣的任何一種相容范數(shù)都有,? (A)?‖A‖,另外, ?? >0, ?一種相容范數(shù), 使 ‖A‖?? (A)+?,,任何兩種矩陣范數(shù)也具有等價性 m ‖A‖?? ‖A‖? ?M ‖A‖? ,? A?Rn?n,矩陣序列的收斂性定義為,§4 線

28、性方程組固有性態(tài)與誤差分析,§4.1 線性方程組的固有性態(tài),考慮線性方程組,其精確解為 x*=(-9800b1+9900b2 , 9900b1-10000b2)T,,若把線性方程組變?yōu)?解為x=(-9800b1+9900b2-19700? , 9900b1-10000b2+19900?)T,可見 x-x*=(-19700? , 19900?)T,解的誤差可能放大到右端項誤差的近2萬倍。

29、,若把線性方程組變?yōu)?因系數(shù)行列式為0，則線性方程組可能無解.,這種由于原始數(shù)據(jù)微小變化而導致解嚴重失真的線性方程組稱為病態(tài)方程組, 相應的系數(shù)矩陣稱為病態(tài)矩陣.,,設(shè)線性方程組 Ax=b,系數(shù)矩陣是精確的, 右端項有誤差?b, 此時記解為x+?x ,則 A(x+?x) =b+?b,于是

30、 Ax+A ?x =b+?b , A?x =?b, ?x =A-1?b,所以 ‖?x ‖= ‖A-1?b‖? ‖A-1‖ ‖?b‖,又由于 ‖b‖= ‖Ax‖? ‖A‖ ‖x‖,因此 ‖?x ‖‖b‖? ‖A‖‖A-1‖‖?b‖‖x‖,即,,可見誤差是否可控制取決于‖A‖ ‖A-1‖的大小。,再設(shè)b是精確的, A有誤差?A,

31、此時記解為x+?x ,則 (A+?A)(x+?x) =b,則有 Ax+A?x +?A(x+?x) =b, A?x +?A(x+?x) =0,所以 ?x = ?A-1?A(x+?x),于是 ‖?x‖? ‖A-1‖‖?A‖‖x+?x‖,也就是,記 Cond(

32、A)=‖A‖‖A-1‖, 稱為方程組Ax=b或矩陣A的條件數(shù)。如果 Cond(A)>>1，則稱A為病態(tài)矩陣。,,經(jīng)常使用的條件數(shù)有 Condp(A)=‖A‖p‖A-1‖p p=1，2，?。,當A為對稱矩陣時，有 Cond2(A)= ‖A‖2‖A-1‖2=|?1|/|?n|其中?1，?n分別是A的按絕對值最大和最小的特征值。,例如，對前面方

33、程組的系數(shù)矩陣A有,Cond1(A)= Cond?(A)= 39601， Cond2(A)?39206,由于計算條件數(shù)運算量較大, 實際計算中若遇到下述情況之一,方程組就有可能是病態(tài)的:,,(1) 矩陣元素間數(shù)量級差很大 ,且無一定規(guī)律;,(2) 矩陣的行列式值相對來說很小 ;,(3) 列主元消去法求解過程中出現(xiàn)量級很小的主元素;,(4) 數(shù)值求解過程中 , 計算解x的剩余向量r=b-Ax已經(jīng)很小, 但x仍不符合要求.,§4.

34、2 預條件和迭代改善,1. 線性方程組的預條件處理,對病態(tài)方程組Ax=b ,考慮線性方程組,其中,,,稱之為預條件方程組, 顯然與原方程組等價.,(1)條件數(shù),(2)方程組Cz=d容易求解。,比Cond(A)明顯小.,對于一般的矩陣A沒有十分有效的方法去選擇預條件矩陣. 當A是對稱正定矩陣時,可取C .,2. 線性方程組的迭代改善,設(shè)已求得方程組Ax=b的近似解x(1), 計算剩余向量,r(1)=b-Ax(1),再

35、求解余量方程組Ax=r(1), 得到解,則x(1)的迭代改善解為:,,,可逆矩陣C稱為預條件矩陣.矩陣C應滿足條件,練習題,第50頁 習題2,2-2(1), 2-6(1)（用二位浮點計算）,,練習題,第50頁 習題2,2-3(1), 2-4, 2-5, 2-8,,練習題,第50頁 習題2,2-10, 2-11, 2-12, 2-14, 2-16, 2-17,,課 堂 練 習,1. 對