python 做 Gauss-Seidal 迭代

#coding:utf-8
import numpy as np
import scipy.linalg as sp


# jacobi diedai
a=np.array([[2,-1,1],[1,1,1],[1,1,-2]])  #AX=B
b=np.array([[1,1,1]])
print 'straight resout=',sp.solve(a,b.T)
print 'jacobi resout=',sp.inv(a).dot(b.T)
l1=sp.tril(a) #上三角矩阵
u1=sp.triu(a) #下三角矩阵
d=l1+u1-a
l=d-l1
u=d-u1




bj=sp.inv(d-l).dot(u) #Gauss-Seidel迭代收敛到充分必要条件就是(D-L）到逆矩阵X （U）的最大特征值小于1
print 'bj=', bj
g=sp.inv(d-l).dot(b.T)
print 'g=', g
lamda=sp.eigvals(bj).max()


x0=np.array([[0,0,0]]).T
print 'x0=' ,x0
if lamda<1: