1、2 线性方程组的误差分析线性方程组的误差分析 /*Error Analysis for Linear system of Equations*/求解求解 时，时，A 和和 的误差对解的误差对解 有何影响？有何影响？bxA bx 设设 A 精确，精确，有误差有误差 ，得到的解为，得到的解为 ，即，即bb xx bbxxA )(bAx 1|1bAx 绝对误差放大因子绝对误差放大因子|xAxAb 又又|1bAx|1bbAAxx 相对误差放大因子相对误差放大因子2 Error Analysis for .bxA 设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即bA xx bxxAA
2、 )(bxxAxxA )()()(1xxAAx|11AAAAAAxxx bxAAxAA )()(xAxAA )(xAxAAIA )(1xAAAAIx 111)(Wait a minute Who said that(I+A 1 A)is invertible?(只要只要 A充分小，使得充分小，使得)1|11 AAAA|1|1|1111AAAAAAAAAAAAxx 是关键是关键的误差放大因子，称为的误差放大因子，称为A的的条件数条件数，记为，记为cond(A),越越 则则 A 越病态，越病态，难得准确解。难得准确解。|1 AA大大2 Error Analysis for .bxA 注注:cond
3、(A)的具体大小与的具体大小与|的取法有关，但相的取法有关，但相对大小一致。对大小一致。cond(A)取决于取决于A，与解题方法无关。，与解题方法无关。|)(1)(|bbAAAAAcondAcondxx 常用条件数有：常用条件数有：cond(A)1cond(A)cond(A)2)(/)(minmaxAAAATT 特别地，若特别地，若 A 对称，则对称，则|min|max)(2 Acond条件数的性质：条件数的性质：A可逆，则可逆，则 cond(A)p 1；A可逆，可逆，R 则则 cond(A)=cond(A)；A正交，则正交，则 cond(A)2=1；A可逆，可逆，R正交，则正交，则 cond
4、(RA)2=cond(AR)2=cond(A)2。2 Error Analysis for .bxA 精确解精确解为为.11 x例：例：97.199.1,98.099.099.01bA计算计算cond(A)2。10000990099009800A 1=解：解：考察考察 A 的特征根的特征根 0)det(AI 000050504.0980050504.121 212)(Acond 39206 1 测试病态程度：测试病态程度：给一个扰动给一个扰动b 3410106.01097.0b，其相对误差为，其相对误差为%01.010513.0|422 bb 此时此时精确解精确解为为 0203.13*x 02
5、03.22*xxx 22|xx 2.0102 200%2 Error Analysis for .bxA 例：例：Hilbert 阵阵 12111131211211nnnnnHcond(H2)=27cond(H3)748cond(H6)=2.9 106cond(Hn)as n 注：注：一般判断矩阵是否病态，并不计算一般判断矩阵是否病态，并不计算A 1，而由经验，而由经验得出。得出。行列式很大或很小（如某些行、列近似相关）；行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；元素间相差大数量级，且无规则；主元消去过程中出现小主元；主元消去过程中出现小主元；特征值相差大数量级。
6、特征值相差大数量级。2 Error Analysis for .bxA 近似解的误差估计及改善：近似解的误差估计及改善：设设 的近似解为的近似解为，则一般有，则一般有bxA*x0*xAbrbrxxx|*|cond(A)误差上限误差上限 改善方法：改善方法：Step 1:近似解近似解 bxA;1xStep 2:;11xAbr Step 3:;111drdA Step 4:;112dxx 若若 可被精确解出，则有可被精确解出，则有 就是精确解了。就是精确解了。1dbAxAbAxx11112)(2x经验表明经验表明：若：若 A 不是非常病态（例如：不是非常病态（例如：），），则如此迭代可达到机器精度
7、；但若则如此迭代可达到机器精度；但若 A 病态，则此算法也病态，则此算法也不能改进。不能改进。1)(Acond HW:p.66#2,#4,#53 Jacobi 法和法和 Gauss-Seidel 法法 /*Jacobi&Gauss-Seidel Iterative Methods*/Jacobi Iterative Method nnnnnnnnnnbxaxaxabxaxaxabxaxaxa.22112222212111212111 nnnnnnnnnnnnbxaxaaxbxaxaaxbxaxaax11112212122211212111.1.1.10 iia写成写成矩阵形式矩阵形式：A=LU
8、DbxULxDbxULDbxA )()(bDxULDx11)(BfJacobi 迭代阵迭代阵bDxULDxkk1)(1)1()(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;
9、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-6Step 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
10、of iterations exceeded);STOP./*unsuccessful*/bxA)0(xiinjijjijiiaXabX 1)0(TOLXXXXiini|0|max|0|1What if aii=0?迭代过程中，迭代过程中，A 的元素的元素不改变，故可以不改变，故可以事先调整事先调整好好 A 使得使得aii 0，否则否则 A不可逆不可逆。必须等必须等X(k)完全计算完全计算好了才能计算好了才能计算X(k+1)，因此，因此需要需要两组向量两组向量存储。存储。A bit wasteful,isnt it?3 Jacobi&Gauss-Seidel Iterative Methods
11、 Gauss-Seidel Iterative Method)(11)(1)(414)(313)(21211)1(1bxaxaxaxaaxknnkkkk )(12)(2)(424)(323)1(12122)1(2bxaxaxaxaaxknnkkkk )(13)(3)(434)1(232)1(13133)1(3bxaxaxaxaaxknnkkkk )(1)1(11)1(33)1(22)1(11)1(nknnnknknknnnknbxaxaxaxaax 只存一组向量即可。只存一组向量即可。写成写成矩阵形式矩阵形式：bDxUxLDxkkk1)()1(1)1()(bxUxLDkk )()1()(bLD
12、xULDxkk1)(1)1()()(BfGauss-Seidel 迭代阵迭代阵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二种方法都存在二种方法都存在收敛性问题收敛性问题。有例子表明：有例子表明：G
13、auss-Seidel法收敛时，法收敛时，Jacobi法可能法可能不收敛；而不收敛；而Jacobi法收敛时，法收敛时，Gauss-Seidel法也可能法也可能不收敛。不收敛。p.76#2 给出了例子。给出了例子。收敛性分析将在下节课讨论。收敛性分析将在下节课讨论。3 Jacobi&Gauss-Seidel Iterative MethodsLab 07.Gauss-Seidel MethodUse the Gauss-Seidel method to solve a given nn linear system with an initial approximation and a given
14、 tolerance TOL.InputThere are several sets of inputs.For each set:The 1st line contains an integer 100 n 0 which is the size of a matrix.n=1 signals the end of file.The following n lines contain the augmented matrix in the following format:The numbers are separated by spaces and new lines.The last l
15、ine of each test case contains a real number TOL,which is the tolerance for|norm,and an integer N 0 which is the maximum number of iteration.bxA 0)0(xnnnnnnbaabaabaa.1222111113 Jacobi&Gauss-Seidel Iterative MethodsOutput /*represents a space*/Each entry of the solution is to be printed as in the C f
16、printf:fprintf(outfile,%12.8fn,x);If the matrix has a zero column,print the message “Matrixhasazero column.No uniquesolutionexists.n”.If the method fails to give a solution after N iterations,print the message“Maximumnumberof iterationsexceeded.n”.If there is an entry of that is out of the range ,pr
17、int the message“No convergence.n”.The outputs of two test cases must be seperated by a blank line.Sample InputSample Input (represents a space)represents a space)31010911027041060.000000001100033101360033040.00000000110001)(kx2,2127127 Sample OutputSample Output (represents a space)represents a space)1.000000001.000000001.00000000 Matrixhasazerocolumn.Nouniquesolutionexists.
