复旦大学Matlab课件--CHAPTER2.pptVIP

Matlab Math线性方程组 Cleve Morler著 陈文斌（wbchen@fudan.edu.cn) 复旦大学2002 * Ax = b 有x = A-1b，但实际上并不显式求A-1 线性方程组 例子： 7x = 21 x = 21/7=3 如果求逆 x = 7-1× 21 = .142857 × 21 = 2.99997 这就需要一次除和一次乘，且精度更低 Backslash运算符 AX = B X = A\B 左除 XA = B X = B/A 右除 3-by-3的例子 算法矩阵表示 单位下三角阵 上三角阵 置换阵 置换阵 PA交换A的行，AP交换A的列 置换向量 p =[4 1 3 2] PA = A(p , :) Px = b x = PTb 置换方程求解 三角阵的求解 Ux = b x = zeros(n, 1); for k = n:-1:1 x(k) = b(k)/U(k,k); i = (1:k-1)’; b(i) = b(i) – x(k)*U(i,k); end x = zeros(n, 1); for k = n:-1:1 j = k+1:n x(k) = (b(k) – U(k,j)*x(j))/U(k,k); end 求出x(n)，然后把x(n)消掉 LU分解 C.F.Gauss: 高斯消去法(GE) 1955~1977: pivot选主元，舍入误差的影响 高斯消去法分为两步：向前消去和向后替换 LU分解例子 线性方程组求解 Ax = b LU = PA Ly = Pb Ux = y x = A\b 舍去和舍入 六位有效数字： 舍去浮点表示（chopped floating-pointing representation） 舍入浮点表示（rounded floating-pointing representation） 为何要选主元？ 5位有效数字 在消元时第二个方程乘上2.5×103 6.001 × 2.5 ×103 =1.50025 ×104 舍入为1.5002 ×104,再加上2.5,再一次舍入，成为 1.5004 ×104 ，则第三个方程为 -0.001x2 +(6)(0.99993) = 6.001 x3= 1.5004 ×104 / (1.5005 ×104 ) = 0.99993 1.5005 ×104 x3 = 1.5004 ×104 x2=-1.5 x1= -0.35 真解： (0, -1, 1)T 部分选主元 如果每次乘数在量上是小于或等于1的，计算出来的解就是精确的。 要使乘数小于或等于1，可以通过部分选主元来实现，选要消去列的最大元素到对角元，交换A和b的行（x不要交换）。 lutx.m : matlab内置lu分解的一个版本 bslashtx.m： matlab内置backslash运算 lugui.m：动态显示LU分解 舍入误差的影响 设计算的解为x*,真解为x = A-1b, 误差为：e = x – x*, 残量为 r = b – Ax* 六位精度算法 (选主元）得到精确解 三位精度算法(选主元） ： 分析 在三位精度的算法中，x2是由两个与舍入误差相同阶的数相除得到的，因此x2可以是任意的数。然后完全任意的x2代入第一个方程得到x1 ，所以第一方程的残量是小的。我们也可以期望第二个方程的残量也是小的，因为矩阵是接近奇异的。 部分主元的高斯消去法可以确保产生小的残量 “小的”舍入误差的阶与三个数量有关：原来矩阵悉数矩阵单元的大小，消去过程中系数矩阵单元的大小和计算解的元素的大小。如果其中一个是“大的”，残量在绝对意义上不必是小的。 即使残量是小的，但不能保证解的误差是小的。它们之间的关系可以部分由矩阵的条件数来刻画 范数 如果Ax = b, 如何测量x关于A和b的敏感性？ x =(1:4)/5 n1 = norm(x,1), n2 = norm(x) ninf = norm(x,inf) *