如果在 doolittle 算法中对角线元素为 0，则除以零误差 [英] Divide by zero error if diagonal elements are 0 in doolittle algorithm

问题描述

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

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

我得到一个除以零错误"因为我的对角线元素为零.是doolittle算法的限制吗?
我能做些什么来解决它?

So, I am trying to create a code which does LU Decomposition of a matrix by doolittle algorithm in python.
However, I am getting a divide by zero error if my diagonal elements are zero.
For example, when i tried solving this as my input matrix:

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

I get a "divide by zero error" as my diagonal elements are zero. Is it a limitation of doolittle algorithm?
What can I do to solve it?

PS:- 我知道为什么我会得到除以零错误"，我想知道的是，我现在能做什么?
我的代码到目前为止:

P.S:- I know why I am getting the "divide by zero error", what I want to know is, what can I do now?
My code So far:

#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 分解 - 或者更好的说法是它的一个或多个绝对对角线值不大于其余(包括零场景)，在该行中 - 您需要首先对其进行置换.

To initiate LU decomposition for a table that is not strictly diagonally dominant - or better put that one or more of its absolute diagonal values is not greater than the sum of the rest (includes the zero scenario), in that row - you would need to firstly permute it.

这种方式允许您交换初始 A 矩阵中的行和印记"通过同时排列其行，将它们放在一个单位矩阵(与 A 的大小相同)上.因此，您最终得到 PA = Pb，它扩展到 P(LU) = Pb，其中 P 是您的单位置换矩阵，其中 A = LU 是您的初始系数矩阵，而 b 是您的初始向量.请参阅下面的代码作为置换矩阵的示例:

This way allows you to exchange rows in your initial A matrix and "imprint" them on a unit matrix (of same size to your A) by simultaneously permuting its rows. Thus you end up having PA = Pb, which extends to P(LU) = Pb where P is your unit permutation matrix where A = LU is your initial coefficient matrix and b is your initial vector. See the code below as an example of a permutation matrix:

#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 分解来自德克萨斯大学，以及吉尔伯特·斯特朗 (Gilbert Strang) 的线性代数及其应用"一书中第 4 版，第 41 页，子章节行交换和置换矩阵".

I also urge you to have a look here at LU factorization from University of Texas, as well as in the book of Gilbert Strang "Linear Algebra and Its Applications" 4th edition, page 41, the subchapter "Row Exchanges and Permutation Matrices".

这篇关于如果在 doolittle 算法中对角线元素为 0，则除以零误差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！