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:
while 1:
x1=bj.dot(x0)+g
print 'x1=',x1
if (x1-x0).all()<1e-10: #判断迭代的停止条件
break
x0=x1                   

print 'Jacobi resout=', x1

print 'real resout', sp.solve(a,b.T)