数值分析--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()