2.1Gauss消去法和LU分解47.所以 det(U22)0,即 U22 E R(n-1)x(n-1)非奇异.由归纳假设可知,存在置换矩阵 P和 P,使得PU22P=L22U22其中L22为单位下三角矩阵，U22为非奇异上三角矩阵.取[10][10P, P=PP[ PP则有0PAPPP0L2110U12u11PL21PPTL22U22PP01010U12PL21PPTL22UPT0P0010U12P2u11PL21L220220会 LU,其中L为单位下三角矩阵，U为非奇异上三角矩阵.口由归纳法可知,结论成立。问题：第k步时，如何选取置换矩阵P和P？1.选取P和P2使得主元为剩下的矩阵中绝对值最大，这种选取方法称为“全主元Gauss消去法”，简称GECP(Gaussian elimination with completepivoting);2.选取P和P2使得主元为第k列中第k到第n个元素中，绝对值最大，这种选取方法称为“部分选主元Gauss消去法，简称GEPP（Gaussianeliminationwithpartialpivoting），此时P2=I，因此也称为列主元Gauss消去法·GECP比GEPP更稳定，但工作量太大，在实际应用中通常使用GEPP算法·GEPP算法能保证L所有的元素的绝对值都不超过1.算法2.8.部分选主元LU分解1:p=1:n：%用于记录置换矩阵2:fori=1ton-ldo%选列主元3:aki= maxi<jSnlajil4:if k~=ithen5:for j=l ton do6:%交换A的第i行与第行tmp=aij,aij=akj,akj=tmp7:end for8:tmp= p(k), p(k) = p(i),p(i) = tmp%更新置换矩阵end if9:
2.1 Gauss 消去法和 LU 分解 · 47 · 所以 det(U22) ̸= 0, 即 U22 ∈ R (n−1)×(n−1) 非奇异. 由归纳假设可知, 存在置换矩阵 P˜ 1 和 P˜ 2 使得 P˜ 1U22P˜ 2 = L˜ 22U˜ 22, 其中 L˜ 22 为单位下三角矩阵, U˜ 22 为非奇异上三角矩阵. 取 P1 = [ 1 0 0 P˜ 1 ] Pˆ 1, P2 = Pˆ 2 [ 1 0 0 P˜ 2 ] , 则有 P1AP2 = [ 1 0 0 P˜ 1 ] [ 1 0 L21 I ] [u11 U12 0 U22] [1 0 0 P˜ 2 ] = [ 1 0 P˜ 1L21 P˜ 1 ] [u11 U12 0 P˜⊺ 1 L˜ 22U˜ 22P˜⊺ 2 ] [1 0 0 P˜ 2 ] = [ 1 0 P˜ 1L21 P˜ 1 ] [1 0 0 P˜⊺ 1 L˜ 22] [u11 U12 0 U˜ 22P˜⊺ 2 ] [1 0 0 P˜ 2 ] = [ 1 0 P˜ 1L21 L˜ 22] [u11 U12P˜ 2 0 U˜ 22 ] ≜ LU, 其中 L 为单位下三角矩阵, U 为非奇异上三角矩阵. 由归纳法可知, 结论成立. □ 问题: 第 k 步时, 如何选取置换矩阵 P1 和 P2？ 1. 选取 P1 和 P2 使得主元为剩下的矩阵中绝对值最大, 这种选取方法称为 “全主元 Gauss 消去法”, 简 称 GECP (Gaussian elimination with complete pivoting); 2. 选取 P1 和 P2 使得主元为第 k 列中第 k 到第 n 个元素中, 绝对值最大, 这种选取方法称为 “部分选 主元 Gauss 消去法”, 简称 GEPP (Gaussian elimination with partial pivoting), 此时 P2 = I, 因此也称为 列主元 Gauss 消去法. • GECP 比 GEPP 更稳定, 但工作量太大, 在实际应用中通常使用 GEPP 算法. • GEPP 算法能保证 L 所有的元素的绝对值都不超过 1. 算法 2.8. 部分选主元 LU 分解 1: p = 1 : n % 用于记录置换矩阵 2: for i = 1 to n − 1 do 3: aki = maxi≤j≤n |aji| % 选列主元 4: if k ∼= i then 5: for j = 1 to n do 6: tmp = aij , aij = akj , akj = tmp % 交换 A 的第 i 行与第 k 行 7: end for 8: tmp = p(k), p(k) = p(i), p(i) = tmp % 更新置换矩阵 9: end if
:48 .第二讲线性方程组直接方法10:forj=i+ltondo%计算L的第i列11:aji=aj/au12:end for13:forj=i+l tondo14:fork=i+l tondo15:%更新A(i+1:n.i+1:najk=ajkaji*ak16:endforend for17:18: end for相应的MATLAB程序如下：MATLAB源代码2.4.部分选主元LU分解function [A,p] = PLU(A)[n,n]=size(A);3p=1:n;4for i-l:n-15[a,k]-max(abs(A(i:n,i)));6ifa==o7error('Error：第%d步的列主元为e!\n'，i);8end9k=k+i-1;10if k~=i11tmp=A(i,:): A(i,:)=A(k,:); A(k,:)=tmp;12tmp=p(i);p(i)=p(k);p(k)=tmp;13end14A(i+1:n,i)=A(i+1:n,i)/A(i,i);15A(i+1:n,i+1:n)=A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n);16endMATLAB代码见PLU.m0.0261.361.5例2.2用部分选主元LU分解求解线性方程组AE=b,其中A要求在3.43-8.525.8运算过程中保留3位有效数字解由于[a21l>[a11l，根据部分选主元LU分解算法，我们需要将第一行与第二行交换，即取P：01，然后计算A=PA的LU分解，可得1-l11=1.00, l21=a21/a11= 5.83×10-3,l22=1.00u11 = a11= 3.43, u12= a12=8.50,u22= a22- l21u12~6.13×10+4.96×10-2~6.13×10
· 48 · 第二讲 线性方程组直接方法 10: for j = i + 1 to n do 11: aji = aji/aii % 计算 L 的第 i 列 12: end for 13: for j = i + 1 to n do 14: for k = i + 1 to n do 15: ajk = ajk − aji ∗ aik % 更新 A(i + 1 : n, i + 1 : n) 16: end for 17: end for 18: end for 相应的 MATLAB 程序如下: ✞ MATLAB 源代码 2.4. 部分选主元 LU 分解 ☎ 1 function [A,p] = PLU(A) 2 [n,n]=size(A); 3 p=1:n; 4 for i=1:n-1 5 [a,k]=max(abs(A(i:n,i))); 6 if a==0 7 error('Error: 第 %d 步的列主元为 0!\n', i); 8 end 9 k=k+i-1; 10 if k~=i 11 tmp=A(i,:); A(i,:)=A(k,:); A(k,:)=tmp; 12 tmp=p(i); p(i)=p(k); p(k)=tmp; 13 end 14 A(i+1:n,i)=A(i+1:n,i)/A(i,i); 15 A(i+1:n,i+1:n)=A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n); 16 end ✝ ✆ MATLAB 代码见 PLU.m 例 2.2 用部分选主元 LU 分解求解线性方程组 Ax = b, 其中 A = [ 0.02 61.3 3.43 −8.5 ] , b = [ 61.5 25.8 ] , 要求在 运算过程中保留 3 位有效数字. 解. 由于 |a21| > |a11|, 根据部分选主元 LU 分解算法, 我们需要将第一行与第二行交换, 即取 P = [ 0 1 1 0] , 然后计算 A˜ = P A 的 LU 分解, 可得 l11 = 1.00, l21 = ˜a21/˜a11 = 5.83 × 10−3 , l22 = 1.00, u11 = ˜a11 = 3.43, u12 = ˜a12 = −8.50, u22 = ˜a22 − l21u12 ≈ 6.13 × 10 + 4.96 × 10−2 ≈ 6.13 × 10