2.1LU分解与Gauss消去法.51.定理2.5也可以利用定理2.1的结论来证明，见习题2.3.事实上，定理2.5中的PL和PB只有一个是必需的置换矩阵的选取问题：第k步时，如何选取置换矩阵P(）和P(K)?选法一。选取P(）和P(K）使得主元为剩下的矩阵中绝对值最大,这种选取方法称为“全主元Gauss消去法”，简称GECP（Gaussianeliminationwithcompletepivoting)；选法二.选取P()和P(k)使得主元为第k列中第k到第n个元素中，绝对值最大,这种选取方法称为“部分选主元Gauss消去法”，简称GEPP（Gaussianeliminationwithpartialpivoting），此时P(k)=I,因此也称为列主元Gauss 消去法.凸（1)GECP比GEPP更稳定，但工作量太大，在实际应用中通常使用GEPP算法(2)GEPP算法能保证L所有的元素的绝对值都不超过1.算法2.8.部分选主元LU分解1:p=1:n%用于记录置换矩阵2:fork =1 to n-1 do3:[amax, ] = ,max [aik]%选列主元，其中1表示主元所在的行e<i<n4:ifl+kthen5:for j = l to n do6:%交换A的第k行与第1行atmp=akj,akj=alj,alj=atmp1end for%%更新置换矩阵Ptmp = p(k), p(k) = p(I), p(I) = Ptmp9:end if10:fori= k +1 to n do11:ak=aik/akk12:for j = k+1 to n do13:ai=aij-aak14:end forend for15:16: end for相应的MATLAB程序见LE_PLU.m0.0261.3例2.2用部分选主元LU分解求解线性方程组AT=b,其中A3.43-8.525.8
2.1 LU 分解与 Gauss 消去法 · 51 · b 定理 2.5 也可以利用定理 2.1 的结论来证明, 见习题 2.3. b 事实上, 定理 2.5 中的 PL 和 PR 只有一个是必需的. 置换矩阵的选取 问题: 第 k 步时, 如何选取置换矩阵 P (k) L 和 P (k) R ？ 选法一. 选取 P (k) L 和 P (k) R 使得主元为剩下的矩阵中绝对值最大, 这种选取方法称为 “全主元 Gauss 消去法”, 简称 GECP (Gaussian elimination with complete pivoting); 选法二. 选取 P (k) L 和 P (k) R 使得主元为第 k 列中第 k 到第 n 个元素中, 绝对值最大, 这种选取 方法称为 “部分选主元 Gauss 消去法”, 简称 GEPP (Gaussian elimination with partial pivoting), 此时 P (k) R = I, 因此也称为列主元 Gauss 消去法. b (1) GECP 比 GEPP 更稳定, 但工作量太大, 在实际应用中通常使用 GEPP 算法. (2) GEPP 算法能保证 L 所有的元素的绝对值都不超过 1. 算法 2.8. 部分选主元 LU 分解 1: p = 1 : n % 用于记录置换矩阵 2: for k = 1 to n − 1 do 3: [amax, l] = max k≤i≤n |aik| % 选列主元, 其中 l 表示主元所在的行 4: if l ̸= k then 5: for j = 1 to n do 6: atmp = akj , akj = alj , alj = atmp % 交换 A的第 k 行与第 l 行 7: end for 8: ptmp = p(k), p(k) = p(l), p(l) = ptmp % 更新置换矩阵 9: end if 10: for i = k + 1 to n do 11: aik = aik/akk 12: for j = k + 1 to n do 13: aij = aij − aikakj 14: end for 15: end for 16: end for 相应的 MATLAB 程序见 LE_PLU.m. 例 2.2 用部分选主元 LU 分解求解线性方程组 Ax = b, 其中 A = " 0.02 61.3 3.43 −8.5 # , b = " 61.5 25.8 # , 要
-52第二讲线性方程组直接方法求在运算过程中保留3位有效数字(板书）解。由于[a21l>[a11l，根据部分选主元LU分解算法，我们需要将第一行与第二行交换，即取[01[3.43-8.5]，然后计算A=PA的LU分解，即设D0.0261.310u1211A= LU9121直接比较等式两边可得u11 = a11 = 3.43, u12 = a12 = 8.5,l21 = a21/u11 ~ 5.83×10-3W22=a22-l21u12~61.3+0.0496~61.3,即1.0003.43-8.50PA~5.83×10-361.31.000解方程组Ly=Pb可得y1 = 25.8,92 ~ 61.2.解方程组Ur=y可得2= y2/u22~0.998,1= (y1 - u12r2)/u11~10.0.口与精确解21=10.2=1相比，数值解具有3位有效数字2.1.5矩阵求逆我们可以通过部分选主元LU分解来计算矩阵的逆.设PA=LU，则A-1 = U-1L-1P,等价于求解下面2n个三角线性方程组Lyi=Pei,Uai=i,i=1,2,...,n.思考：也可以分别计算L-1和U-1，然后相乘.哪种方法划算？2.1.6分块LU分解为了提高计算效率，实际计算中通常采用分块LU分解，即1[A1A12..AipUiiU12UipL2112U22. U2pA21 A22* A2pA:= LU,:....Up]Ap1Ap2... AppLp1..Lp.p-1Ip其中Ai是非奇异的方阵
· 52 · 第二讲 线性方程组直接方法 求在运算过程中保留 3 位有效数字. (板书) 解. 由于 |a21| > |a11|, 根据部分选主元 LU 分解算法, 我们需要将第一行与第二行交换, 即取 P = " 0 1 1 0# , 然后计算 A˜ = P A = " 3.43 −8.5 0.02 61.3 # 的 LU 分解, 即设 A˜ = LU = " 1 0 l21 1 # "u11 u12 0 u22# . 直接比较等式两边可得 u11 = ˜a11 = 3.43, u12 = ˜a12 = −8.5, l21 = ˜a21/u11 ≈ 5.83 × 10−3 , u22 = ˜a22 − l21u12 ≈ 61.3 + 0.0496 ≈ 61.3, 即 P A ≈ " 1.00 0 5.83 × 10−3 1.00# "3.43 −8.50 0 61.3 # . 解方程组 Ly = P b 可得 y1 = 25.8, y2 ≈ 61.2. 解方程组 Ux = y 可得 x2 = y2/u22 ≈ 0.998, x1 = (y1 − u12x2)/u11 ≈ 10.0. 与精确解 x1 = 10, x2 = 1 相比, 数值解具有 3 位有效数字. □ 2.1.5 矩阵求逆 我们可以通过部分选主元 LU 分解来计算矩阵的逆. 设 P A = LU, 则 A −1 = U −1L −1P, 等价于求解下面 2n 个三角线性方程组 Lyi = P ei , Uxi = yi , i = 1, 2, . . . , n. 思考：也可以分别计算 L −1 和 U −1 , 然后相乘. 哪种方法划算? 2.1.6 分块 LU 分解 为了提高计算效率, 实际计算中通常采用 分块 LU 分解, 即 A = A11 A12 · · · A1p A21 A22 · · · A2p . . . . . . . . . Ap1 Ap2 · · · App = I1 L21 I2 . . . . . . . . . Lp1 · · · Lp,p−1 Ip U11 U12 · · · U1p U22 · · · U2p . . . . . . Upp = LU, 其中 Aii 是非奇异的方阵
2.1LU分解与Gauss消去法.53.与定理2.1相类似，对于分块LU分解，我们有下面的结论定理2.6(分块LU分解的存在性和唯一性）矩阵A存在唯一分块LU分解的充要条件是A的所有顺序分块主子矩阵都非奇异（留作练习）几点注记日分块LU分解不是LU分解，若A存在LU分解，则一定存在分块LU分解，反之不成立日分块LU分解与普通LU分解的总运算量是一样的，但分块LU分解可以借助3级BLAS运算，因此效率更高.日在计算分块LU分解过程中需要求解一系列小规模的线性方程组（以Ai为系数矩阵），可以采用普通（选主元）LU分解方法求解
2.1 LU 分解与 Gauss 消去法 · 53 · 与定理 2.1 相类似, 对于分块 LU 分解, 我们有下面的结论. 定理 2.6 (分块 LU 分解的存在性和唯一性) 矩阵 A 存在唯一分块 LU 分解的充要条件是 A 的 所有顺序分块主子矩阵都非奇异. (留作练习) 几点注记 分块 LU 分解不是 LU 分解, 若 A 存在 LU 分解, 则一定存在分块 LU 分解, 反之不成立. 分块 LU 分解与普通 LU 分解的总运算量是一样的, 但分块 LU 分解可以借助 3 级 BLAS 运算, 因此效率更高. 在计算分块 LU 分解过程中需要求解一系列小规模的线性方程组 (以 Aii 为系数矩阵), 可 以采用普通 (选主元) LU 分解方法求解
:54.第二讲线性方程组直接方法2.2特殊方程组的求解如果系数矩阵具有一定的特殊结构或性质，则可以充分利用这些特殊结构或性质来设计高效的数值算法.本节考虑以下特殊方程组的求解：·对称正定矩阵·对称不定矩阵·三对角矩阵·带状矩阵·Toeplitz矩阵2.2.1对称正定线性方程组考虑线性方程组Ar=b其中AERnxn对称正定的定理2.7(Cholesky分解）设AERnxn对称正定，则存在唯一的对角线元素为正的下三角矩阵L,使得A= LLT.(板书)该分解称为Cholesky分解证明.首先证明存在性，用数学归纳法当n=1时，由A的对称正定性可知a11>0.取l11=Va11即可，假定结论对所有不超过n-1阶的对称正定矩阵都成立.设AERnxn是n阶对称正定，则A可分解为0Va1ia11A12Va111A:A1-A2I1Af2A2210Y0其中A22=A22－AI2A12/a11.由于A对称正定，所以也对称正定，故A22是n一1阶0A22对称正定矩阵.根据归纳假设，存在唯一的对角线元素为正的下三角矩阵L，使得A22=LLT.令[OL=0易知，L是对角线元素均为正的下三角矩阵，且VanioVa1i[]LLT 1=A3-AI1由归纳法可知，对任意对称正定实矩阵A，都存在一个对角线元素为正的下三角矩阵L，使得A= LLT.口唯一性可以采用反证法，留做练习
· 54 · 第二讲 线性方程组直接方法 2.2 特殊方程组的求解 如果系数矩阵具有一定的特殊结构或性质, 则可以充分利用这些特殊结构或性质来设计高效 的数值算法. 本节考虑以下特殊方程组的求解: • 对称正定矩阵 • 对称不定矩阵 • 三对角矩阵 • 带状矩阵 • Toeplitz 矩阵 2.2.1 对称正定线性方程组 考虑线性方程组 Ax = b 其中 A ∈ R n×n 对称正定的. 定理 2.7 (Cholesky 分解) 设 A ∈ R n×n 对称正定, 则存在唯一的对角线元素为正的下三角矩阵 L, 使得 A = LLT . 该分解称为 Cholesky 分解. (板书) 证明. 首先证明存在性, 用数学归纳法. 当 n = 1 时, 由 A 的对称正定性可知 a11 > 0. 取 l11 = √ a11 即可. 假定结论对所有不超过 n − 1 阶的对称正定矩阵都成立. 设 A ∈ R n×n 是 n 阶对称正定, 则 A 可分解为 A = " a11 A12 AT 12 A22# = " √ a11 0 √ 1 a11 AT 12 I # "1 0 0 A˜ 22# " √ a11 0 √ 1 a11 AT 12 I #T , 其中 A˜ 22 = A22 − AT 12A12/a11. 由于 A 对称正定, 所以 " 1 0 0 A˜ 22# 也对称正定, 故 A˜ 22 是 n − 1 阶 对称正定矩阵. 根据归纳假设, 存在唯一的对角线元素为正的下三角矩阵 L˜, 使得 A˜ 22 = L˜L˜T. 令 L = " √ a11 0 √ 1 a11 AT 12 I # "1 0 0 L˜ # = " √ a11 0 √ 1 a11 AT 12 L˜ # . 易知, L 是对角线元素均为正的下三角矩阵, 且 LLT = " √ a11 0 √ 1 a11 AT 12 I # "1 0 0 L˜ # "1 0 0 L˜T # " √ a11 0 √ 1 a11 AT 12 I #T = A. 由归纳法可知, 对任意对称正定实矩阵 A, 都存在一个对角线元素为正的下三角矩阵 L, 使得 A = LLT . 唯一性可以采用反证法, 留做练习. □