第二讲线性方程组直接方法: 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)