§2线性方程组的误差分析 / Error Analysis for Linear system of Equations * 求解Ax=b时,A和b的误差对解x有何影响? >设A精确,误差δb,得到的解为x+δx,即 绝对误差放大因子 A(x+8x)=b+8b →6x=A8b→6x<业4m 相对误差放大因子 又‖b=‖A‖s‖A‖·‖‖→ b‖ bxl‖l lSb b
§2 线性方程组的误差分析 /* Error Analysis for Linear system of Equations */ 求解 Ax b 时，A 和 的误差对解 有何影响？ = b x ➢ 设 A 精确， b 有误差 ，得到的解为 ，即 b x x + A x x b b ( + ) = + x A b −1 = || || || || || || 1 x A b − 绝对误差放大因子 || b || || Ax || || A|| || x || 又 = || || || || || || 1 b A x || || || || || || || || || || || || 1 b b A A x x − 相对误差放大因子
8 2 Error Analysis for Ax=b 设b精确 A|·A‖是关键x+Sx,即 的误差放大因子,称为 的条件数,记为cond(A4), 越大则A越病 难得准确解。 A(x+)+S4(x (A+d4)x+(4+d4)G=b =-4-S4(+&c →(4+d4)a=-4x l ac A4(I+A4)=-64x s‖A‖·G4‖ lx+&x‖ &c=-(+A SA) A SAx 4|-1A31D 4‖ A‖l (要高充分使得 sd#+4<1 is Invertib ‖4t·才田 ‖<‖A-‖‖4‖ A c‖-1-‖A-l·‖4‖ 4| A‖l
§2 Error Analysis for Ax b . = ➢ 设 b 精确，A有误差 ，得到的解为 ，即 A x x + A A x x b ( + )( + ) = A x x A x x b ( + )+ ( + ) = ( ) 1 x A A x x = − + − || || || || || || || || || || || || || || || || 1 1 A A A A A A x x x = + − − A A x A A x b ( + ) + ( + ) = A A x Ax ( + ) = − A I A A x Ax + = − − ( ) 1 x I A A A Ax 1 1 1 ( ) − − − = − + Wait a minute … Who said that ( I + A−1 A ) is invertible? (只要 A充分小，使得 || || || || || || 1 ) 1 1 − − A A A A || || || || 1 || || || || || || || || || || || || 1 || || || || || || || || || || || || 1 1 1 1 A A A A A A A A A A A A x x − = − − − − − 是关键 的误差放大因子，称为 A的条件数，记为cond (A) , 越 则 A 越病态， 难得准确解。 || || || || −1 A A 大
82 Error Analysis for Ax=b 注:cond(4)的昊体大小与‖·‖的取法有关,但相对 大小—致。 cond(4)取决于A,与解题方法无关 cond(a) dA‖l.‖b‖l ‖-1-cond(A川&A‖∥A‖A‖‖b‖ 纔常用条件数有: om(n)eom(42com(4)2=3(44-(4小 特别地,若A对称,则cond(42 max m 条件数的性质: 国A可逆,则cond(A)≥1; 国A可逆,a∈R则cond(aA)=cond(4 国A正交,则con(4)2=1; 国A可逆,R正交,则cond(RA)2=cond(4R)2=cond(4)2
§2 Error Analysis for Ax b . = 注: cond (A) 的具体大小与 || · || 的取法有关，但相对 大小一致。 cond (A) 取决于A，与解题方法无关。 + − || || || || || || || || 1 ( )|| || || || ( ) || || || || b b A A cond A A A cond A x x 常用条件数有： cond (A)1 cond (A) cond (A)2 ( )/ ( ) max A A min A A T T = 特别地，若 A 对称，则 min | | max | | ( ) 2 cond A = 条件数的性质： A可逆，则 cond (A)p 1； A可逆， R 则 cond ( A) = cond (A) ； A正交，则 cond (A)2=1； A可逆，R正交，则cond (RA)2 = cond (AR)2 = cond (A)2
82 Error Analysis for Ax= b 099 199 例:A= b 0990.98 1.97 精确解为x 98009900 计算cond(4)2 9900 10000 解:考察A的特征根 de(-A)=0→孔 1.980050504 2=-0.000050504 cond (a) 39206>>1 12 测试病态程度: 给b一个扰动Sb 0.97×104 ,其相对误差为 0.106×10 lb‖ ≈0.513×10-200% 2.0203
§2 Error Analysis for Ax b . = 精确解为 . 1 1 x = 例： = = 1.97 1.99 , 0.99 0.98 1 0.99 A b 计算cond (A)2 。 − − 9900 10000 9800 9900 A−1 = 解：考察 A 的特征根 det(I − A) = 0 0.000050504 1.980050504 2 1 = − = = 2 1 2 ( ) cond A 39206 >> 1 测试病态程度： 给 b 一个扰动 − = − − 3 4 0.106 10 0.97 10 b ，其相对误差为 0.513 10 0.01% || || || || 4 2 2 − b b 此时精确解为 − = 1.0203 3 x* − = − = 2.0203 2 x x * x 2 2 || || || || x x 2.0102 > 200%
82 Error Analysis for Ax= b 112 例: Hilbert阵H +1 2n cond(H2)oo= 27 cond(H3)≈748 eom(H)2=2.9×105comd(Hn)→∞asm→∞ 注:一般判断矩阵是否病态,并不计算A1,而由经验 得出。 G行列式很大或很小(如某些行、列近似相关) 元素间相差大数量级,且无规则; G主元消去过程中出现小主元; G特征值相差大数量级
§2 Error Analysis for Ax b . = 例：Hilbert 阵 = + 2 −1 1 1 1 1 3 1 2 1 1 2 1 1 n n n n Hn cond (H2 ) = 27 cond (H3 ) 748 cond (H6 ) = 2.9 106 cond (Hn )→ as n → 注：一般判断矩阵是否病态，并不计算A−1，而由经验 得出。 行列式很大或很小（如某些行、列近似相关）； 元素间相差大数量级，且无规则； 主元消去过程中出现小主元； 特征值相差大数量级
8 2 Error Analysis for Ax=b 近似解的误差估计及改善: 设Ax=b的近似解为x则一般有P=b-Ax*≠0 x-x <(cond (A) HW: p 66 狄改善方法: #2,#4,#5 S1:A=b→近似解x1; Step 2: =b-A, 若d可被精确解出,则有 Sep3:Ad1=→l1, x2=x, +A(6-Ax=ab Step 4: x2=x,+d,i 2就是精确解了。 经验表明:若A不是非常病态(例如:cond(A)<1), 则如此迭代可达到机器精度;但若A病态,则此算法也不 能改进
§2 Error Analysis for Ax b . = ➢ 近似解的误差估计及改善： 设 Ax b 的近似解为 ，则一般有 = x* * 0 r = b − Ax b r x x x || || || || || * || − cond (A) 误差上限 改善方法： Step 1: Ax = b 近似解 ; x1 Step 2: ; 1 b Ax1 r = − Step 3: ; 1 1 d1 Ad r = Step 4: ; x2 x1 d1 = + 若 可被精确解出，则有 就是精确解了。 d1 x x A b Ax A b 1 1 1 2 1 ( ) − − = + − = 2 x 经验表明：若 A 不是非常病态（例如： ）， 则如此迭代可达到机器精度；但若 A 病态，则此算法也不 能改进。 cond (A) 1 HW: p.66 #2, #4, #5
§3 Jacobi法和 Gauss- Seidel法 k Jacobi Gauss-Seidelliterative Methods >Jacobi lterative method -(-a1 1 x+b aux,+au2x2 +.+anxn=b a,,x、n…飞M七=b,n≠0 x,= (-a b2) 1 a 十…十a…x ( .x,+ 写成矩阵形式: Ax=b冷(D+L+U)x=b U 分D=-(L+U)x+b →x=-D(L+U)x+Db B Jacobi迭代阵 x (k+1) -D(L+U)x(k) +D b
§3 Jacobi 法和 Gauss - Seidel 法 /* Jacobi & Gauss-Seidel Iterative Methods */ ➢ Jacobi Iterative Method + + + = + + + = + + + = n n nn n n n n n n a x a x a x b a x a x a x b a x a x a x b ... ... ... ... ... ... ... 1 1 2 2 21 1 22 2 2 2 11 1 12 2 1 1 ( ) ( ) ( ) = − − − + = − − − + = − − − + n nn− n− n nn n n n n n a x a x b a x a x a x b a x a x a x b a x 1 1 1 1 2 1 1 2 2 2 2 2 1 2 2 1 1 1 1 1 ... 1 ... ... ... ... ... 1 ... 1 0 ii a 写成矩阵形式： A = L U D Dx L U x b Ax b D L U x b = − + + = + + = ( ) ( ) x D L U x D b 1 1 ( ) − − = − + + B f Jacobi 迭代阵 x D L U x D b k k ( 1) 1 ( ) 1 ( ) + − − = − + +
83 Jacobi Gauss-Seidel Iterative Methods Algorithm: Jacobi iterative Method Solve ax=b given an initial approximation x(o Input: the number of equations and unknowns n; the matrix entries a[[; the entries bl l; the initial approximation nn l. tolerance TOL maximum number ofiteratiop 迭代过程中,A的元麦 Output: approximate solution Step 1 set k=1; Ⅺ(不改变,故 a bit wasteful, Step 2 whi steps 3-6 a≠0 sn’tit? 必须等w完全计算 好了才能计算x,因此aG 需要两组向量存储。 Impute xk * Step4If‖X-X0|l X-XO, TOL then Output (XI D; STOP: /* successful Step 5 For i=1,.,n Set X0[sXI /update X0*/ Step 6 Set k ++ Step 7 Output(Maximum number of iterations exceeded); STOP. / unsuccessful *
§3 Jacobi & Gauss-Seidel Iterative Methods Algorithm: Jacobi Iterative Method Solve given an initial approximation . Input: the number of equations and unknowns n; the matrix entries a[ ][ ]; the entries b[ ]; the initial approximation X0[ ]; tolerance TOL; maximum number of iterations Nmax. Output: approximate solution X[ ] or a message of failure. Step 1 Set k = 1; Step 2 While ( k Nmax) do steps 3-6 Step 3 For i = 1, …, n Set ; /* compute xk */ Step 4 If then Output (X[ ]); STOP; /* successful */ Step 5 For i = 1, …, n Set X 0[ ] = X [ ]; /* update X0 */ Step 6 Set k ++; Step 7 Output (Maximum number of iterations exceeded); STOP. /* unsuccessful */ Ax b = (0) x ii n j j i i ij j i a b a X X = − = 1 ( 0 ) X X Xi X i TOL i n − = − || 0 || max | 0 | 1 What if aii = 0? 迭代过程中，A的元素 不改变，故可以事先调整好A 使得 aii 0，否则 A不可逆。 必须等X(k)完全计算 好了才能计算X(k+1)，因此 需要两组向量存储。 A bit wasteful, isn’t it?
83 Jacobi Gauss-Seidel Iterative Methods Gauss-Seidel iterative method (k+1) (k) 122 134 aun,xm+bu) x41=1(a1x+-a<只存一组向量即可 (k+1) (k+1) (k+1) 24X 4~4 3nx+b3) (k+1) .,x (k+1) (k+1) n1~1 n2~2 anx+1-…-m1x1+bn) 写成矩阵形式:x(k+"=-D'(L+)+U)+D 冷(D+L)x(+=-L(k)+b (k+1 (D+L)UxA)+(D+L)b Gauss-Seidel B 迭代阵
§3 Jacobi & Gauss-Seidel Iterative Methods ➢ Gauss - Seidel Iterative Method ( ) 1 1 ( ) 1 ( ) 1 4 4 ( ) 1 3 3 ( ) 1 2 2 1 1 ( 1) 1 a x a x a x a x b a x k n n k k k k = − − − − − + + ( ) 1 2 ( ) 2 ( ) 2 4 4 ( ) 2 3 3 ( 1) 2 1 1 2 2 ( 1) 2 a x a x a x a x b a x k n n k k k k = − − − − − + + + ( ) 1 3 ( ) 3 ( ) 3 4 4 ( 1) 3 2 2 ( 1) 3 1 1 3 3 ( 1) 3 a x a x a x a x b a x k n n k k k k = − − − − − + + + + ( ) 1 ( 1) 1 1 ( 1) 3 3 ( 1) 2 2 ( 1) 1 1 ( 1) n k nn n k n k n k n nn k n a x a x a x a x b a x = − − − − − + + − − + + + + … … … … 只存一组向量即可。 写成矩阵形式： x D Lx Ux D b k k k ( 1) 1 ( 1) ( ) 1 ( ) + − + − = − + + D L x Ux b k k + = − + ( +1) ( ) ( ) x D L Ux D L b k k ( 1) 1 ( ) 1 ( ) ( ) + − − = − + + + B f Gauss-Seidel 迭代阵
83 Jacobi Gauss-Seidel Iterative Methods A mathematician about his colleague: He made a lot of mistakes but he made them in a good direction. I tried to copy this, but I found U out that it is very difficult to make good mistakes. 注:二种方法都存在收敛性问题。 有例子表明: Gauss-Seidel法收敛时, Jacobi法可能 不收敛)而 Jacobi法收敛时, Gauss-Seidel法也可能 不收敛。 p76#2给出了例子。 收敛性分析将在下节课讨论
A mathematician about his colleague: " He made a lot of mistakes, but he made them in a good direction. I tried to copy this, but I found out that it is very difficult to make good mistakes. " §3 Jacobi & Gauss-Seidel Iterative Methods 注：二种方法都存在收敛性问题。 有例子表明：Gauss-Seidel法收敛时，Jacobi法可能 不收敛；而Jacobi法收敛时， Gauss-Seidel法也可能 不收敛。 p.76 #2 给出了例子。 收敛性分析将在下节课讨论