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. ]
