如果在doolittle算法中对角元素为0，则除以零误差

问题描述

因此，我正在尝试创建一个代码，该代码通过python中的doolittle算法对矩阵进行LU分解。
但是，如果对角元素为零，我将得到除以零的错误。
例如，当我尝试解决此问题作为输入矩阵时：

A = [[0 2 8 6] [0 0 1 2] [0 1 0 1] [3 7 1 0]]

因为对角元素为零，所以我得到了“除以零误差”。 是doolittle算法的局限性吗？
我该怎么解决？

P.S：-我知道为什么会出现“除以零误差”的情况，我想知道的是，我现在该怎么办？
到目前为止，我的代码：

#Function to read the file and getting numbers in a Matrix
def file_opener(file_name):
    Matrix = []
    with open(file_name) as file:
        M_string = file.readlines()
    for line in M_string:
        Matrix.append(list(map(lambda i: int(i),line.split(" "))))
    return Matrix

#Function to create Identity Matrix
def identity(n):
    I = [[0 for y in range(n)] for x in range(n)]
    for i in range(n):
        I[i][i] = 1
    return I


#Function for LU decomposition
#Both L and U are directly kept together in A
#Note:- Diagonal elements of L should be 1 but are not stored in A
def luDecompose(A,n):
    for i in range(n):
        # Upper Triangle Matrix (i is row index,j is column index,k is summation index)
        for j in range(i,n):
            # Summation part
            sum = 0
            for k in range(i):
                if(i==k):
                    sum += A[k][j]  # Since diagonal elements of Lower matrix is 1
                else:
                    sum += A[i][k]*A[k][j]
        
            A[i][j] = A[i][j] - sum
        
        # Lower Triangle Matrix (j is row index,i is column index,k is summation index)
        for j in range(i+1,n):
            # Summation part
            sum = 0
            for k in range(i):
                if(j==k):
                    sum += A[k][i]  # Since diagonal elements of Lower matrix is 1
                else:
                    sum += A[j][k]*A[k][i]
            A[j][i] = (A[j][i] - sum)/A[i][i]


#Function for Forward-Backward Substitution
#Works in complimentary to LU Decomposition function
def forwardBackwardSubstitution(A,B,m,n):
    #Forward Substitution
    for j in range(n):
        for i in range(m):
            sum = 0
            for k in range(i):
                sum += A[i][k]*B[k][j]
            B[i][j] = B[i][j] -  sum

    #Backward Substitution
    for j in range(n):
        for i in range(m-1,-1,-1):
            sum = 0
            for k in range(i+1,m):
                sum += A[i][k]*B[k][j]
            B[i][j] = (B[i][j] - sum)/A[i][i]


# The solving process starts here
A = file_opener("Matrix-A2.txt")
sizeA = len(A[0])

#Creating Identity matrix
I = identity(sizeA)

luDecompose(A,sizeA)
forwardBackwardSubstitution(A,I,sizeA,sizeA)
print(I)

解决方法

要为不是严格对角占优的表启动 LU 分解 - 或者更好地将其一个或多个绝对对角线值不大于其余的总和（包括零场景），在该行中 - 您需要先对它进行置换。

这种方式允许您交换初始 A 矩阵中的行，并通过同时置换其行将它们“印刻”在单位矩阵（与 A 的大小相同）上。因此，您最终得到 PA = Pb，它扩展到 P(LU) = Pb，其中 P 是您的单位置换矩阵，其中 A = LU 是您的初始系数矩阵，而 b 是您的初始向量。请参阅下面的代码作为置换矩阵的示例：

#Permutation Matrix
import numpy as np

def permute(matrix):
    n = matrix.shape[0]
    U = matrix.copy()
    P = np.identity(n)
    for k in range(n):
        index = np.argmax(abs(U[k:,k]))
        index = index + k 
        if index != k:
            P = np.identity(n)
            P[[index,k],k:n] = P[[k,index],k:n]

    return P

我还敦促您查看 LU factorization from University of Texas 以及 Gilbert Strang “线性代数及其应用”第 4 版第 41 页的子章节“行交换和置换矩阵” .

,

所以，我明白了。 上面的矩阵不能直接进行分解。 因此，我们必须进行一些行交换，并将交换信息存储在另一个数组中。然后可以分解矩阵。