import numpy as np#分解矩阵A至L和Udef LU_mat(A):L=np.zeros(np.shape(A))U=np.eye(len(A))L[0,0]=A[0,0]U[0,1]=A[0,1]/A[0,0]for i in range(1,A.shape[0]-1):L[i,i-1]=A[i,i-1]L[i,i]=A[i,i]-A[i,i-1]*U[i-1,i]U[i,i+1]=A[i,i+1]/(A[i,i]-A[i,i-1]*U[i-1,i])n=A.shape[0]-1L[n,n-1]=A[n,n-1]L[n,n]=A[n,n]-A[n,n-1]*U[n-1,n]print(L,U)return L,Udef Solve(A,f):L,U=LU_mat(A)y=np.zeros((f.shape[0],1))x=np.zeros((f.shape[0],1))y[0,0]=f[0,0]/L[0,0]for i in range(1,L.shape[0]):y[i,0]=(f[i,0]-L[i,i-1])/L[i,i]x[U.shape[0]-1]=y[U.shape[0]-1,0]for j in range(U.shape[0]-2,-1,-1):x[j,0]=y[j,0]-U[j,j+1]*x[j+1,0]return xA = np.array([[4,1,0],[0,3,2],[0,-1,3]])f = np.array([[2],[9],[8]])print("原系数矩阵A:")print(A, "\n")print("f:")print(f, "\n")print("最终求解结果:")print(Solve(A, f))
