线性方程组的误差分析

max
Step 3 For i = 1, …, n 必须等X(k)完全计算 n b ∑( 好了才能计算X(k+1)，i 因此aij X 0 j ) j ≠i Set X = 需要两组向量存储。 两组向量存储 j 需要两组向量存储。=1 i
ii
aii
; /* compute xk */
Step 4 If || X X 0 ||∞ = max | X i X 0 i | < TOL then Output (X[ ]); 1≤ i ≤ n 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 */
Jacobi Iterative Method
a11 x1 + a12 x 2 + ... + a1n x n = b1 a21 x1 + a 22 x 2 + ... + a 2 n xn = b2 aii ≠ 0 ... ... ... ... an1 x1 + a n 2 x 2 + ... + a nn x n = bn
… … …
…
( ( ( ( ( xnk +1) = 1 (an1 x1k +1) an 2 x2k +1) an3 x3k +1) L ann1 xnk +1) + bn ) 1 ann v v ( k +1) v ( k +1) v (k ) 1 1 写成矩阵形式 矩阵形式： 写成矩阵形式： x = D ( Lx + Ux ) + D b v (k +1) v (k ) v (D + L)x = U x + b v v ( k +1 ) 1 v ( k ) 1 x = ( D + L ) Ux + ( D + L ) b v Gauss-Seidel B f 迭代阵
0 . 106 × 10 v 3 || δ b || 2 v 4 v ≈ 0 .513 × 10 < 0 .01 % 此时精确wk.baidu.com为 x* = 此时精确解 精确解为 1.0203 || b || 2
δx = x * x = 2 .0203
v
v
v
2
v || δ x ||2 ≈ 2.0102 > 200% v || x ||2
常用条件数有： 常用条件数有： cond (A)1
max | λ | 特别地， 对称， 特别地，若 A 对称，则 cond ( A) 2 = min | λ |
cond (A)∞
cond (A)2
= λmax ( AT A) / λmin ( AT 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 。 可逆， 正交 正交， 可逆
§2 线性方程组的误差分析
/* Error Analysis for Linear system of Equations */
v v v v 有何影响？ 求解 A x = b 时，A 和 b 的误差对解 x 有何影响？ v v v v 精确， 设 A 精确，b 有误差 δ b ，得到的解为 x + δ x ，即
2.9 × 106
注：一般判断矩阵是否病态，并不计算A1，而由经验 一般判断矩阵是否病态，并不计算 得出。 得出。 行列式很大或很小（如某些行、列近似相关）； 行列式很大或很小（如某些行、列近似相关）； 元素间相差大数量级，且无规则； 元素间相差大数量级，且无规则； 主元消去过程中出现小主元； 主元消去过程中出现小主元； 特征值相差大数量级。 特征值相差大数量级。
v v v v A( x + δ x ) = b + δ b
v v 1 δ x=A δb
绝对误差放大因子
v v 1 || δ x || ≤ || A || || δ b ||
1 || A || v ≤ v || x || || b ||
v v v || b || = || Ax || ≤ || A || || x || 又
§3 Jacobi & Gauss-Seidel Iterative Methods
Gauss - Seidel Iterative Method ( ( ( ( ( x1k +1) = 1 ( a12 x2k ) a13 x3k ) a14 x4k ) L a1n xnk ) + b1 ) a11 ( ( ( ( ( 只存一组向量即可。 x2k +1) = 1 ( a21 x1k +1) a23 x3k ) 只存一组向量即可k ) + b2 ) a24 x4k ) L a2n xn。 a22 ( ( ( ( ( x3k +1) = 1 (a31 x1k +1) a32 x2k +1) a34 x4k ) L a3n xnk ) + b3 ) a33
相对误差放大因子
v v || δ x || || δ b || 1 v ≤ || A || || A || v || x || || b ||
§2
v v. Error Analysis for Ax = b
v v v || A || || A1 是关键 b 精确，A有误差 δ A||，得到的解为 x + δ x ，即 设 精确， 有误差 的误差放大因子， 的误差放大因子，称为v v v A的A + δ A，记为δ x ) = b , 条件数， x + 的条件数 )(记为cond (A) ( 越病态， 越大 则 A 越病态， 难得准确解。 v v v v v v v 难得准确解。 v ( A + δ A ) x + ( A + δ A )δ x = b A( x + δ x ) + δ A( x + δ x ) = b v v v v v 1 ( A + δ A )δ x = δ A x δx = A δA( x + δx ) v v 1 v A( I + A δ A )δ x = δ A x || δx || 1 v v ≤ || A || || δA || v v || x + δx || δ x = ( I + A 1δ A) 1 A 1δ Ax || δA || 1 = || A || || A || (只要δ A充分小，使得 只要 充分小， Wait a minute … || A || Who 1δ A that ( I+ || || δ A || < 1 ) || A said || ≤ || A 1 A 1 δ )
§3 Jacobi & Gauss-Seidel Iterative Methods
Algorithm: Jacobi Iterative Method
v v v Solve Ax = b given an initial approximation x ( 0.) 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. 迭代过程中， 迭代过程中，A 的元素 Output: approximate solution X[ ] or a message of failure. 不改变，故可以事先调整 事先调整好 不改变，故可以事先调整好 A 使得 A bit wasteful, Step 1 Set k = 1; aii What if aA不可逆。 不可逆。 ≠ 0，否则 isn’t it? ， 不可逆 = 0? Step 2 While ( k ≤ N ) do steps 3-6
§2
v v. Error Analysis for Ax = b
1 A= 例： 0 . 99
v 1 . 99 0 . 99 , b = 1 . 97 0 . 98
计算cond (A)2 。 计算 解：考察 A 的特征根 det( λ I A ) = 0
cond ( A)2 =
§2
1 1 2 例：Hilbert 阵 H n = 1 n
1 2 1 3
1 n +1
1 2 n 1
1 n
v v. Error Analysis for Ax = b
cond (H2)∞ = 27 cond (H6)∞ =
cond (H3)∞ ≈ 748 cond (Hn)∞→ ∞ as n → ∞
x = 1 ( a x ... a x + b ) 12 2 1n n 1 1 a11 1 ( a x ... a x + b ) x2 = 21 1 2n n 2 a 22 ... ... ... ... 1 x n = a ( a n1 x1 ... a nn 1 x n 1 + bn ) nn
§2
v v. Error Analysis for Ax = b
近似解的误差估计及改善： 近似解的误差估计及改善： v v v v v v v 设 Ax = b 的近似解为 x *，则一般有 r = b A x * ≠ 0 v v v || x x* || cond (A) || r || v ≤ v HW: p.66 || x || b #2, #4, #5 改善方法： 改善方法： v v v x1 ; Step 1: Ax = b 近似解 v v v v 误差上限 r1 = b A x 1 ; Step 2: 可被精确解出， 若 d 1 可被精确解出，则有 v v v v v v v v 1 1 Step 3: A d 1 = r1 d 1 ; x2 = x1 + A (b Ax1 ) = A b v v v v x2 就是精确解了。 Step 4: x 2 = x 1 + d 1 ; 就是精确解了。
§2
v v. Error Analysis for Ax = b
注:
cond (A) 的具体大小与 || || 的取法有关，但相对 的取法有关， 大小一致。 大小一致。 cond (A) 取决于 ，与解题方法无关。 取决于A，与解题方法无关。
v v || δ A || || δ b || || δ x || cond ( A ) v v ≤ || A || + || b || || x || 1 cond ( A ) || δ A || || A ||
ε 经验表明： 不是非常病态（例如： 经验表明：若 A 不是非常病态（例如： cond ( A)∞ < 1）， 则如此迭代可达到机器精度； 病态， 则如此迭代可达到机器精度；但若 A 病态，则此算法也 不能改进。 不能改进。
§3 Jacobi 法和 Gauss - Seidel 法
/* Jacobi & Gauss-Seidel Iterative Methods */
写成矩阵形式： 写成矩阵形式： 矩阵形式
A=
D
U
v v Ax = b
v v (D + L + U )x = b v v v Dx = ( L + U ) x + b v v v 1 1 x = D (L + U )x + D b
L
Jacobi 迭代阵
B
v f
v v ( k +1) v (k ) 1 1 x = D (L + U )x + D b
1 精确解为 v 精确解为 x = . 1 9900 9800 1 = A 9900 10000
λ1 ≈ 39206 >> 1 λ2
λ1 = 1 .980050504 λ 2 = 0 .000050504
测试病态程度： 测试病态程度： v v 0 . 97 × 10 4 给 b 一个扰动 δ b = ，其相对误差为 3
|| A || || A 1 || v || δ x || || A 1 || || δ A || || A || v ≤ = || x || 1 || A 1 || || δ A || 1 || A || || A 1 || || δ A || || A ||
is invertible? || || δ A