Python - 求解具有LU分解的线性方程组的数值问题

Question

信息

操作系统：曼扎罗 Linux

开发环境：Spyder4(Python 3.7(

库：Numpy

问题

你好

我编写了一些函数，按照以下三种方法求解线性方程组：

1. LU分解
2. 雅可比法
3. 高斯-赛德尔法

该程序运行完美。然而，LU分解的结果困扰着我。 例如，如果我的矩阵 A 和向量 b 是

A = np.array([[10., -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]])
b = np.array([6., 25., -11., 15.])

那么我的结果是

1. LU：[ 0.95034091 1.67840909 -0.9125 1.875 ]
2. 雅各比：[ 1. 2. -1. 1.]
3. 高斯-塞德尔：[ 1. 2. -1. 1.]

如您所见，LU 给我的结果略有不同。是舍入错误还是截断错误的问题？任何建议或帮助将不胜感激。

干杯！

<小时 />

编辑1：

我修复了这个问题，它只是 LU 分解过程中代码中的错误。感谢大家的反馈。

共享代码如下：

def lu_decomposition(matrix_in, b_in, n):
lower = [[0 for x in range(n)]
for y in range(n)];
upper = [[0 for x in range(n)]
for y in range(n)];
# Doolittle's Method
# Decomposing matrix into upper and lower matrices
for i in range(n):
# Upper triangle
for k in range(i, n):
# Sigma from j to i for each row k of U
sum = 0
for j in range(i):
sum += lower[i][j] * upper[j][k]
# Evaluate U for row k
upper[i][k] = matrix_in[i][k] - sum

# Lower Triangle
for k in range(i,n):
if(i == k): # Entry of a diagonal element
lower[i][i] = 1 
else:
# Sigma from j to i for each column k of L
sum = 0
for j in range(i):
sum += (lower[k][j] * upper[j][i])
# Evaluate L for column k
lower[k][i] = int( (matrix_in[k][i] - sum)/ upper[i][i])

# Perform forward substitution Ly=b
y = [0 for x in range(n)]
lower_inv = np.linalg.inv(lower)
y = np.dot(lower_inv, b_in)
# Perform back substitution Ux=y
x_sol = [0 for x in range(n)]
upper_inv = np.linalg.inv(upper)
x_sol = np.dot(upper_inv, y)
# printing results
# setw is for displaying nicely 
print("Lower TriangularttUpper Triangular"); 
# Displaying the result : 
for i in range(n): 

# Lower 
for j in range(n): 
print(lower[i][j], end = "t");  
print("", end = "t"); 
# Upper 
for j in range(n): 
print(upper[i][j], end = "t"); 
print(""); 

print("Here's the solution for your matrix: ")
print(x_sol)

Answer 1

代码中有很多问题需要解决：

• Python不需要;要结束一行，python知道行何时结束。蟒蛇只使用;将两行代码放在同一物理行上。比如"x=7;打印(x("。
• 如果在列表推导式中不使用变量，则习惯上使用 _。例如，[0 表示范围 (10( 中的 _]。当然，python有更好的方法"[0]*10"。
• 我注意到您导入了 numpy，但您没有对 numpy 数组使用操作。这将大大改进您的代码(并使其更快(并使其更易于阅读。例如，可以一次写入和读取整行或整列。"矩阵[0，：]"(第一行(，"矩阵[：，0]"(第一行(。
没有必要在 python 中包含对象的长度，
• 因为 python 会自动存储其对象的长度，这些长度始终可以使用内置的 len(( 进行检索。
• 对你的代码来说并不重要(因为你可能是为了练习而这样做(，但你可能已经知道，lu 分解已经存在，例如 scipy.linalg.lu((。事实上，快速检查 linalg.lu(A，True(会发现你的L和U都有错误。
• 看到你手动生成 L 和 U 矩阵，然后继续在 L 和 U 上使用 np.linalg.inv((，这真的很奇怪。如果你愿意使用 np.linalg.inv(( 函数，那么你的答案是一行，"np.linalg.inv(A( @ b"。通常，人们会更容易找到L和你手动求解X。找到L和你，然后使用numpy的反函数，有点违背了目的。
• 虽然它有时会有所帮助，但 python 不需要您在创建对象之前在内存中开辟空间。Python 会自动管理您的内存创建和删除(无需创建空的零列表(。
• np.dot(( 只是在 Python 中访问 "@" 运算符的手动方式。

一些例子：

import numpy as np
from scipy import linalg
def lu_dec_solver_v1(A,b):
return np.linalg.inv(A) @ b
def lu_dec_solver_v2(A,b):
L, U  = linalg.lu(A,True)
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)
def lu_dec_solver_v3(A,b):
U = A.copy()
L = np.identity(len(A))
for n in range(0,len(A)-1):
for m in range(n+1,len(A)):
L[m,n]  = U[m,n]/U[n,n]
U[m,:] += -L[m,n]*U[n,:]
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)

Answer 2

如果你想自己实现LU分解，这是一个不需要任何外部库的代码。我把它附加到这里，以防你只是想复制它

def Abs(x):
'''
Input x  and you'll get the abs.
'''
#Not the best way, but it is an easy way to not use numpy
return (x**2)**0.5

def ind_max(row,N):
'''
Find the index of the maximum of a list (row) of lentgth N.
'''
_in=0
_max=row[0]
i=0
while i<N:#the end of the row should be included (convension in how I use LUP..)
if row[i]>_max:
_max=row[i]
_in=i
i+=1

return _in

def Sum(List,N):
'''
Calculates the sum of a List of size N
'''
s=0
for i in range(N):
s+=List[i]
return s

def index_swap(A,index_1,index_2):
'''
index_swap takes an array and interchanges 
A[index_1] with A[index_2].

Example:
A=[0,1,2,3]
index_swap(A,0,2)
A
>>[2, 1, 0, 3]
'''

tmp=A[index_1]
A[index_1]=A[index_2]
A[index_2]=tmp


def apply_permutations_vector(A,P,N):
'''
Applies the permutations given by P from LUP
to a list A of length N, and returns the result.
Example:
A=[1,2,5,8,3]
P=[2,4,0,3,1]
apply_permutations_vector(A,P,5)
>>[5, 3, 1, 8, 2]
'''
#that is, you make a list like this (P basically gives you the indices of A):)

Ap=[A[ P[i] ] for i in range(N)]
return Ap

def apply_permutations_matrix(M,P,N):
'''
Applies the permutations given by P from LUP
to a N*N array M of length N, and returns the result.

M=[
[ 0,  2,  2 , 3 , 5],
[-3, -1,  1 , 5 , 9],
[ 1, -1,  1 , 4 , 7],
[ 1, -1,  1 , 0 , 2],
[ 1, -1,  1 , 0 , 3]
]
P=[2,0,1,4,3]
apply_permutations_matrix(M,P,5)
>>[
[ 1, -1, 1, 4, 7],
[ 0,  2, 2, 3, 5],
[-3, -1, 1, 5, 9],
[ 1, -1, 1, 0, 3],
[ 1, -1, 1, 0, 2]
]
'''
Mp=[ [M[ P[i] ][j]for j in range(N)] for i in range(N) ]

return Mp


def LUP(M,N,_tiny=1e-20):
U=[  [ M[i][j] for j in range(N)] for i in range(N) ]
L=[  [ 0 if i!=j else 1 for j in range(N)] for i in range(N) ]
#this is the "permutation vector". if it is e.g. [2 1 0 3] it means you make 0<->2
P=[  i for i in range(N) ]

for k in range(1,N):
for i in range(k,N):
#find the index of the maximum in column
_col=[Abs(U[_r][k-1]) for _r in range(k-1,N)]

#find the index of the maximum of _col
# notice that the length of _col is N-(k-1)
len_col=N-(k-1)
pivot=ind_max( _col ,len_col) + k - 1 #convert the index of _col (it has a length of len_col) to the index of  a row of U   

##################################################
#if you remove it, then you get a lot of infinities
#it has to do with the fact that if U[pivot][k-1] <_tiny , then U[k-1][k-1] will be a zero, 
#L[i][k-1] explodes. 
#You are allowed to skip this i, then, because if U[pivot][k-1] <_tiny , then all U[i][k-1] are small!
#Check that this is true by  uncommenting print(_col)
if Abs(U[pivot][k-1]) < _tiny  :         
#print(_col)
break
###################################################
#if the maximum is not at k-1, swap!
if pivot != k-1 : 
# Permute rows k-1 and pivot in U

index_swap(P,k-1,pivot)

tmpU=[U[k-1][_r] for _r in range(k-1,N)]

#print(U)
for _r in range(k-1,N):
U[k-1][_r]=U[pivot][_r]
#print(U)
for _r in range(k-1,N):
U[pivot][_r]=tmpU[_r-(k-1)]#again we have to convert the index of tmpU
#print(U)
#print("=========================")
tmpL=[L[k-1][_r] for _r in range(k-1)]
#print(L)
for _r in range(k-1):
L[k-1][_r]=L[pivot][_r]
#print(L)
for _r in range(k-1):
L[pivot][_r]=tmpL[_r]

#print(L)
#print("========================")

L[i][k-1]=U[i][k-1]/U[k-1][k-1]


for j in range(k-1,N):
U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j]

return L,U,P