华东师范大学：《矩阵计算》课程教学资源（课堂讲义）第二讲 线性方程组直接方法

第二讲线性方程组直接方法一般来说，求解线性方程组的数值方法可以分为两类：直接法与选代法.本章介绍直接法，即Gauss消去法.直接法相对比较稳定，因此在工程计算中很受欢迎.但由于运算量是O（n3），当问题规模较大时，时间会很长（这里n表示未知量的个数）.目前，直接法主要用于小规模或中等规模线性方程组的数值求解2.1Gauss消去法和LU分解2.1.1LU分解考虑线性方程组Ar= b,(2.1)其中AERnxn非奇异，bER"为给定的右端项.Gauss消去法本质上就是对系数矩阵A进行LU分解，即将A分解成两个矩阵的乘积A = LU,(2.2)其中L是单位下三角矩阵，U为非奇异上三角矩阵.这个分解就称为LU分解假定矩阵A存在LU分解（2.2)，则方程组（2.1）就转化为求解下面两个三角方程组Ly = b,Ur=y显然，这两个方程组都非常容易求解基于LU分解的Gauss消去法描述如下：算法2.1.Gauss消去法1：将A进行LU分解：A=LU，其中L为单位下三角矩阵，U为非奇异上三角矩阵2:利用向前回代,求解Ly=b,即得y=L-1b3:利用向后回代，求解U=9,即得=U-1y=（LU)-1b=A-1b我们知道，当系数矩阵A非奇异时，方程组（2.1）总是存在唯一解但是，并不是每个非奇异矩阵都存在LU分解定理2.1（LU分解的存在性和唯一性）设AERnXn.则存在唯一的单位下三角矩阵L和非奇异上三角矩阵U，使得A=LU的充要条件是A的所有顺序主子矩阵A=A(1：k,1：k)都非奇异k=1,2,...,n.37

