Is there a linear algebra libra开发者_运维技巧ry that implements iterative Gauss-Seidel to solve linear systems? Or maybe a preconditioned gradient solver?
Thanks
EDIT: In the end I used a kind of crude but correct way to solve it. As i had to create the matrix A (for Ax=b) anyway, I partitioned the matrix as
A = M - N
with
M = (D + L) and N = -U
where D is the diagonal, L is the lower triangular section, and U the upper triangular section. Then
Pinv = scipy.linalg.inv(M)
x_k_1 = np.dot(Pinv,np.dot(N,x_k)) + np.dot(Pinv,b)
Also did some convergence tests. It works.
I know this is old but, I haven't found any pre existing library in python for gauss - seidel. Instead I created my own little function that with the help of a permutation matrix as seen in another answer of mine permutation matrix will produce the solution (x vector) for any square matrix, including those with zeros on the diagonal.
def gauss_seidel(A, b, tolerance, max_iterations, x):
#x is the initial condition
iter1 = 0
#Iterate
for k in range(max_iterations):
iter1 = iter1 + 1
print ("The solution vector in iteration", iter1, "is:", x)
x_old = x.copy()
#Loop over rows
for i in range(A.shape[0]):
x[i] = (b[i] - np.dot(A[i,:i], x[:i]) - np.dot(A[i,(i+1):], x_old[(i+1):])) / A[i ,i]
#Stop condition
#LnormInf corresponds to the absolute value of the greatest element of the vector.
LnormInf = max(abs((x - x_old)))/max(abs(x_old))
print ("The L infinity norm in iteration", iter1,"is:", LnormInf)
if LnormInf < tolerance:
break
return x
After re searching I have found the following that could be used as a ready to go package, but haven't used it myself numerica PyPI
from sage.all import *
a=matrix(QQ,[[12,3,-5],[1,5,3],[3,7,13]])
b=matrix(QQ,[[1],[28],[76]])
x=[]
r=len(a[0])
i=0
while(i<r):
li=raw_input('give approximate solution :')
h=eval(li)
x.append(h)
i=i+1
def gausseidel():
tol=0.00001;l=0;itern=10
fnd=0
while(l<itern and fnd==0):
j=0
while(j<r):
temp=b[j,0]
k=0
while(k<r):
if (j!=k):
temp-=a[j,k]*x[k]
k=k+1
c=temp/a[j,j]
if (abs(x[j]-c)<=tol):
fnd=1
break
x[j]=c
j=j+1
l=l+1
if l==itern:
print "Iterations are over"
for p in range(r):
print 'x',p,'=',float(x[p])
gausseidel()
精彩评论