擾動(dòng)誤差分析-精選.

1、第2章線性代數(shù)方程組數(shù)值解法I :直接法 考慮方程組 Ax b ( A Rn n非奇異，b Rn且b 0 ) (261) 設(shè)A有誤差 A,b有誤差b，則因此引起解x有誤差，即有擾動(dòng)方程組 (A A)(x x) b b 現(xiàn)在來(lái)研究如何通過(guò) A和b對(duì)x的影響作出估計(jì)。 定理2.6.1 設(shè) 方程組(2.6.1)中A,b分別有擾動(dòng)A，b，因而解向量有誤差x; 又A足夠小，使得 A1 A 1，則有誤差估計(jì)式 A 11A 1 Ia1! II A (A (A 證明 (A)x ( A)( x) 1 b A 1( A)x A 1( A)( x) 兩邊取范數(shù)有 INI |Alb IIAl 州 (1 A1 A x

2、A1 I AMI A1 1 ( 1 All Air III A) x| 得到| x b Ax) 又注意到 Ax b 有 A b從而得到_x b，故上述不等式左邊乘以一x，右 -i1，第二項(xiàng)乘以，并從括號(hào)中提出| A，則得(263) x b 定理的結(jié)果實(shí)際包含兩種特殊情形： 邊圓括號(hào)第一項(xiàng)乘以 (1) A精確，即 A 0，b有擾動(dòng)b，從而A(x x) b b A有擾動(dòng)A, b精確，即b 0,這時(shí)(A A)(x x) b x k|A A x 1 A1|A A 當(dāng)IA 1| A很小時(shí)，上式可近似表示為 x x 2 條件數(shù)與病態(tài)方程組 定義2.6.1設(shè)A為非奇異矩陣，稱數(shù) A1 A即 cond(A)A

3、A| 為矩陣A的條件數(shù)。 矩陣A的條件數(shù)的一些基本性質(zhì)： (1) 任何非奇異矩陣A，對(duì)任一算子范數(shù)| |均有 cond(A) 1 cond (A) cond( A ) cond ( A) cond(A) 根據(jù)定義 A 2. max(ATA),可得 cond (A)2 A12A max (A A) min (AT A) 若U為正交陣，即uu T I，則 con d(U)21 又A非奇異，則cond (A)2 cond(AU )2 cond (UA)2 (4)設(shè)1與n為A按模最大和最小的特征值，則 con d(A) n 特別地，若A At (即A對(duì)稱)，則cond(A) 若A對(duì)稱正定，則 cond

4、(A)2 證明略 定理 262 (事后誤差估計(jì)) 設(shè)方程組 Ax b, A非奇異，b 0 , x是精確解, x是近似解, 剩余向量 r b Ax，則有估計(jì)式 1 r cond(A) | b cond(A) r ibi 證明: b，得 r b Ax Ax Ax A(x x)，從而 x x A 1r，于是 八1 x x A r 因Ax 1 ，又由x :，于是得估計(jì)式右端 rA cond(A) b b 類似地，由上述r A(x x)，得IHI |a| ，或A ，由x A 1b得 A1 b，綜合兩式得估計(jì)式兩端 1 _ IHI |r|_ cond(A) b| IaIIa1 b|x 例 2.6.1 用平

5、方根法解線性方程組 _ 10787 一 r廠 _ 32 7565 23 86109 33 -75910 - -31 粗看來(lái)沒(méi)什么特別，可以驗(yàn)證其特確解為 現(xiàn)虛假設(shè)方程組右端項(xiàng)可能有諼總 比如設(shè)b有撫動(dòng)Bb,即 (32+0,1? 23-0. 1, 33+0. 1, 31-0. 1)T 再解方程組，劇有擾動(dòng)解 = (黑 2廠口 6,1 5t-Ll)r 可見b的微小擾動(dòng)bb對(duì)解X有相當(dāng)大的影響6x.為此，檢查方程組的條件 數(shù).采用2-條件數(shù)由A算出牯征值 入1 = 30/887入嚴(yán)3”58入嚴(yán)厲8431 入嚴(yán)Q. 01015 注意方程組是一個(gè)對(duì)稱組*可得方程組條件數(shù) A 30.2887 叭尸廠麗嚴(yán)亦

6、 a 可見方釋組相當(dāng)病態(tài)+上述出St的規(guī)象實(shí)厲有因，解方程組時(shí)注意檢瓷 條件數(shù)，碉有助于避免盲目地計(jì)算遠(yuǎn)里.由5b =( 0, 1 ,-Q. 1,0, 1廠0, 1) 5x=( 8. 2 廠13. 6r3. 5廠2 1)T 可算出 I b II 5=60, 025 于是強(qiáng)實(shí)際相對(duì)誤差 II 5 b II =0* 21x1 =2 I 5 x II ,16, 397 世=8,198 II x II n 2 II 5xtt3=16, 397 而按定理牯計(jì) Lllb 2984 sJ2LL. I Xll2 務(wù)切l(wèi)bl2 0L 2 9,943 由此看出,理論分析與實(shí)際計(jì)算基本一致. 3 事后誤差估計(jì) 定理

7、2.6.2設(shè) 方程組Ax b , A非奇異，A非奇異，b 0, x是精確解，x 是近似解，剩余向量 r b Ax,則 有估計(jì)式 例題講解2 例題 2.1 對(duì)方程組 Ax=b, A 非奇異不一定能作順序 Gauss 消去過(guò)程， 或者 說(shuō)， A 非奇異不一定有 LU 分解。 證 這類命題只需舉一個(gè)反例即可。 反例要盡可能簡(jiǎn)單， 這里可考慮 2 階、 元素為 1 或 0。構(gòu)造反例一般要經(jīng)過(guò)多次“失敗 - 修正”過(guò)程。 10 考慮 A= ，易見 A 非奇異。顯然，順序 Gauss 消去過(guò)程的 01 第 1 步就不能進(jìn)行。從 LU 分解來(lái)看，設(shè)有 A=LU 分解，則有 1 0 = 1 b c = b c

8、 0 1 a 1 d ab ac d 于是有 b=0 和 ab=1 同時(shí)成立，這就自相矛盾了。 例題 2.2 設(shè)方程組 0.001 2.000 3.000 x1 1.000 1.000 3.712 4.623 x2 = 2.000 2.000 1.072 5.643 x3 3.000 試手算(或輔以計(jì)算器)分別用 (1) 順序Guass消去法 (2) 列主元Gauss消去法 ( 3)直接三角分解法 求解，要求計(jì)算中取 4位有效數(shù)字， 最后結(jié)果舍入成 3位有效數(shù)字。 上述方程組 用 4 位浮點(diǎn)數(shù)進(jìn)行計(jì)算而舍入為 3 位有效數(shù)字的精確解為： *T x* ( 0.490 0.0510 0.368)T

9、 解 (1)順序Gauss消去過(guò)程計(jì)算有 1.000 l21 1000 0.001 2.000 3.000 A b 1.000 3.712 4.623 1.000 第 1 步消元 21 l31 2000 2.000 1.072 5.643 1.000 0.001 2.000 3.000 1.000 2004 3005 1002 第 2 步消元l32 1.977 4001 6006 2006 0.001 2.000 3.000 1.000 2004 3005 1002 5015 2006 回代過(guò)程求解并舍入得 x 0.399,0.0998, 0.400, T 3)令 A=LU ，用 0.001

10、2.000 3.000 1000 2004 3005 2000 1.997 1 5015 由 Ly=b 解得 1.001, 1.002, 2006 T 有 Ux=y 解得 0.400, 0.0998, 0.400 T (2)列主兀Gauss消去過(guò)程計(jì)算有 0.001 2.000 3.000 1.000 Ab 1.000 3.712 4.623 2.000 選主元 a31 2.000換行 E1E 2.000 1.072 5.643 3.000 2.000 1.072 5.643 3.000 第1步消元 l21 0.5000 1.000 3.712 4.623 2.000 l31 0.0005 0

11、.001 2.000 3.000 1.000 2.000 1.072 5.643 3.000 3.716 1.802 0.5000 選主元,仍為 a22 2.001 3.003 1.002 2.000 1.072 5.643 3.000 3.716 1.802 0.5000 第 2 步消元 l32 0.6300 1.868 0.6870 回代過(guò)程求解并舍入 0.0513, 0.368 T x 0.490, Doolittle 分解計(jì)算得 1 與精確解 x* 比較可見，列主元 Gauss 消去法的解精確度最高；其他兩種方法的 解精確度較差，但彼此接近(這正好符合兩種方法實(shí)質(zhì)是一樣的情況) 。 例

12、題2.3 Gauss消去法的一種自然而又簡(jiǎn)單的改進(jìn)是所謂Gauss-Jordan消 (k)(k) 去法。先考慮順序Gauss消去法過(guò)程的情形。它只是把 a這一列中a下面的 kkkk (k)(k) 元素消為0,而Gauss-Jordan消去過(guò)程則把a(bǔ)這一列元素的a以外全部消去為 kkkk (k) 0,并且約化a kk=1。為此，需進(jìn)行n步消元，第n列也消為只剩下一個(gè)元素1, 其余均為0 (因此， n工0在這里也是必要的)。這樣一來(lái)，不用回代過(guò)程, 方程組的解就在b的位置上?，F(xiàn)在依據(jù)上述導(dǎo)引，做 (1) 試用Gauss-Jordan消肖去法解方程組 111 x. 2x2 = 0 211X3 (2)

13、 試給出Gauss-Jordan算法的核心部分。 (3) 由上述兩點(diǎn)能得出什么思考？ 解(1)用增廣矩陣的演變來(lái)描述求解過(guò)程 11 111111 1 2 2 00 11 1 21110 313 100210 02 01110102 00260013 于是有x= (2, 2, 3) (2) Gauss-Jordan消肖去法算法的核心部分為： 對(duì)k=1,2,，n,做 a akj J -(j=k，k+1,，n，n+1) akk 對(duì)k=1, 2,，nA i工k,做 aij J a j-ajkakj(j=k+1,n,n+1) 輸出解 Xi =ai,n 1 (i=1,2,n) (4) 從上述做法可引發(fā)如下

14、思考：Gauss-Jordan消去過(guò)程自然也可考慮 加上選主元和換行的技巧；就計(jì)算量而言，Gauss-Jordan消去過(guò)程 3 顯然要大(可推導(dǎo)其乘除運(yùn)算次數(shù)為級(jí)，而順序Gauss消去過(guò) 2 3 程乘除運(yùn)算次數(shù)為 級(jí))。因此，用Gauss-Jordan消去法解方程組 3 并不經(jīng)濟(jì)。但它有什么用途呢？這又引出一個(gè)新的思考。 例題2.4設(shè)L Rn*n為非奇異下三角陣，試 (1)寫出求解方程組Lx= f的計(jì)算公式; (2) 統(tǒng)計(jì)上述求解過(guò)程的乘除法次數(shù)； (3) 給出求L 1的計(jì)算公式。 解(1)這是解下三角方程組的遞推公式: x1 f1 /l11 i1 x2 (filij xj)/lii (i 2

15、,3, ,n) j1 (2) 乘除法運(yùn)算，第一式有1次，第二次i=2時(shí)有2次，i = n時(shí)有n次, 故上述公式的乘除運(yùn)算次數(shù)有 1 + 2+ n = n(n 1)/2 (次) (3) 因L非奇異，L1存在，且由矩陣知識(shí)知L1也為下三角陣，記L1 (Vj), 由 LL 1 I ，即 l11 V11 1 LL 1 l21 l22 V21 V22 1 ln1 l n2 l nn Vn1 Vn2 Vnn 1 考慮 I 的 對(duì) 角 線 元 素 ， 顯 然 有 li iVii 1(i 1,2,n), 于 1/lii (i 1,2,n)。 考慮 I 對(duì)角線以下的元素，即對(duì) i 2,3, ,n; j 1,2,

16、 ,i 1 ，則有 n lik Vkj k1 i lik Vkj k1 lik Vkj lijVij 0 于是得 Vij(lik Vkj ) /l ii (i 2,3,n; j 1,2, ,i 1) 是得 Vii 例題2.5設(shè)A= (aj)Rn*n對(duì)稱，其順序主子式厶i工0(i=1,2，,n),試 (1) 證明分解A=LDL存在唯一，其中L為單位下三角陣，D為對(duì)角陣； (2) 寫出利用此分解求解方程組 Ax=b的步驟(這稱為改進(jìn)的平方根法); ( 3) 用改進(jìn)的平方跟法解方程組 3 3 5 x1 10 3 5 9 x2 16 5 9 17 x3 30 解 （ 1）由 A 的順序主子式不為零，故

17、存在唯一分解式 A=LU ， L 為單位下三角 陣，U為上三角陣。改寫 A=LU=LDU 0 , D為對(duì)角陣，U0為單位上三角陣。由 A對(duì)稱，有A=AT= （LDU0）T=U； （DLt），又由分解的唯一性即得Uj=L或 U0=Lt，于是代入可得 A=LU=LDU o=LDLt存在唯一。 （2）將 A=LDL T 代入 Ax=b 得 LDL Tx=b，令 DLTx=y （即 LTx=D-1y）， 則Ly=b，改進(jìn)平方根法如下： 1）將A直接分解為A=LDL T，即求出L，D ； 2）求解下三角方程組Ly=b，得y; 3）求解上三角方程組LTx=D-1y，得x。 （ 3）令 A=LDL T 3

18、3 51 d1 3 5 9 l21 1 d2 5 9 17 l31 l32 1 d3 由矩陣乘法并對(duì)比等式兩邊得 d13 d2 2 d3=2/3 l 21 l32 2 l315/3 1 y1 10 解 Ly=b 即 1 1 y2 16 5/3 21 y3 30 1 1 5/3 x1 1/3 解 LTx=D-1y 即 1 2 x2 1x3 1 l 21 l31 1 l32 1 10 得 y= 6 4/3 10 1/ 216 得 x=（ 1 ， -1 ， 2） T 3/2 4/3 可見改進(jìn)平方根法比原始平方根法避免了開方運(yùn)算; 另一方面， 改進(jìn)平方根法對(duì) 滿足 LU 分解條件的對(duì)稱矩陣（不一定要正

19、定）都適用。 例題 2.6 已知方程組 Ax=b 的系數(shù)矩陣形如： 其中 A 的 n-1 階、主子距陣， *為其實(shí)數(shù)，。試設(shè)計(jì)一個(gè)求解上述方程組基于追趕 的解決方案 解 把 A 記為 A= Bn 1 T u 其中 Bn-1 為 A 的 n-1 階（對(duì)三角型）主子距陣， V, U R.對(duì)Bn 1用追趕法可 作三角分解 11 2 d1 1 2 2 2 1 d2 2 Bn 1= p2 1 =P-Q n2 n2n2 1 n2 n 1 n1 dn 2 再對(duì) A 作三角分解， A= Bn 1 P 0 Q Z PQ P WT 1 0 dn = WT Q WZ T d 其中Z RN-1）由下三角方程組 Pz=

20、v求出，W R（N-1）由方程組WT Q=uT（即下三角 方程組QtW=U）求出，dn由WT+dn= n算出。由此，求解上述方程組的解決方 案如下： （1）對(duì) n 1用追趕法求出 P， Q； （2）解下三角方程組Pz=v,求出z; （ 3）解下三角方程組 QtW=U 求出 W； （ 4）計(jì)算 dn = n-W T z; P0 5）前推解下三角方程組 WPT 10 y=b 求出 y; 6）回代解上三角方程組 Z dn x=y 求出 x; aibi xi c2 a2 b2 x2 例題 2.7 設(shè)三對(duì)角方程組 Ax=b, 其中 A=c3 a3b3 x3 x= cn i an i bn i xn i

21、cnan xn d1 d2 d= d3 dn 1 dn 試建立A的順序主子式k(k=1,2,n)的逆推公式 (2) 試設(shè)計(jì)求解上述三對(duì)角方程組的另一種有別于一般追趕法 新算法。 (基于 LU 分解 )的 解 補(bǔ)充記 0 =1 ，則顯然有 1=a1 2 =a1a2 b1c2 =a2 1 b1c2 般地，猜想 k = ak k i bk iCk k2(k=3,4,n) 用歸納法證明：當(dāng) k=2 時(shí)，上式成立。現(xiàn)假設(shè)有 則由 k i=ak i k 2 bk 2ck i k 3 aibi c2 a2 b2 k= =ak k ibk ickk 2 ck 2ak 2bk 2 ck i ak ibk i c

22、kak 可知上述猜想成立。綜上所述，可得遞推公式： 1a1 k ak k 1 bk 1ck k 2 (k 23., n) (2)由Ax=d的第一個(gè)方程可得 b1d1 X1X2 即得x1由x2表示的關(guān)系式。將它代入第二個(gè)方程c2x1 a1x2 b2x3 d2，又可得 X2由X3表示的關(guān)系式如此繼續(xù)代入，直到第n 1個(gè)方程，于是有如下形式 的關(guān)系式 XkfkXk1 gk(k=2,3,n 1) 其中fk和gk待定。現(xiàn)將形式如(2)的Xk 1關(guān)系式代入Ax=d的第 k(k=2,3, n 1) 個(gè)方程CkXk 1 akXkbkXk 1 dk，則有 Ck f k 1 Xk Ckgk 1akXk bkXk

23、1 dk 整理可得 Xk bk Ck f k 1 ak Xk 1 d k Ck g k 1 Ck fk 1 ak (k=2,3,n 1) (3) 與關(guān)系式(2)比較得 bk k Ck ak gk dk ck g k 1 ck fk 1 ak (k=2,3,n 1) 注意到方程組中的記號(hào), 0 , bm 0。由中取k=1,可得f1 色4 a1 a1 代入(1)得 X1仏2 g1，可見(2)式對(duì)k=1也成立；又由中取k=n,則可得fn 0 , gn蟲Cn1，而由于方程組Ax=d的第n個(gè)方程CnXn1 anXn dn可以寫成 an fn 1 an CnXn 1 anXn bn Xn 1 dn，其中g(shù) 0， Xn 1 0，從而上述的代入可直至第n 個(gè)方程，也就是說(shuō)(2)和(3)式對(duì)k n也成立，于是有Xn gn。綜上所述，可得解 三對(duì)角方程組的新算法如下： b1d