线性代数 - 线性方程组的 LU 分解解法
线性代数 - 线性方程组的 LU 分解解法
flyfish
线性代数 - 线性方程组的原始解法(高斯消元法)
线性代数 - LU分解(LU-Factorization、LU Decomposition)
线性方程组
方程组:{2x1+x2+x3=1(1)4x1+3x2+3x3=2(2)8x1+7x2+9x3=6(3)即 Ax=b,A=[211433879],b=[126]\text{方程组:}\begin{cases} 2x_1 + x_2 + x_3 = 1 \quad (1) \\ 4x_1 + 3x_2 + 3x_3 = 2 \quad (2) \\ 8x_1 + 7x_2 + 9x_3 = 6 \quad (3) \end{cases} \quad \text{即}\ \mathbf{Ax}=\mathbf{b},\ \mathbf{A}=\begin{bmatrix}2&1&1\\4&3&3\\8&7&9\end{bmatrix},\ \mathbf{b}=\begin{bmatrix}1\\2\\6\end{bmatrix} 方程组:⎩⎨⎧2x1+x2+x3=1(1)4x1+3x2+3x3=2(2)8x1+7x2+9x3=6(3)即 Ax=b, A=248137139, b=126
LU分解
原始矩阵 A=[211433879]\mathbf{A} = \begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix}A=248137139
一、LU分解目标
将方阵A\mathbf{A}A分解为两个特殊矩阵的乘积:
L\mathbf{L}L:单位下三角矩阵(主对角线元素全为1,主对角线以上元素全为0,仅下方有非零元素);
U\mathbf{U}U:上三角矩阵(主对角线以下元素全为0,仅主对角线及上方有非零元素);
最终满足 A=LU\mathbf{A} = \mathbf{L}\mathbf{U}A=LU。
二、步骤1:初始化L\mathbf{L}L和U\mathbf{U}U
初始化L\mathbf{L}L:因A\mathbf{A}A是3阶矩阵,L\mathbf{L}L为3阶单位下三角矩阵(对角线填1,其余位置先填0):
L=[100010001]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} L=100010001
初始化U\mathbf{U}U:将A\mathbf{A}A的元素直接复制到U\mathbf{U}U中(后续通过消元逐步修改U\mathbf{U}U,使其成为上三角矩阵):
U=[211433879]\mathbf{U} = \begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix} U=248137139
三、步骤2:第一阶段消元——消去U\mathbf{U}U第2、3行的“第1列元素”(消去x1x_1x1的系数)
目标:让U\mathbf{U}U第2行第1列、第3行第1列的元素变为0(主对角线下方第1列全为0),同时将“消元用的倍数”记录到L\mathbf{L}L的对应位置。
1. 处理U\mathbf{U}U的第2行(消去第2行第1列元素)
- 确定主元:U\mathbf{U}U第1行第1列的元素 u11=2u_{11} = 2u11=2(主元是消元的“基准”,需非零)。
- 计算消元倍数kkk:要让U\mathbf{U}U第2行第1列的元素(当前为4)变为0,需用“第2行 = 第2行 - k×k×k×第1行”,其中:
k=U第2行第1列元素主元u11=42=2k = \frac{\text{U第2行第1列元素}}{\text{主元}u_{11}} = \frac{4}{2} = 2 k=主元u11U第2行第1列元素=24=2 - 记录倍数到L\mathbf{L}L:这个k=2k=2k=2是“用第1行消第2行”的倍数,对应L\mathbf{L}L的“第2行第1列”(因L\mathbf{L}L的lijl_{ij}lij表示“用第j行消第i行的倍数”,此处i=2,j=1i=2,j=1i=2,j=1),因此更新L\mathbf{L}L:
L=[100210001]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ \boxed{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} L=120010001 - 更新U\mathbf{U}U第2行:按“第2行 = 第2行 - 2×第1行”逐元素计算:
- 第2行第1列:4−2×2=04 - 2×2 = 04−2×2=0(目标达成,变为0);
- 第2行第2列:3−2×1=13 - 2×1 = 13−2×1=1;
- 第2行第3列:3−2×1=13 - 2×1 = 13−2×1=1;
因此,U\mathbf{U}U更新为:
U=[211011879]\mathbf{U} = \begin{bmatrix} 2 & 1 & 1 \\ \boxed{0} & 1 & 1 \\ 8 & 7 & 9 \end{bmatrix} U=208117119
2. 处理U\mathbf{U}U的第3行(消去第3行第1列元素)
- 仍用主元u11=2u_{11}=2u11=2:要让U\mathbf{U}U第3行第1列的元素(当前为8)变为0。
- 计算消元倍数kkk:
k=U第3行第1列元素主元u11=82=4k = \frac{\text{U第3行第1列元素}}{\text{主元}u_{11}} = \frac{8}{2} = 4 k=主元u11U第3行第1列元素=28=4 - 记录倍数到L\mathbf{L}L:此k=4k=4k=4是“用第1行消第3行”的倍数,对应L\mathbf{L}L的“第3行第1列”(i=3,j=1i=3,j=1i=3,j=1),更新L\mathbf{L}L:
L=[100210401]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ \boxed{4} & 0 & 1 \end{bmatrix} L=124010001 - 更新U\mathbf{U}U第3行:按“第3行 = 第3行 - 4×第1行”逐元素计算:
- 第3行第1列:8−4×2=08 - 4×2 = 08−4×2=0(目标达成,变为0);
- 第3行第2列:7−4×1=37 - 4×1 = 37−4×1=3;
- 第3行第3列:9−4×1=59 - 4×1 = 59−4×1=5;
因此,U\mathbf{U}U更新为:
U=[211011035]\mathbf{U} = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ \boxed{0} & 3 & 5 \end{bmatrix} U=200113115
第一阶段消元后状态:
L=[100210401]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 0 & 1 \end{bmatrix}L=124010001,U\mathbf{U}U第1列下方全为0(已完成第1列消元)。
四、步骤3:第二阶段消元——消去U\mathbf{U}U第3行的“第2列元素”(消去x2x_2x2的系数)
目标:让U\mathbf{U}U第3行第2列的元素变为0(主对角线下方第2列全为0),此时U\mathbf{U}U将成为上三角矩阵。
1. 处理U\mathbf{U}U的第3行
- 确定新主元:U\mathbf{U}U第2行第2列的元素 u22=1u_{22} = 1u22=1(因第1列已消完,当前消第2列,主元为第2行第2列)。
- 计算消元倍数kkk:要让U\mathbf{U}U第3行第2列的元素(当前为3)变为0,需用“第3行 = 第3行 - k×k×k×第2行”,其中:
k=U第3行第2列元素主元u22=31=3k = \frac{\text{U第3行第2列元素}}{\text{主元}u_{22}} = \frac{3}{1} = 3 k=主元u22U第3行第2列元素=13=3 - 记录倍数到L\mathbf{L}L:此k=3k=3k=3是“用第2行消第3行”的倍数,对应L\mathbf{L}L的“第3行第2列”(i=3,j=2i=3,j=2i=3,j=2),更新L\mathbf{L}L:
L=[100210431]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & \boxed{3} & 1 \end{bmatrix} L=124013001 - 更新U\mathbf{U}U第3行:按“第3行 = 第3行 - 3×第2行”逐元素计算:
- 第3行第2列:3−3×1=03 - 3×1 = 03−3×1=0(目标达成,变为0);
- 第3行第3列:5−3×1=25 - 3×1 = 25−3×1=2;
(第3行第1列已为0,无需计算)
因此,U\mathbf{U}U最终更新为上三角矩阵:
U=[211011002]\mathbf{U} = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & \boxed{0} & 2 \end{bmatrix} U=200110112
第二阶段消元后状态:
L\mathbf{L}L是完整的单位下三角矩阵,U\mathbf{U}U是完整的上三角矩阵,LU分解初步完成。
五、步骤4:验证分解正确性——证明A=LU\mathbf{A} = \mathbf{L}\mathbf{U}A=LU
通过矩阵乘法计算LU\mathbf{L}\mathbf{U}LU,看结果是否等于原始矩阵A\mathbf{A}A。
矩阵乘法规则:L\mathbf{L}L的第iii行 × U\mathbf{U}U的第jjj列 = 结果矩阵的第iii行第jjj列元素(对应元素相乘再相加)。
1. 计算LU\mathbf{L}\mathbf{U}LU的第1行
L\mathbf{L}L第1行:[1,0,0][1, 0, 0][1,0,0] × U\mathbf{U}U各列:
- 第1行第1列:1×2+0×0+0×0=21×2 + 0×0 + 0×0 = 21×2+0×0+0×0=2;
- 第1行第2列:1×1+0×1+0×0=11×1 + 0×1 + 0×0 = 11×1+0×1+0×0=1;
- 第1行第3列:1×1+0×1+0×2=11×1 + 0×1 + 0×2 = 11×1+0×1+0×2=1;
结果第1行:[2,1,1][2, 1, 1][2,1,1](与A\mathbf{A}A第1行一致)。
2. 计算LU\mathbf{L}\mathbf{U}LU的第2行
L\mathbf{L}L第2行:[2,1,0][2, 1, 0][2,1,0] × U\mathbf{U}U各列:
- 第2行第1列:2×2+1×0+0×0=42×2 + 1×0 + 0×0 = 42×2+1×0+0×0=4;
- 第2行第2列:2×1+1×1+0×0=32×1 + 1×1 + 0×0 = 32×1+1×1+0×0=3;
- 第2行第3列:2×1+1×1+0×2=32×1 + 1×1 + 0×2 = 32×1+1×1+0×2=3;
结果第2行:[4,3,3][4, 3, 3][4,3,3](与A\mathbf{A}A第2行一致)。
3. 计算LU\mathbf{L}\mathbf{U}LU的第3行
L\mathbf{L}L第3行:[4,3,1][4, 3, 1][4,3,1] × U\mathbf{U}U各列:
- 第3行第1列:4×2+3×0+1×0=84×2 + 3×0 + 1×0 = 84×2+3×0+1×0=8;
- 第3行第2列:4×1+3×1+1×0=74×1 + 3×1 + 1×0 = 74×1+3×1+1×0=7;
- 第3行第3列:4×1+3×1+1×2=4+3+2=94×1 + 3×1 + 1×2 = 4 + 3 + 2 = 94×1+3×1+1×2=4+3+2=9;
结果第3行:[8,7,9][8, 7, 9][8,7,9](与A\mathbf{A}A第3行一致)。
六、最终分解结果
A=LU=[100210431]⏟L(单位下三角矩阵)×[211011002]⏟U(上三角矩阵)\mathbf{A} = \mathbf{L}\mathbf{U} = \underbrace{\begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{bmatrix}}_{\mathbf{L}(单位下三角矩阵)} \times \underbrace{\begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{bmatrix}}_{\mathbf{U}(上三角矩阵)} A=LU=L(单位下三角矩阵)124013001×U(上三角矩阵)200110112
将系数矩阵 A=[211433879]\mathbf{A} = \begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix}A=248137139 分解为 L=[100210431]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{bmatrix}L=124013001(单位下三角矩阵)和 U=[211011002]\mathbf{U} = \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{bmatrix}U=200110112(上三角矩阵)。
应用于方程组
接下来,我们用这两个矩阵求解原始线性方程组 Ax=b\mathbf{Ax} = \mathbf{b}Ax=b(其中 b=[126]\mathbf{b} = \begin{bmatrix} 1 \\ 2 \\ 6 \end{bmatrix}b=126,方程组为 {2x1+x2+x3=14x1+3x2+3x3=28x1+7x2+9x3=6\begin{cases} 2x_1 + x_2 + x_3 = 1 \\ 4x_1 + 3x_2 + 3x_3 = 2 \\ 8x_1 + 7x_2 + 9x_3 = 6 \end{cases}⎩⎨⎧2x1+x2+x3=14x1+3x2+3x3=28x1+7x2+9x3=6),分两步解三角方程组:先解 Ly=b\mathbf{Ly} = \mathbf{b}Ly=b,再解 Ux=y\mathbf{Ux} = \mathbf{y}Ux=y。
第一步:解下三角方程组 Ly=b\mathbf{Ly} = \mathbf{b}Ly=b(前向替换)
下三角矩阵的特点是“主对角线以下有元素,以上全为0”,因此从第一个未知数开始,逐个代入计算(前向替换)。
设 y=[y1y2y3]\mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}y=y1y2y3,代入 L\mathbf{L}L 和 b\mathbf{b}b 得方程组:
{1⋅y1+0⋅y2+0⋅y3=1(1)2⋅y1+1⋅y2+0⋅y3=2(2)4⋅y1+3⋅y2+1⋅y3=6(3)\begin{cases} 1 \cdot y_1 + 0 \cdot y_2 + 0 \cdot y_3 = 1 \quad (1) \\ 2 \cdot y_1 + 1 \cdot y_2 + 0 \cdot y_3 = 2 \quad (2) \\ 4 \cdot y_1 + 3 \cdot y_2 + 1 \cdot y_3 = 6 \quad (3) \end{cases} ⎩⎨⎧1⋅y1+0⋅y2+0⋅y3=1(1)2⋅y1+1⋅y2+0⋅y3=2(2)4⋅y1+3⋅y2+1⋅y3=6(3)
计算过程:
- 由方程(1):y1=1y_1 = 1y1=1(直接得出,因为其他项系数为0);
- 代入方程(2):2×1+y2=2⟹y2=2−2=02 \times 1 + y_2 = 2 \implies y_2 = 2 - 2 = 02×1+y2=2⟹y2=2−2=0;
- 代入方程(3):4×1+3×0+y3=6⟹y3=6−4=24 \times 1 + 3 \times 0 + y_3 = 6 \implies y_3 = 6 - 4 = 24×1+3×0+y3=6⟹y3=6−4=2。
最终得到 y=[102]\mathbf{y} = \begin{bmatrix} 1 \\ 0 \\ 2 \end{bmatrix}y=102。
第二步:解上三角方程组 Ux=y\mathbf{Ux} = \mathbf{y}Ux=y(后向替换)
上三角矩阵的特点是“主对角线以上有元素,以下全为0”,因此从最后一个未知数开始,反向代入计算(后向替换)。
设 x=[x1x2x3]\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}x=x1x2x3,代入 U\mathbf{U}U 和 y\mathbf{y}y 得方程组:
{2⋅x1+1⋅x2+1⋅x3=1(4)0⋅x1+1⋅x2+1⋅x3=0(5)0⋅x1+0⋅x2+2⋅x3=2(6)\begin{cases} 2 \cdot x_1 + 1 \cdot x_2 + 1 \cdot x_3 = 1 \quad (4) \\ 0 \cdot x_1 + 1 \cdot x_2 + 1 \cdot x_3 = 0 \quad (5) \\ 0 \cdot x_1 + 0 \cdot x_2 + 2 \cdot x_3 = 2 \quad (6) \end{cases} ⎩⎨⎧2⋅x1+1⋅x2+1⋅x3=1(4)0⋅x1+1⋅x2+1⋅x3=0(5)0⋅x1+0⋅x2+2⋅x3=2(6)
计算过程:
- 由方程(6):2x3=2⟹x3=12x_3 = 2 \implies x_3 = 12x3=2⟹x3=1(直接得出,其他项系数为0);
- 代入方程(5):x2+1=0⟹x2=0−1=−1x_2 + 1 = 0 \implies x_2 = 0 - 1 = -1x2+1=0⟹x2=0−1=−1;
- 代入方程(4):2x1+(−1)+1=1⟹2x1=1⟹x1=0.52x_1 + (-1) + 1 = 1 \implies 2x_1 = 1 \implies x_1 = 0.52x1+(−1)+1=1⟹2x1=1⟹x1=0.5(或 12\frac{1}{2}21)。
最终解与验证
方程组 Ax=b\mathbf{Ax} = \mathbf{b}Ax=b 的解为:
x=[x1x2x3]=[0.5−11]\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 0.5 \\ -1 \\ 1 \end{bmatrix} x=x1x2x3=0.5−11
验证(代入原方程组,确保正确):
第一个方程:2×0.5+(−1)+1=1−1+1=12 \times 0.5 + (-1) + 1 = 1 - 1 + 1 = 12×0.5+(−1)+1=1−1+1=1(与右边相等);
第二个方程:4×0.5+3×(−1)+3×1=2−3+3=24 \times 0.5 + 3 \times (-1) + 3 \times 1 = 2 - 3 + 3 = 24×0.5+3×(−1)+3×1=2−3+3=2(与右边相等);
第三个方程:8×0.5+7×(−1)+9×1=4−7+9=68 \times 0.5 + 7 \times (-1) + 9 \times 1 = 4 - 7 + 9 = 68×0.5+7×(−1)+9×1=4−7+9=6(与右边相等)。
