LU分解求线性方程组的解

       LU分解是矩阵分解的一种，可以将一个矩阵分解为一个上三角矩阵和一个下三角矩阵的乘积。

LU分解可以用来求逆矩阵，解线性方程组等。本文将介绍LU分解求线性方程组的解。

    1.定义

        如果A是一个方阵，A的LU分解是将A分解成如下形式:

                         

其中L,U分别为下三角矩阵和上三角矩阵。

    2.例子

      对于如下矩阵A，对A进行LU分解

                                       

       首先将矩阵第一列对角线上元素A11下面的元素通过矩阵初等行变换变为0，

                                    

    

       然后再将矩阵第二列对角线上元素A22 下面的元素通过矩阵初等行变换变为0。                                            

                                                                                

则得到的上三角矩阵就是U。这个时候，L也已经求出来了。通过将下三角形主对角线上的元素

都置为1，乘数因子放在下三角相应的位置（放在消元时将元素变为0的那个元素的位置），就

可以得到下三角矩阵L。如下:

                                                         

 

对于L的构造，举个例子。如将第一列的元素2变为0时，第二行减去第一行乘以2，于是A21

就变成了0。这个乘数因子将元素A21变成了0，对应的，下三角矩阵L中对应位置的元素L21就为

乘数因子2。其它的与之类似。

       3.LU分解程序实现(java实现)

            通过上面举的例子，我们应该对LU分解的过程有了一个大致的了解。接下来可以看看程序

是怎么实现LU分解的，进一步加深对LU分解的了解。

[java] view plain copy

1. /** 
2.      * Get matrix L and U. list.get(0) for L, list.get(1) for U 
3.      * @param a - Coefficient matrix of the equations 
4.      * @return matrix L and U, list.get(0) for L, list.get(1) for U 
5.      */  
6.     private static List<double[][]> decomposition(double[][] a) {  
7.         final double esp = 0.000001;          
8.         double[][] U = a;  
9.         double[][] L = createIdentiyMatrix(a.length);  
10.           
11.         for(int j=0; j<a[0].length - 1; j++) {             
12.             if(Math.abs(a[j][j]) < esp) {  
13.                 throw new IllegalArgumentException("zero pivot encountered.");  
14.             }  
15.               
16.             for(int i=j+1; i<a.length; i++) {  
17.                 double mult = a[i][j] / a[j][j];   
18.                 for(int k=j; k<a[i].length; k++) {  
19.                     U[i][k] = a[i][k] - a[j][k] * mult;  
20.                 }  
21.                 L[i][j] = mult;  
22.             }  
23.         }  
24.         return Arrays.asList(L, U);  
25.     }  

        上面函数中出现的 createIdentiyMatrix 函数就是创建一个 a.length()行a.length()列的单位矩阵。

而

           Math.abs(a[j][i]) < esp

的作用是当主元为0时，经典的LU分解算法不再适用，后面将进一步谈论这个。

 

      4.LU分解解线性方程组

           通过上面的介绍，我们已经知道，一个方阵A可以分解成 A=LU的形式(这里假设矩阵A能够进行LU分解)。

对于一个线性方程组 Ax=b，则由 A=LU 有 LUx = b。

           为了求出x,我们可以先将Ux看成一个整体V(V=UX)，通过求解线性方程组 LV=b 得到V,即Ux,

再通过求解线性方程组 Ux=V 即可求出 x。

           看到这里，你可能会觉得这样求解很麻烦。但是，别忘了L和U分别是下三角矩阵和上三角矩阵，

求解释很容易，不需要通过高斯消去法等求线性方程组的算法来求解。

           首先来看一下 LV=b 求解V的程序代码：        

[java] view plain copy

1. /** 
2.      * Get U multiply X 
3.      * @param a - Coefficient matrix of the equations 
4.      * @param b - right-hand side of the equations 
5.      * @param L - L of LU Decomposition 
6.      * @return U multiply X 
7.      */  
8.     private static double[] getUMultiX(double[][] a, double[] b, double[][] L) {  
9.         double[] UMultiX = new double[a.length];  
10.         for(int i=0; i<a.length; i++) {  
11.             double right_hand = b[i];  
12.             for(int j=0; j<i; j++) {  
13.                 right_hand -= L[i][j] * UMultiX[j];  
14.             }  
15.             UMultiX[i] = right_hand / L[i][i];  
16.         }  
17.         return UMultiX;  
18.     }  

         然后是 Ux=V 求解的程序代码:

[java] view plain copy

1. /** 
2.      * Get solution of the equations 
3.      * @param a - Coefficient matrix of the equations 
4.      * @param U - U of LU Decomposition 
5.      * @param UMultiX - U multiply X 
6.      * @return Equations solution 
7.      */  
8.     private static double[] getSolution(double[][] a, double[][] U,  
9.             double[] UMultiX) {  
10.         double[] solutions = new double[a[0].length];  
11.         for(int i=U.length-1; i>=0; i--) {  
12.             double right_hand = UMultiX[i];  
13.             for(int j=U.length-1; j>i; j--) {  
14.                 right_hand -= U[i][j] * solutions[j];  
15.             }  
16.             solutions[i] = right_hand / U[i][i];  
17.         }  
18.         return solutions;  
19.     }  

       这两个求解过程 是不是很简单 ？    

       如果觉得整个LU分解求解方程组的解过程 还没有连接起来的话，可以看看下面整个程序的完整代码。

