LU分解(LU Factorization)计算方法(手算+MATLAB),关于置换矩阵(Permutation Matrix),部分主元消去法(Partial Pivoting) |
您所在的位置:网站首页 › 矩阵可以做lu分解的条件 › LU分解(LU Factorization)计算方法(手算+MATLAB),关于置换矩阵(Permutation Matrix),部分主元消去法(Partial Pivoting) |
背景:
求解一些列具有相同系数矩阵的线性方程,如: A x = b 1 Ax=b_1 Ax=b1, A x = b 2 Ax=b_2 Ax=b2,… A x = b p Ax=b_p Ax=bp等。当矩阵 A A A可逆时,可以先求出矩阵 A A A的逆 A − 1 A^{-1} A−1,再计算 A − 1 b 1 A^{-1}b_1 A−1b1, A − 1 b 2 A^{-1}b_2 A−1b2等。 但是,也可以用LU分解法来解这一系列方程:先使用初等行变换化简解出 A x = b 1 Ax=b_1 Ax=b1,并同时得到矩阵 A A A的LU分解,剩下的方程使用LU分解法求解即可。 方法:设矩阵 A A A是 m × n m \times n m×n的矩阵,它可以使用行变换(不进行行的交换)化简为阶梯阵,则矩阵 A A A可以写成 A = L U A=LU A=LU的形式。 L L L是 m × m m\times m m×m的单位下三角矩阵(即主对角线元素全为1的下三角矩阵); U U U是 m × n m\times n m×n的上三角阶梯形矩阵。 用途当 A = L U A=LU A=LU时,方程 A x = b Ax=b Ax=b可以写成 L ( U x ) = b L(Ux)=b L(Ux)=b,令 U x = y Ux=y Ux=y,则有: L y = b U x = y Ly=b \qquad Ux=y Ly=bUx=y 即先求解 L y = b Ly=b Ly=b得到 y y y,再求解 U x = y Ux=y Ux=y得到 x x x,因为矩阵 L L L和 U U U都是三角矩阵,所以求解上述两个方程比直接求解 A x = b Ax=b Ax=b要简单。 例题(1)对矩阵 A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix} A=⎣⎡−510153−814−92⎦⎤进行LU分解。 解: A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix} A=⎣⎡−510153−814−92⎦⎤有3行,所以, L L L矩阵应该为 3 × 3 3\times 3 3×3的矩阵。 先求 U U U矩阵: 先把矩阵 A A A一步步消成上三角矩阵就能得到 U U U矩阵: A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] → [ − 5 3 4 0 − 2 − 1 0 10 14 ] → [ − 5 3 4 0 − 2 − 1 0 0 9 ] = U A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&10&14\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=U A=⎣⎡−510153−814−92⎦⎤→⎣⎡−5003−2104−114⎦⎤→⎣⎡−5003−204−19⎦⎤=U 再求 L L L矩阵: 把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵,然后每列除以此列的主元即得到 L L L矩阵: [ − 5 0 0 10 − 2 0 15 10 9 ] → [ 1 0 0 − 2 1 0 − 3 − 5 1 ] = L \begin{bmatrix}-5&0&0\\10&-2&0\\15&10&9\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}=L ⎣⎡−510150−210009⎦⎤→⎣⎡1−2−301−5001⎦⎤=L (-5, 10, 15这一列每个元素除以主元-5;-2,10这一列每个元素除以主元-2;9这一列每个元素除以主元9,即得到 L L L矩阵) 验算: L U = [ 1 0 0 − 2 1 0 − 3 − 5 1 ] [ − 5 3 4 0 − 2 − 1 0 0 9 ] = [ − 5 3 4 10 − 8 − 9 15 1 2 ] = A LU=\begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}\begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}=A LU=⎣⎡1−2−301−5001⎦⎤⎣⎡−5003−204−19⎦⎤=⎣⎡−510153−814−92⎦⎤=A 所以上述分解是正确的。 (2)对矩阵 A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix} A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤进行LU分解。 解: A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix} A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤有4行,所以, L L L矩阵应该为 4 × 4 4\times 4 4×4的矩阵, L L L矩阵的第一列应该是矩阵 A A A的第一列除以其主元元素(2)的结果: L = [ 1 0 0 0 − 2 1 0 0 1 1 0 − 3 1 ] L=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&\space&1&0\\-3&\space&\space&1\end{bmatrix} L=⎣⎢⎢⎡1−21−301 001 0001⎦⎥⎥⎤ 先求 U U U矩阵: 先把矩阵 A A A一步步消成上三角矩阵就能得到 U U U矩阵: A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 − 9 − 3 − 4 10 0 12 4 12 − 5 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 4 7 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 0 5 ] = U A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&-9&-3&-4&10\\0&12&4&12&-5\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&4&7\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=U A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤→⎣⎢⎢⎡200043−912−11−3452−412−2−310−5⎦⎥⎥⎤→⎣⎢⎢⎡20004300−11005224−2−317⎦⎥⎥⎤→⎣⎢⎢⎡20004300−11005220−2−315⎦⎥⎥⎤=U 再求 L L L矩阵: 把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵,然后每列除以此列的主元就得到 L L L矩阵: [ 2 0 0 0 − 4 3 0 0 2 − 9 2 0 − 6 12 4 5 ] → [ 1 0 0 0 − 2 1 0 0 1 − 3 1 0 − 3 4 2 1 ] = L \begin{bmatrix}2&0&0&0\\-4&3&0&0\\2&-9&2&0\\-6&12&4&5\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}=L ⎣⎢⎢⎡2−42−603−91200240005⎦⎥⎥⎤→⎣⎢⎢⎡1−21−301−3400120001⎦⎥⎥⎤=L (2,-,4, 2, -6这一列每个元素除以主元2;3,-9,12这一列每个元素除以主元3;2,4这一列每个元素除以主元2,5这一列除以主元5,即得到 L L L矩阵) 验算: L U = [ 1 0 0 0 − 2 1 0 0 1 − 3 1 0 − 3 4 2 1 ] [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 0 5 ] = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] = A LU=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}\begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}=A LU=⎣⎢⎢⎡1−21−301−3400120001⎦⎥⎥⎤⎣⎢⎢⎡20004300−11005220−2−315⎦⎥⎥⎤=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤=A 所以上述分解是正确的。 关于MATLAB中的LU分解函数MATLAB中提供了LU分解的函数lu( ),有两种用法,一种是返回三个矩阵(L,U和置换矩阵P),一种是只返回L和U矩阵。例如: F = 1 3 0 4 4 8 1 2 3 >> [L1 U1 P] = lu(F) //使用MATLAB内置的lu( )函数求矩阵F的LU分解,返回三个矩阵L,U和P的用法 L1 = 1.0000 0 0 0.2500 1.0000 0 0.2500 0.5000 1.0000 U1 = 4 4 8 0 2 -2 0 0 2 P = 0 1 0 1 0 0 0 0 1 >> [L2 U2] = lu(F) //使用MATLAB内置的lu( )函数求矩阵F的LU分解,返回2个矩阵L和U的用法 L2 = 0.2500 1.0000 0 1.0000 0 0 0.2500 0.5000 1.0000 U2 = 4 4 8 0 2 -2 0 0 2对于上述矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡141342083⎦⎤,手动LU分解为: L U = [ 1 0 0 4 1 0 1 1 8 1 ] [ 1 3 0 0 − 8 8 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F LU=\begin{bmatrix}1&0&0\\4&1&0\\1&\frac{1}{8}&1\end{bmatrix}\begin{bmatrix}1&3&0\\0&-8&8\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F LU=⎣⎡1410181001⎦⎤⎣⎡1003−80082⎦⎤=⎣⎡141342083⎦⎤=F 上面MATLAB第一种方法(L1,U1和P三个矩阵)算出来的结果: L 1 U 1 = [ 1 0 0 1 4 1 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 4 4 8 1 3 0 1 2 3 ] ≠ F L_1U_1=\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix}\neq F L1U1=⎣⎡141410121001⎦⎤⎣⎡4004208−22⎦⎤=⎣⎡411432803⎦⎤=F MATLAB第二种方法(L2和U2 两个矩阵)算出来的结果: L 2 U 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F L_2U_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F L2U2=⎣⎡411411021001⎦⎤⎣⎡4004208−22⎦⎤=⎣⎡141342083⎦⎤=F 可见输出两个矩阵的方法得到LU矩阵是正确的,第一种方法得到的 L 1 U 1 ≠ F L_1U_1\neq F L1U1=F。两种方法得到的LU矩阵都和手动算的结果不一样。 原因: 第一种方法:输入[L1 U1 P] = lu(F)命令的时候,MATLAB首先生成 P F = [ 4 4 8 1 3 0 1 2 3 ] PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} PF=⎣⎡411432803⎦⎤,其中 P P P为置换矩阵(Permutation Matrix), P = [ 0 1 0 1 0 0 0 0 1 ] P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix} P=⎣⎡010100001⎦⎤,所以,这种方法求得的分解是 L 1 U 1 = P F L_1U_1=PF L1U1=PF而不是 L 1 U 1 = F L_1U_1=F L1U1=F。 第二种方法:输入[L2 U2] = lu(F)命令的时候,算出来的 U 2 U_2 U2矩阵和第一种方法一样(但是和手算的不一样);算出来的 L 2 L_2 L2矩阵其实是 L 2 = P T L 1 L_2=P^TL_1 L2=PTL1,即: L 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] = P T L 1 = [ 0 1 0 1 0 0 0 0 1 ] [ 1 0 0 1 4 1 0 1 4 1 2 1 ] L_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}=P^TL_1=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix} L2=⎣⎡411411021001⎦⎤=PTL1=⎣⎡010100001⎦⎤⎣⎡141410121001⎦⎤ 第二种方法算出来的 L 2 U 2 = F L_2U_2=F L2U2=F是成立的。 关于置换矩阵(Permutation Matrix)上MATLAB计算LU分解时,两种方法都用到了置换矩阵 P P P。 置换矩阵是(0,1)矩阵,即只有0和1两种元素的矩阵,其作用是让矩阵的行进交换的矩阵。置换矩阵 P P P是改变相应阶数的单位矩阵得到的,例如对单位矩阵 I I I交换第一第二行得到 P 21 = [ 0 1 0 1 0 0 0 0 1 ] P_{21}=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix} P21=⎣⎡010100001⎦⎤,那么得到的置换矩阵 P 21 P_{21} P21的功能就是,使得和它相乘的矩阵也交换第一第二行。 例如: P 21 = [ 0 1 0 1 0 0 0 0 1 ] P_{21}=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix} P21=⎣⎡010100001⎦⎤对于上述矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡141342083⎦⎤,让 P 21 P_{21} P21乘以 F F F就是把 F F F的第二行和第一行交换: P 21 F = [ 0 1 0 1 0 0 0 0 1 ] [ 1 3 0 4 4 8 1 2 3 ] = [ 4 4 8 1 3 0 1 2 3 ] P_{21}F=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} P21F=⎣⎡010100001⎦⎤⎣⎡141342083⎦⎤=⎣⎡411432803⎦⎤ 类似的,还有: P 32 = [ 1 0 0 0 0 1 0 1 0 ] P_{32}=\begin{bmatrix}1&0&0\\0&0&1\\0&1&0\end{bmatrix} P32=⎣⎡100001010⎦⎤ P 31 = [ 0 0 1 0 1 0 1 0 0 ] P_{31}=\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\end{bmatrix} P31=⎣⎡001010100⎦⎤ 关于部分主元消去法(Partial Pivoting)部分主元消去法在求解线性系统时常被用到(如LU分解和使用A\b求 A x = b Ax=b Ax=b等),在计算机的浮点数计算中能提高精确度。 例如,上述第一种方法求LU分解时,对于矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡141342083⎦⎤,MATLAB首先生成 P F = [ 4 4 8 1 3 0 1 2 3 ] PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} PF=⎣⎡411432803⎦⎤,把 F F F的第一第二行进行交换,让第一列最大的元素4放在主元位置。 类似的,使用部分主元消去法时,每一行主元位置的元素与其下面的元素中最大的元素放在主元位置(使用行交换将最大的元素换到主元位置),后面使用高斯消元法,这就是部分主元消去法。所以,上面求解中MATLAB使用部分主元消去法时,先使用置换矩阵 P P P来调整主元使之最大,然后再求解。 部分主元消去法就是多了一个把最大的值调到主元位置的操作,却能提高计算精度。 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |