1 # coding:utf8 2 import numpy as np 3 4 def lu(mat): 5 r,c=np.shape(mat) 6 s=min(r,c) 7 for k in range(s): 8 x=1.0/mat[k][k] # 将后续除法变成乘法 9 for i in range(k+1,r): 10 mat[i][k]=mat[i][k]*x # L[1:][0]*U[0][0]=A[1:][0]；A[0][:]=mat[0][:] 11 for i in range(k+1,r): 12 for j in range(k+1,c): 13 # U[1][2]*L[1][1]=A[1][2]-U[0][2]*L[1][0]；L[1][1]=1 14 # L[2][1]*U[1][1]=A[2][1]-U[0][1]*L[2][0]；下一个k时mat[i][j]/mat[k][k](i>j) 15 mat[i][j]=mat[i][j]-mat[k][j]*mat[i][k] 16 return mat,c 17 18 def solve(A,b): 19 mat,n=lu(A) # LU合并 20 print mat # [[16, 4, 8], [0.25, 4.0, -6.0], [0.5, -1.5, 9.0]] 21 Z= np.zeros(n) # L*Z=b U*X=Z 22 X= np.zeros(n) 23 Z[0]=b[0] 24 for i in range(1,n): 25 sumup=0 26 for tmp in range(0,i): 27 sumup+=mat[i][tmp]*Z[tmp] 28 Z[i]=(b[i]-sumup) 29 X[n-1]=Z[n-1]/mat[n-1][n-1] 30 for i in range(n-2,-1,-1): 31 sumup=0 32 for tmp in range(i+1,n): 33 sumup+=mat[i][tmp]*X[tmp] 34 X[i]=(Z[i]-sumup)/mat[i][i] 35 return X 36 37 A=[[16,4,8],[4,5,-4],[8,-4,22]] 38 y=[-4,3,10] 39 print "The result of the fomula is:"+str(solve(A,y)) # [-2.25 4. 2. ]