线性方程组矩阵分解Doolittle

YandexOj

已于 2022-04-19 13:28:51 修改

阅读量151

点赞数

文章标签：矩阵线性代数算法

于 2021-11-27 16:17:04 首次发布

本文链接：https://blog.csdn.net/weixin_47989730/article/details/121577464

版权

矩阵分解

import numpy as np
def LU_decomposition(A):
    n=len(A[0])
    L = np.zeros([n,n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i]=1
        if i==0:
            for j in range(0,n):
                U[0][j]=A[0][j]
                L[j][0]=A[j][0]/U[0][0]
        else:
                for j in range(i, n):#U
                    temp=0
                    for k in range(0, i):
                        temp = temp+L[i][k] * U[k][j]
                    U[i][j]=A[i][j]-temp
                for j in range(i+1, n):#L
                    temp = 0
                    for k in range(0, i ):
                        temp = temp + L[j][k] * U[k][i]
                    L[j][i] = (A[j][i] - temp)/U[i][i]
    return L,U
 
if __name__ == '__main__': 
    #A=np.random.randint(1,10,size=[3,3])  #注意A的顺序主子式大于零
    n = int(input('输入矩阵长度:'))
    list = [[]for i in range(n)]
    for i in range(n):
        line = input().split(' ')
        for j in range(len(line)):
            list[i].append(int(line[j]))
    A=np.array(list)
    L,U=LU_decomposition(A)
    print(L,'\n',U)

解题

import numpy as np
def LU_decomposition(A):
    n=len(A[0])
    L = np.zeros([n,n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i]=1
        if i==0:
            for j in range(0,n):
                U[0][j]=A[0][j]
                L[j][0]=A[j][0]/U[0][0]
        else:
                for j in range(i, n):#U
                    temp=0
                    for k in range(0, i):
                        temp = temp+L[i][k] * U[k][j]
                    U[i][j]=A[i][j]-temp
                for j in range(i+1, n):#L
                    temp = 0
                    for k in range(0, i ):
                        temp = temp + L[j][k] * U[k][i]
                    L[j][i] = (A[j][i] - temp)/U[i][i]
    return L,U
def cauculate(L,U,b):
    Y=np.matmul(np.linalg.inv(L),b)
    X=np.matmul(np.linalg.inv(U),Y)
    return X

if __name__ == '__main__': 
    #A=np.random.randint(1,10,size=[3,3])  #注意A的顺序主子式大于零
    n = int(input('输入矩阵长度:'))
    list = [[]for i in range(n)]
    for i in range(n):
        line = input().split(' ')
        for j in range(len(line)):
            list[i].append(int(line[j]))
    A=np.array(list)
    print("输入矩阵b")
    blist = input().split(' ')
    b=np.array(blist,dtype=np.double)
    L,U=LU_decomposition(A)
    X=cauculate(L,U,b)
    print(L,'\n',U)
    print(X)