第二讲 线性方程组直接方法 一般来说, 求解线性方程组的数值方法可以分为两类: 直接法与迭代法. 本章介绍直接法, 即 Gauss 消去法. 直接法相对比较稳定, 因此在工程计算中很受欢迎. 但由于运算量是 O(n 3 ), 当问题规模较大 时, 时间会很长 (这里 n 表示未知量的个数). 目前, 直接法主要用于小规模或中等规模线性方程组的数 值求解. 2.1 Gauss 消去法和 LU 分解 2.1.1 LU 分解 考虑线性方程组 Ax = b, (2.1) 其中 A ∈ R n×n 非奇异, b ∈ R n 为给定的右端项. Gauss 消去法本质上就是对系数矩阵 A 进行 LU 分解, 即将 A 分解成两个矩阵的乘积 A = LU, (2.2) 其中 L 是单位下三角矩阵, U 为非奇异上三角矩阵. 这个分解就称为 LU 分解. 假定矩阵 A 存在 LU 分解 (2.2), 则方程组 (2.1) 就转化为求解下面两个三角方程组 { Ly = b, Ux = y. 显然, 这两个方程组都非常容易求解. 基于 LU 分解的 Gauss 消去法描述如下: 算法 2.1. Gauss 消去法 1: 将 A 进行 LU 分解: A = LU, 其中 L 为单位下三角矩阵, U 为非奇异上三角矩阵; 2: 利用向前回代, 求解 Ly = b, 即得 y = L −1 b; 3: 利用向后回代, 求解 Ux = y, 即得 x = U −1y = (LU) −1 b = A−1 b. 我们知道, 当系数矩阵 A 非奇异时, 方程组 (2.1) 总是存在唯一解. 但是, 并不是每个非奇异矩阵都 存在 LU 分解. 定理 2.1 (LU 分解的存在性和唯一性) 设 A ∈ R n×n. 则存在唯一的单位下三角矩阵 L 和非奇异上 三角矩阵 U, 使得 A = LU 的充要条件是 A 的所有顺序主子矩阵 Ak = A(1 : k, 1 : k) 都非奇异, k = 1, 2, . . . , n. 37

:38.第二讲线性方程组直接方法证明.必要性：设A11是A的k阶顺序主子矩阵，将A=LU写成分块形式LnU12A11A12LnUiiL11U11U12L21U11L21U12+L22U22A21A22L210U22可得A11=LiiUi1.由于L11和Uii均非奇异，所以A11也非奇异充分性：用归纳法当n=1时，结论显然成立假设结论对n－1阶矩阵都成立，即对任意n－1阶矩阵，如果其所有的顺序主子矩阵都非奇异，则存在LU分解考虑n阶的矩阵A,写成分块形式21其中A11ER(n-1)x(n-1)是A的n-1阶顺序主子矩阵.由归纳假设可知,A11存在LU分解,即存在单位下三角矩阵L11和非奇异上三角矩阵U11使得A1 = LiuUi.令L21=A21U-l,U12=LiA12,U22=A22-L21U12则[A11.A12[LuUn[L1L11U120U11ALU1A21A22L2iUi1U22+L21U12L210U22易知U非奇异，所以A存在LU分解，下面证明唯一性设A存在两个不同的LU分解A=LU=LU其中L和工为单位下三角矩阵，U和ü为非奇异上三角矩阵.则有L-1L=UU-1该等式左边为下三角矩阵，右边为上三角矩阵，所以只能是对角矩阵，由于单位下三角矩阵的逆仍然是单位下三角矩阵，所以L-1L的对角线元素全是1，故L-1i= I,即i=L,U=U.口由归纳法可知，结论成立

· 38 · 第二讲 线性方程组直接方法 证明. 必要性: 设 A11 是 A 的 k 阶顺序主子矩阵, 将 A = LU 写成分块形式 [ A11 A12 A21 A22] = [ L11 0 L21 L22] [U11 U12 0 U22] = [ L11U11 L11U12 L21U11 L21U12 + L22U22] . 可得 A11 = L11U11. 由于 L11 和 U11 均非奇异, 所以 A11 也非奇异. 充分性: 用归纳法. 当 n = 1 时, 结论显然成立. 假设结论对 n − 1 阶矩阵都成立, 即对任意 n − 1 阶矩阵, 如果其所有的顺序主子矩阵都非奇异, 则 存在 LU 分解. 考虑 n 阶的矩阵 A, 写成分块形式 A = [ A11 A12 A21 A22] , 其中 A11 ∈ R (n−1)×(n−1) 是 A 的 n − 1 阶顺序主子矩阵. 由归纳假设可知, A11 存在 LU 分解, 即存在单 位下三角矩阵 L11 和非奇异上三角矩阵 U11 使得 A11 = L11U11. 令 L21 = A21U −1 11 , U12 = L −1 11 A12, U22 = A22 − L21U12, 则 A = [ A11 A12 A21 A22] = [ L11U11 L11U12 L21U11 U22 + L21U12] = [ L11 0 L21 1 ] [U11 U12 0 U22] ≜ LU. 易知 U 非奇异, 所以 A 存在 LU 分解. 下面证明唯一性. 设 A 存在两个不同的 LU 分解: A = LU = L˜U, ˜ 其中 L 和 L˜ 为单位下三角矩阵, U 和 U˜ 为非奇异上三角矩阵. 则有 L −1L˜ = UU˜ −1 , 该等式左边为下三角矩阵, 右边为上三角矩阵, 所以只能是对角矩阵. 由于单位下三角矩阵的逆仍然是 单位下三角矩阵, 所以 L −1L˜ 的对角线元素全是 1, 故 L −1L˜ = I, 即 L˜ = L, U˜ = U. 由归纳法可知, 结论成立. □

2.1Gauss消去法和LU分解·39 .2.1.2LU分解的实现给定一个矩阵a11α12aina21a22a2nERnxnA=:..[anl...an2ann我们可以通过矩阵初等变换来构造A的LU分解，·第一步：假定a11≠0,构造矩阵[100...012110..0Qil01.0其中l31Li =lil:.i=2,3,...,n.a11.[ni00.1易知L的逆为0001121100Lil =l3101...0:..-ln1001...用L-1左乘A,并将所得到的矩阵记为A(1),则a11α12ain(1)(1)0022A(1) = L-"A.....a0即左乘L-1后，A的第一列中除第一个元素外其它都变为0.·第二步：类似地，我们可以将上面的操作作用在A(1）的子矩阵A(1)(2：n.2：n）上，将其第一列除第一个元素外都变为0.也就是说，假定a≠0.构造矩阵100700100al)..13210C其中L2 =1.2 ...：..01n20...1

2.1 Gauss 消去法和 LU 分解 · 39 · 2.1.2 LU 分解的实现 给定一个矩阵 A =        a11 a12 · · · a1n a21 a22 · · · a2n . . . . . . an1 an2 · · · ann        ∈ R n×n . 我们可以通过矩阵初等变换来构造 A 的 LU 分解. • 第一步: 假定 a11 ̸= 0, 构造矩阵 L1 =          1 0 0 · · · 0 l21 1 0 · · · 0 l31 0 1 · · · 0 . . . . . . ln1 0 0 · · · 1          , 其中 li1 = ai1 a11 , i = 2, 3, . . . , n. 易知 L1 的逆为 L −1 1 =          1 0 0 · · · 0 −l21 1 0 · · · 0 −l31 0 1 · · · 0 . . . . . . −ln1 0 0 · · · 1          . 用 L −1 1 左乘 A, 并将所得到的矩阵记为 A(1) , 则 A (1) = L −1 1 A        a11 a12 · · · a1n 0 a (1) 22 · · · a (1) 2n . . . . . . . . . 0 a (1) n2 · · · a (1) nn        . 即左乘 L −1 1 后, A 的第一列中除第一个元素外其它都变为 0. • 第二步: 类似地, 我们可以将上面的操作作用在 A(1) 的子矩阵 A(1)(2 : n, 2 : n) 上, 将其第一列除 第一个元素外都变为 0. 也就是说, 假定 a (1) 22 ̸= 0, 构造矩阵 L2 =          1 0 0 · · · 0 0 1 0 · · · 0 0 l32 1 · · · 0 . . . . . . . . . 0 ln2 0 · · · 1          , 其中 li2 = a (1) i2 a (1) 22 , i = 3, 4, . . . , n

:40 .第二讲线性方程组直接方法用L-1左乘A(1)并将所得到的矩阵记为A(2)，则[a11a12aina13.e号ageen000A(2) = L-1 A(1) = L-1L-1 A =...:aea2]Lo0....依此类推，假定ak-1)0(k=3,4,n-1),则我们可以构造一系列的矩阵L3,L4....，Ln-1,使得[a11α12α13ain.ago0...00L---. L?"L-'A -..........alm-1)]L 000...为一个上三角矩阵.我们将这个上三角矩阵记为U，并记n000..10l21...013113210...L = L,L2-.- Ln-1 -(2.3)......"..[n1In3...In21则可得A= LU,这就是A的LU分解将上面的过程写成算法，描述如下：算法2.2.LU分解1:SetL=I,U=0%将L设为单位矩阵，U设为零矩阵2:fork=1ton-1do3:fori=k+1tondo4:%计算L的第K列li=ai/akkend for5:6:for j=k ton do7:%计算U的第人行ukj = akjend for8:for i=k+1 to n do9:for j=k+1 to ndo10:11:%更新A(k+1:n,k+1:n)aij = aij -likukj12:end for

· 40 · 第二讲 线性方程组直接方法 用 L −1 2 左乘 A(1) , 并将所得到的矩阵记为 A(2) , 则 A (2) = L −1 2 A (1) = L −1 2 L −1 1 A =          a11 a12 a13 · · · a1n 0 a (1) 22 a (1) 23 · · · a (1) 2n 0 0 a (2) 33 · · · a (2) 3n . . . . . . . . . . . . 0 0 a (2) n3 · · · a (2) nn          . • 依此类推, 假定 a (k−1) kk ̸= 0 (k = 3, 4, . . . , n − 1), 则我们可以构造一系列的矩阵 L3, L4, . . . , Ln−1, 使得 L −1 n−1 · · ·L −1 2 L −1 1 A =          a11 a12 a13 · · · a1n 0 a (1) 22 a (1) 23 · · · a (1) 2n 0 0 a (2) 33 · · · a (2) 3n . . . . . . . . . . . . 0 0 0 · · · a (n−1) nn          为一个上三角矩阵. 我们将这个上三角矩阵记为 U, 并记 L = L1L2 · · ·Ln−1 =          1 0 0 · · · 0 l21 1 0 · · · 0 l31 l32 1 · · · 0 . . . . . . . . . ln1 ln2 ln3 · · · 1          , (2.3) 则可得 A = LU, 这就是 A 的 LU 分解. 将上面的过程写成算法, 描述如下: 算法 2.2. LU 分解 1: Set L = I, U = 0 % 将 L 设为单位矩阵, U 设为零矩阵 2: for k = 1 to n − 1 do 3: for i = k + 1 to n do 4: lik = aik/akk % 计算 L 的第 k 列 5: end for 6: for j = k to n do 7: ukj = akj % 计算 U 的第 k 行 8: end for 9: for i = k + 1 to n do 10: for j = k + 1 to n do 11: aij = aij − likukj % 更新 A(k + 1 : n, k + 1 : n) 12: end for

2.1Gauss消去法和LU分解41.13:end for14:end forGauss消去法的运算量由算法2.2可知，LU分解的运算量（含加减乘除）为n-17121+Z(n -i+2(n-i) =n3+0(n2)13j=i+1k=i+1i=1i=i+1n3 +0(n2)由于回代过程的运算量为O(n2)，所以Gauss消去法的总运算量为G十评价算法的一个主要指标是执行时间，但这依赖于计算机硬件和编程技巧等，因此直接给出算法执行时间是不太现实的.所以我们通常是统计算法中算术运算（加减乘除）的次数在数值算法中，大多仅仅涉及加减乘除和开方运算，一般地，加减运算次数与乘法运算次数具有相同的量级而除法运算和开方运算次数具有更低的量级十为了尽可能地减少运算量，在实际计算中，数，向量和矩阵做乘法运算时的先后执行次序为：先计算数与向量的乘法，然后计算矩阵与向量的乘法，最后才计算矩阵与矩阵的乘法矩阵L和U的存储当A的第i列被用于计算L的第i列后，在后面的计算中不再被使用.同样地，A的第i行被用于计算U的第行后，在后面的计算中也不再被使用，因此，为了节省存储空间，我们可以在计算过程中将L的第i列存放在A的第i列，将U的第i行存放在A的第i行，这样就不需要另外分配空间存储L和U.计算结束后，A的上三角部分为U，其绝对下三角部分为L的绝对下三角部分（L的对角线全部为1，不需要存储).此时算法可以描述为：算法2.3.LU分解（用A存储L和U)l:fork=1ton-1do2:fori=k+1tondo3:aik = aik/akk4:for j=k+1 to n do5:aij=aj-aiakj6:endforend for7:8:endfor十根据指标的循环次序，算法2.3也称为KIJ型LU分解.在实际计算中，我们一般不建议使用这个算法.因为对于指标k的每次循环，都需要更新A的第k十1至第n行.这种反复读取数据的做法会使得计算效率大大降低

2.1 Gauss 消去法和 LU 分解 · 41 · 13: end for 14: end for Gauss 消去法的运算量 由算法 2.2 可知, LU 分解的运算量 (含加减乘除) 为 n∑−1 i=1   ∑n j=i+1 1 + ∑n j=i+1 ∑n k=i+1 2   = n∑−1 i=1 ( n − i + 2(n − i) 2 ) = 2 3 n 3 + O(n 2 ). 由于回代过程的运算量为 O(n 2 ), 所以 Gauss 消去法的总运算量为 2 3 n 3 + O(n 2 ). † 评价算法的一个主要指标是执行时间, 但这依赖于计算机硬件和编程技巧等, 因此直接给出算法 执行时间是不太现实的. 所以我们通常是统计算法中算术运算 (加减乘除) 的次数. 在数值算法 中, 大多仅仅涉及加减乘除和开方运算. 一般地, 加减运算次数与乘法运算次数具有相同的量级, 而除法运算和开方运算次数具有更低的量级. † 为了尽可能地减少运算量, 在实际计算中, 数, 向量和矩阵做乘法运算时的先后执行次序为: 先计 算数与向量的乘法, 然后计算矩阵与向量的乘法, 最后才计算矩阵与矩阵的乘法. 矩阵 L 和 U 的存储 当 A 的第 i 列被用于计算 L 的第 i 列后, 在后面的计算中不再被使用. 同样地, A 的第 i 行被用于 计算 U 的第 i 行后, 在后面的计算中也不再被使用. 因此, 为了节省存储空间, 我们可以在计算过程中将 L 的第 i 列存放在 A 的第 i 列, 将 U 的第 i 行存放在 A 的第 i 行, 这样就不需要另外分配空间存储 L 和 U. 计算结束后, A 的上三角部分为 U, 其绝对下三角部分为 L 的绝对下三角部分 (L 的对角线全部为 1, 不需要存储). 此时算法可以描述为： 算法 2.3. LU 分解 (用 A 存储 L 和 U) 1: for k = 1 to n − 1 do 2: for i = k + 1 to n do 3: aik = aik/akk 4: for j = k + 1 to n do 5: aij = aij − aikakj 6: end for 7: end for 8: end for † 根据指标的循环次序, 算法 2.3 也称为 KIJ 型 LU 分解. 在实际计算中, 我们一般不建议使用这个 算法. 因为对于指标 k 的每次循环, 都需要更新 A 的第 k + 1 至第 n 行. 这种反复读取数据的做 法会使得计算效率大大降低

第二讲线性方程组直接方法: 42对于按行存储的数据结构，我们一般采用后面介绍的IKI型LU分解MATLAB源代码2.1.LU分解的MATLAB代码（KII型）function A=mylu(A)2n=size(A,1);3for k=1:n-14if A(k,k) ==5fprintf('Error:A(%d,%d)=0!\n',k,k);6return;7end8for i=-k+1:n9A(i,k)=A(i,k)/A(k,k);10for j-k+1:n11A(i,j)=A(i,j)-A(i,k)*A(k,j);12end13end14end为了充分利用MATLAB的向量运算优势，提高运算效率，上面的程序可改写为MATLAB源代码2.2.LU分解（KII型）function Amylu(A)2n=size(A,1);3for k=1:n-14ifA(k,k)=05fprintf('Error:A(%d,%d)=o!\n',k,k);6return;7end8A(k+1:n,k)=A(k+1:n,k)/A(k,k);9A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);10end2.1.3IKJ型LU分解如果数据是按行存储的，如CIC++，我们一般采用下面的IKI型LU分解算法2.4.LU分解(IKJ型)l:fori=2tondo2:fork=ltoi-ldo3:aik=ai/akk4:forj=k+ltondo5:aij=aij -aikakj

· 42 · 第二讲 线性方程组直接方法 对于按行存储的数据结构, 我们一般采用后面介绍的 IKJ 型 LU 分解. MATLAB 源代码 2.1. LU 分解的 MATLAB 代码 (KIJ 型) ✞ ☎ 1 function A = mylu(A) 2 n=size(A,1); 3 for k=1:n-1 4 if A(k,k) == 0 5 fprintf('Error: A(%d,%d)=0!\n', k, k); 6 return; 7 end 8 for i=k+1:n 9 A(i,k)=A(i,k)/A(k,k); 10 for j=k+1:n 11 A(i,j)=A(i,j)-A(i,k)*A(k,j); 12 end 13 end 14 end ✝ ✆ 为了充分利用 MATLAB 的向量运算优势, 提高运算效率, 上面的程序可改写为 MATLAB 源代码 2.2. LU 分解 (KIJ 型) ✞ ☎ 1 function A = mylu(A) 2 n=size(A,1); 3 for k=1:n-1 4 if A(k,k) == 0 5 fprintf('Error: A(%d,%d)=0!\n', k, k); 6 return; 7 end 8 A(k+1:n,k)=A(k+1:n,k)/A(k,k); 9 A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); 10 end ✝ ✆ 2.1.3 IKJ 型 LU 分解 如果数据是按行存储的, 如 C/C++, 我们一般采用下面的 IKJ 型 LU 分解. 算法 2.4. LU 分解 (IKJ 型) 1: for i = 2 to n do 2: for k = 1 to i − 1 do 3: aik = aik/akk 4: for j = k + 1 to n do 5: aij = aij − aikakj

2.1Gauss消去法和LU分解43 -end for6:end for7:8:end for上述算法可以用下图来描述AccessedbutnotmodifiedAccessedandmodifiedNotaccessed如果数据是按列存储的，如FORTRAN或MATLAB，则该怎样设计算法？2.1.4待定系数法计算LU分解我们也可以利用待定系数法来实现矩阵的LU分解.假设A存在LU分解，即A=LU，或1aina11a12a13u12u13uinu11..l211a2122a2nu22u23a23u2n131132a3132a33a3nu33...u2n.....+·[in]ann]In2Unn[anlan2an3通过比较等式两边的元素来计算L和U中的各元素的值.具体计算过程如下：(1）比较等式两边的第一行，可得uj=alj,j=1,2,...,n.再比较等式两边的第一列，可得aii =liu1 → lii = aii/u1l, i= 2,3,--,n.(2）比较等式两边的第二行，可得a2j=l21uj→u2j=a2j-l21uij，j=2,3,..,n.再比较等式两边的第二列，可得ai2=li1u12+li2u22→li1=(ai2-liiu12)/u22,i=3,4,...,n

2.1 Gauss 消去法和 LU 分解 · 43 · 6: end for 7: end for 8: end for 上述算法可以用下图来描述. ♣ 如果数据是按列存储的, 如 FORTRAN 或 MATLAB, 则该怎样设计算法? 2.1.4 待定系数法计算 LU 分解 我们也可以利用待定系数法来实现矩阵的 LU 分解. 假设 A 存在 LU 分解, 即 A = LU, 或          a11 a12 a13 · · · a1n a21 a22 a23 · · · a2n a31 a32 a33 · · · a3n . . . . . . . . . an1 an2 an3 · · · ann          =          1 l21 1 l31 l32 1 . . . . . . ln1 ln2 · · · ln,n−1 1                   u11 u12 u13 · · · u1n u22 u23 · · · u2n u33 · · · u2n . . . . . . unn          . 通过比较等式两边的元素来计算 L 和 U 中的各元素的值. 具体计算过程如下: (1) 比较等式两边的第一行, 可得 u1j = a1j , j = 1, 2, . . . , n. 再比较等式两边的第一列, 可得 ai1 = li1u11 ⇒ li1 = ai1/u11, i = 2, 3, . . . , n. (2) 比较等式两边的第二行, 可得 a2j = l21u1j ⇒ u2j = a2j − l21uij , j = 2, 3, . . . , n. 再比较等式两边的第二列, 可得 ai2 = li1u12 + li2u22 ⇒ li1 = (ai2 − li1u12)/u22, i = 3, 4, . . . , n

: 44 .第二讲线性方程组直接方法(3）以此类推，第k步时，比较等式两边的第行，可得ukj=akj-(lkiuij+...+lk,k-1uk-1j),j=k,k+l,...,n.比较等式两边的第列，可得lik=(aik-liuik-...-li,k-iuk-1,k)/ukk,i=k+l,k+2,...,n直到第n步，即可计算出L和U的所有元素同样，我们可以利用A来存储L和U.算法描述如下：算法2.5.LU分解(待定系数法或Doolittle法)I:fork=ltondok-12:i=k.k+1.....n>ak= akjakiaij,i=1k-113:i=k+1,k+2,...,naikaikaijajkakkj=l4:endfor相应的MATLAB程序为：MATLAB源代码2.3.待定系数法LU分解functionA=mylu2(A)2[n,n]=size(A);3for k=1:n4A(k,k)=A(k,k)-A(k,1:k-1)*A(1:k-1,k):5if (A(k,k)==0)6fprintf('Error:A(%d,%d)=0!In', i,i); return;7end8A(k,k+1:n)=A(k,k+1:n)-A(k,1:k-1)*A(1:k-1,k+1:n):9A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k);10end2.1.5三角方程求解得到A的LU分解后，我们最后需要用回代法求解两个三角方程组Ly=b,Ua=y,其中L是单位下三角矩阵，U为非奇异上三角矩阵下面是关于一般下三角方程组的求解算法（行存储方式）

· 44 · 第二讲 线性方程组直接方法 (3) 以此类推, 第 k 步时, 比较等式两边的第 k 行, 可得 ukj = akj − (lk1u1j + · · · + lk,k−1uk−1,j ), j = k, k + 1, . . . , n. 比较等式两边的第 k 列, 可得 lik = (aik − li1u1k − · · · − li,k−1uk−1,k)/ukk, i = k + 1, k + 2, . . . , n. 直到第 n 步, 即可计算出 L 和 U 的所有元素. 同样, 我们可以利用 A 来存储 L 和 U. 算法描述如下: 算法 2.5. LU 分解 (待定系数法或 Doolittle 法) 1: for k = 1 to n do 2: akj = akj − k ∑−1 i=1 akiaij , j = k, k + 1, . . . , n 3: aik = 1 akk  aik − k ∑−1 j=1 aijajk  , i = k + 1, k + 2, . . . , n 4: end for 相应的 MATLAB 程序为: ✞ MATLAB 源代码 2.3. 待定系数法 LU 分解 ☎ 1 function A = mylu2(A) 2 [n,n]=size(A); 3 for k=1:n 4 A(k,k)=A(k,k)-A(k,1:k-1)*A(1:k-1,k); 5 if (A(k,k)==0) 6 fprintf('Error: A(%d,%d)=0!\n', i,i); return; 7 end 8 A(k,k+1:n)=A(k,k+1:n)-A(k,1:k-1)*A(1:k-1,k+1:n); 9 A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k); 10 end ✝ ✆ 2.1.5 三角方程求解 得到 A 的 LU 分解后, 我们最后需要用回代法求解两个三角方程组 Ly = b, Ux = y, 其中 L 是单位下三角矩阵, U 为非奇异上三角矩阵. 下面是关于一般下三角方程组的求解算法 (行存储方式)

2.1Gauss消去法和LU分解45 .算法2.6.向前回代求解Ly=b（行存储方式）1: 91 = b /l112:fori=2:ndoforj=l:i-1do3:4:bi=bi-lijyi5:end for6:yi=bi/li7:end for如果数据是按列存储的，则采用列存储方式效率会高一些，下面是按列存储方式求解上三角方程组.算法2.7.向后回代求解U=y(列存储方式)l: En = Yn/unn2: fori= n-1:-1:1doforj=l:ido3:4:Yj = yj - Ti+1ui+1.5end for5:6:Ti=yi/ui7:end for+这两个算法的运算量均为n2+O（n）以上两个算法都是向后稳定的（componentwisebackwardstable)[41]2.1.6选主元LU分解在LU分解算法2.2中,我们称ak-1)为主元.如果ak-1)=0,则算法就无法进行下去.即使ak-1)不为零，但如果[ak-1)I的值很小，由于舍人误差的原因，也可能会给计算结果带来很大的误差。此时我们就需要通过选主元来解决这个问题61.361.50.02例2.1用LU分解求解线性方程组Ar=b,其中A要求在运算过程中3.4325.8-8.5保留3位有效数字解，根据LU分解算法2.2，我们可得l11 = 1.00, l21=a21/a11=1.72×102,l22= 1.00u11=α11=2.00×10-2,u12=a12=6.13×10U22=22- l21u12~-8.5-1.05×104~-1.05×104

2.1 Gauss 消去法和 LU 分解 · 45 · 算法 2.6. 向前回代求解 Ly = b (行存储方式) 1: y1 = b1/l11 2: for i = 2 : n do 3: for j = 1 : i − 1 do 4: bi = bi − lijyj 5: end for 6: yi = bi/lii 7: end for 如果数据是按列存储的, 则采用列存储方式效率会高一些. 下面是按列存储方式求解上三角方程 组. 算法 2.7. 向后回代求解 Ux = y (列存储方式) 1: xn = yn/unn 2: for i = n − 1 : −1 : 1 do 3: for j = 1 : i do 4: yj = yj − xi+1ui+1,j 5: end for 6: xi = yi/uii 7: end for † 这两个算法的运算量均为 n 2 + O(n). 以上两个算法都是向后稳定的 (componentwise backward stable) [41]. 2.1.6 选主元 LU 分解 在 LU 分解算法 2.2 中, 我们称 a (k−1) kk 为主元. 如果 a (k−1) kk = 0, 则算法就无法进行下去. 即使 a (k−1) kk 不为零, 但如果 |a (k−1) kk | 的值很小, 由于舍入误差的原因, 也可能会给计算结果带来很大的误差. 此时我 们就需要通过选主元来解决这个问题. 例 2.1 用 LU 分解求解线性方程组 Ax = b, 其中 A = [ 0.02 61.3 3.43 −8.5 ] , b = [ 61.5 25.8 ] , 要求在运算过程中 保留 3 位有效数字. 解. 根据 LU 分解算法 2.2, 我们可得 l11 = 1.00, l21 = a21/a11 = 1.72 × 102 , l22 = 1.00, u11 = a11 = 2.00 × 10−2 , u12 = a12 = 6.13 × 10, u22 = a22 − l21u12 ≈ −8.5 − 1.05 × 104 ≈ −1.05 × 104

: 46 .第二讲线性方程组直接方法即[2.00×10-21.0006.12 × 1001.72×102-1.05×1041.00解方程组Ly=b可得91=6.15×10,92=b2-l213/1~-1.06×104解方程组Ur=y可得T2=92/u22~1.01,1=(y1-u12*2)/u11~-0.413/u11～-20.7口易知，方程的精确解为T1=10.0和a2=1.00.我们发现T1的误差非常大.导致这个问题的原因就是a11太小，用它做主元时会放大舍人误差.所以我们需要通过置换矩阵来选主元首先介绍置换矩阵的一些基本性质引理2.1设PERnxn为置换矩阵，XERnxn为任意矩阵，则(1PX相当于将X的行进行置换；XP相当于将X的列进行置换；(2)P-1=PT,即P是正交矩阵;(3) det(P) = ±1;(4）置换矩阵的乘积仍然是置换矩阵定理2.2（选主元LU分解的存在性）设AERnXn非奇异，则存在置换矩阵Pi，P2，以及单位下三角矩阵L和非奇异上三角矩阵U,使得PiAP2=LU.其中P和P中只有一个是必需的证明.用归纳法当n=1时.取P=P=L=1U=A即可假设结论对n-1成立。设AERnxn是n阶非奇异矩阵，则A至少存在一个非零元，取置换矩阵P1和P2使得a11A12PIAP,:A21A22其中 a11 ± 0, A22 E R(n-1)x(n-1), 记U12=A12,L21=A21/a11,U22=A22-L21U12.ui1 = a11,则有a11A12U12PAP2U22A2A21两边取行列式可得0 det(PIAP2) = dea11 -det(U22)

· 46 · 第二讲 线性方程组直接方法 即 A ≈ [ 1.00 0 1.72 × 102 1.00] [2.00 × 10−2 6.12 × 10 0 −1.05 × 104 ] . 解方程组 Ly = b 可得 y1 = 6.15 × 10, y2 = b2 − l21y1 ≈ −1.06 × 104 . 解方程组 Ux = y 可得 x2 = y2/u22 ≈ 1.01, x1 = (y1 − u12 ∗ x2)/u11 ≈ −0.413/u11 ≈ −20.7 □ 易知, 方程的精确解为 x1 = 10.0 和 x2 = 1.00. 我们发现 x1 的误差非常大. 导致这个问题的原因就 是 |a11| 太小, 用它做主元时会放大舍入误差. 所以我们需要通过置换矩阵来选主元. 首先介绍置换矩阵的一些基本性质. 引理 2.1 设 P ∈ R n×n 为置换矩阵, X ∈ R n×n 为任意矩阵, 则 (1) P X 相当于将 X 的行进行置换; XP 相当于将 X 的列进行置换; (2) P −1 = P ⊺ , 即 P 是正交矩阵; (3) det(P) = ±1; (4) 置换矩阵的乘积仍然是置换矩阵. 定理 2.2 (选主元 LU 分解的存在性) 设 A ∈ R n×n 非奇异, 则存在置换矩阵 P1, P2, 以及单位下三角 矩阵 L 和非奇异上三角矩阵 U, 使得 P1AP2 = LU. 其中 P1 和 P2 中只有一个是必需的. 证明. 用归纳法. 当 n = 1 时, 取 P1 = P2 = L = 1, U = A 即可. 假设结论对 n − 1 成立. 设 A ∈ R n×n 是 n 阶非奇异矩阵, 则 A 至少存在一个非零元, 取置换矩阵 Pˆ 1 和 Pˆ 2 使得 Pˆ 1APˆ 2 = [ a11 A12 A21 A22] , 其中 a11 ̸= 0, A22 ∈ R (n−1)×(n−1) . 记 u11 = a11, U12 = A12, L21 = A21/a11, U22 = A22 − L21U12. 则有 [ 1 0 L21 I ] [u11 U12 0 U22] = [ a11 A12 A21 A22] = Pˆ 1APˆ 2. 两边取行列式可得 0 ̸= det(P1APˆ 2) = det ([ 1 0 L21 I ]) · det ([u11 U12 0 U22]) = a11 · det(U22)