我通过以下链接实现了矩阵LU分解的标准方程/算法:(1)和(2)
这将返回如下所示的方形矩阵的LU分解。
然而,我的问题是——它也给出了一个Divide by Zero warning
代码如下:
import numpy as np
def LUDecomposition (A):
L = np.zeros(np.shape(A),np.float64)
U = np.zeros(np.shape(A),np.float64)
acc = 0
L[0,0]=1
for i in np.arange(len(A)):
for k in range(i,len(A)):
for j in range(0,i):
acc += L[i,j]*U[j,k]
U[i,k] = A[i,k]-acc
for m in range(k+1,len(A)):
if m==k:
L[m,k]=1
else:
L[m,k] = (A[m,k]-acc)/U[k,k]
acc=0
return (L,U)
A = np.array([[-4, -1, -2],
[-4, 12, 3],
[-4, -2, 18]])
L, U = LUDecomposition (A)
我哪里做错了
似乎您可能在第一个内部级别
for
循环方面出现了一些缩进错误:U
必须在L
之前进行计算;您还没有正确计算求和项acc
,并且没有正确地将L
的对角项设置为1。在进行一些其他语法修改后,您可以按如下方式重写函数:与scipy.linalg.lu相比,这一次为矩阵
A
提供了正确的输出,作为可靠的参考:注意:与scipy lapack getrf算法不同,这个Doolittle实现不包括旋转,只有当
scipy.linalg.lu
返回的置换矩阵P
是一个单位矩阵时,这两个比较才是正确的,即,scipy没有执行任何置换,这确实是矩阵A
的情况。scipy算法中确定的置换矩阵旨在优化生成矩阵的条件数,以减少舍入误差。最后,您只需验证A = LU
,如果因子分解正确,则始终如此:然而,就数值效率和精度而言,我不建议您使用自己的函数来计算LU分解。希望这有帮助
相关问题 更多 >
编程相关推荐