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]