数值分析--python--LU分解法

**

1、LU分解法介绍

**
LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法（Doolittle algorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方) 左右。

2、示例代码

    import numpy as np
    
    a = [[2, 2, 3], [4, 7, 7], [-2, 4, 5]]
    b = [[3], [1], [-7]]
    unit_matrix=[[1,0,0],[0,1,0],[0,0,1]]
    
     #合并两个矩阵
    def augmented_matrix(x, y):
        return [aa + bb for aa, bb in zip(x, y)]
    
    def LU_decomposition():
        a_unit_matrix=augmented_matrix(a,unit_matrix)
        print(np.array(a_unit_matrix))
        # 交换主元
        for j in range(len(a_unit_matrix)):
            if a_unit_matrix[j][j] == 0:
                for w in range(len(a_unit_matrix)):
                    if a_unit_matrix[w][j] != 0:
                        l = a_unit_matrix[w]
                        a_unit_matrix[w] = a_unit_matrix[j]
                        a_unit_matrix[j] = l
                    break
                    
        #化为上三角矩阵
        j = 0
        while j < len(a_unit_matrix):
            line = a_unit_matrix[j]
            lest_0 = []
            if j==0:
                for x0 in line:
                    x = x0
                    lest_0.append(x)
            else:
                for x0 in line:
                    x = x0/a_unit_matrix[0][0]
                    lest_0.append(x)
                a_unit_matrix[j] = lest_0
            flag = j + 1
            while flag < len(a_unit_matrix):
                lest_1=[]
                temp1 = a_unit_matrix[flag][j]/a_unit_matrix[j][j]
                i = 0
                for x1 in a_unit_matrix[flag]:
                    x = x1-(temp1 * lest_0[i])
                    lest_1.append(x)
                    i=i+1
                a_unit_matrix[flag] = lest_1
                flag = flag + 1
            j = j + 1
        print("将A矩阵与单位矩阵合并，然后化为下三角矩阵为：")
    
    
        for i in range(len(a_unit_matrix)):
            m = 1 / a_unit_matrix[i][i+len(a)]
            for j in range(len(a_unit_matrix[0])):
                a_unit_matrix[i][j]=m*a_unit_matrix[i][j]
        print(np.array(a_unit_matrix))
        # 分解为L的逆矩阵和U矩阵
        L=[]
        U=[]
        for i in range(len(a_unit_matrix)):
            L.append([])
            U.append([])
        for i in range(len(a_unit_matrix)):
            for j in range(len(a_unit_matrix[i])):
                if j<len(a_unit_matrix[i])/2:
                    U[i].append(a_unit_matrix[i][j])
                else:
                    L[i].append(a_unit_matrix[i][j])
    
        print("L矩阵为：")
        print(np.linalg.inv(np.array(L)))
        print("U矩阵为：")
        print(np.array(U))
    
        # 求D矩阵
        D=np.dot(L,b)
        print("D矩阵为：")
        print(D)
    
        # 通过U的逆矩阵求X矩阵
        U_inverse = np.linalg.inv(np.array(U))
        X = np.dot(U_inverse, D)
        print("X矩阵为：")
        print(X)
    
    LU_decomposition()