[java] view plain copy

1. import java.util.Arrays;  
2. import java.util.List;  
3.   
4. public class LUDecomposition {  
5.     /** 
6.      * Get solutions of the equations 
7.      * @param a - Coefficient matrix of the equations 
8.      * @param b - right-hand side of the equations 
9.      * @return solution of the equations 
10.      */  
11.     public static double[] solve(double[][] a, double[] b) {          
12.         List<double[][]> LAndU = decomposition(a);  //LU decomposition  
13.         double[][] L = LAndU.get(0);  
14.         double[][] U = LAndU.get(1);  
15.         double[] UMultiX = getUMultiX(a, b, L);       
16.         return getSolution(a, U, UMultiX);        
17.     }  
18.   
19.     /** 
20.      * Get solution of the equations 
21.      * @param a - Coefficient matrix of the equations 
22.      * @param U - U of LU Decomposition 
23.      * @param UMultiX - U multiply X 
24.      * @return Equations solution 
25.      */  
26.     private static double[] getSolution(double[][] a, double[][] U,  
27.             double[] UMultiX) {  
28.         double[] solutions = new double[a[0].length];  
29.         for(int i=U.length-1; i>=0; i--) {  
30.             double right_hand = UMultiX[i];  
31.             for(int j=U.length-1; j>i; j--) {  
32.                 right_hand -= U[i][j] * solutions[j];  
33.             }  
34.             solutions[i] = right_hand / U[i][i];  
35.         }  
36.         return solutions;  
37.     }  
38.   
39.     /** 
40.      * Get U multiply X 
41.      * @param a - Coefficient matrix of the equations 
42.      * @param b - right-hand side of the equations 
43.      * @param L - L of LU Decomposition 
44.      * @return U multiply X 
45.      */  
46.     private static double[] getUMultiX(double[][] a, double[] b, double[][] L) {  
47.         double[] UMultiX = new double[a.length];  
48.         for(int i=0; i<a.length; i++) {  
49.             double right_hand = b[i];  
50.             for(int j=0; j<i; j++) {  
51.                 right_hand -= L[i][j] * UMultiX[j];  
52.             }  
53.             UMultiX[i] = right_hand / L[i][i];  
54.         }  
55.         return UMultiX;  
56.     }  
57.       
58.     /** 
59.      * Get matrix L and U. list.get(0) for L, list.get(1) for U 
60.      * @param a - Coefficient matrix of the equations 
61.      * @return matrix L and U, list.get(0) for L, list.get(1) for U 
62.      */  
63.     private static List<double[][]> decomposition(double[][] a) {  
64.         double[][] U = a;  
65.         double[][] L = createIdentityMatrix(a.length);  
66.           
67.         for(int j=0; j<a[0].length - 1; j++) {             
68.             if(a[j][j] == 0) {  
69.                  throw new IllegalArgumentException("zero pivot encountered.");  
70.              }  
71.               
72.             for(int i=j+1; i<a.length; i++) {  
73.                  double mult = a[i][j] / a[j][j];   
74.                 for(int k=j; k<a[i].length; k++) {  
75.                      U[i][k] = a[i][k] - a[j][k] * mult;  
76.                  }  
77.                 L[i][j] = mult;  
78.             }  
79.         }  
80.         return Arrays.asList(L, U);  
81.     }  
82.   
83.     private static double[][] createIdentityMatrix(int row) {  
84.         double[][] identityMatrix = new double[row][row];  
85.         for(int i=0; i<identityMatrix.length; i++) {  
86.             for(int j=i; j<identityMatrix[i].length; j++) {  
87.                 if(j == i) {   
88.                     identityMatrix[i][j] = 1;  
89.                 } else {  
90.                     identityMatrix[i][j] = 0;  
91.                 }  
92.             }  
93.         }  
94.         return identityMatrix;  
95.     }  
96. }               

      5. LU分解的不足及改进

         经典的LU分解算法当方阵中主元(主对角线上的元素) 出现0时，上面介绍的经典LU分解算法将失效，

上面的算法中也已经体现出来了。不过，我们可以在A=LU分解的基础上做出比较小的改动，就可以

使这个算法在上述情况下来能适用。随人 PA=LU 方法也可以解决这一问题，但是计算的耗费较大，

我们可以 将  decomposition 函数中的 对主元是否为0进行判断的 if 语句

                               if(a[j][j] == 0) {......}

         改为

                               if(a[j][j] == 0) { a[j][j] = 1e-20; }

            这样就可以方便的解决主元为0的问题，而且不需要额外的计算。

       6.参考文献

        1.  维基百科，http://zh./zh-cn/LU%E5%88%86%E8%A7%A3

        2.  Numerical Analysis, 2nd edition,  Timothy Sauer.