import numpy as np
import math
A=np.array([[1,2,3],[2,5,2],[3,1,5]])  # np.mat创建矩阵,np.arry创建数组
b=np.array([14,18,20])
n=len(A)
def LU_mat(A):   #分解矩阵A的函数，得到L和U:
    L=np.eye(len(A))
    U=np.zeros(np.shape(A))
    for k in range(n):
        U[0,0]=A[0,0]
        if k==0:
            for j in range(1,n):
                U[0,j]=A[0,j]
                L[j,0]=A[j,0]/U[0,0] 
        else:
            for j in range(k,n):
                m=0
                for r in range(k):
                    m+=L[k,r]*U[r,j]
                U[k,j]=A[k,j]-m
            for i in range(k+1,n):
                m=0
                for r in range(k):
                    m+=L[i,r]*U[r,k]
                L[i,k]=(A[i,k]-m)/U[k,k]
    print("L={}".format(L),"U={}".format(U))
    return L,U



def LU_fun(A,b):  #定义一个LU分解函数，自变量是A,b
    m1,m2=A.shape  # m,n分别代表矩阵A的行数和列数
    
    if A[0,0]==0:
        print("no answer")     
    if m1<m2:
        print("这是一个解空间")
    else:
        L,U=LU_mat(A)
        y=b
        for k in range(1,n):
            m=0
            for r in range(k):
                m+=L[k,r]*y[r]
            y[k]-=m
        print("y={}".format(y))
        x=y
        x[n-1]=y[n-1]/U[n-1,n-1]
        for i in range(n-2,-1,-1):
            m=0
            for k in range(i+1,n):
                m+=U[i,k]*x[k]
            x[i]=(x[i]-m)/U[i,i]
    print("x={}".format(x))
    return x

LU_fun(A,b)

运行结果

L=[[ 1.  0.  0.]
   [ 2.  1.  0.]
   [ 3. -5.  1.]] 
U=[[  1.   2.   3.]
   [  0.   1.  -4.]
   [  0.   0. -24.]]
y=[ 14 -10 -72]
x=[1 2 3]