用于高斯-赛德尔迭代求解的Python库？

Question

有没有线性代数的库可以用来实现迭代的高斯-赛德尔方法来解决线性方程组？或者有没有什么预处理的梯度求解器？

谢谢！

补充：最后我用了一种简单但正确的方法来解决这个问题。因为我本来就需要创建矩阵A（用于Ax=b），所以我把矩阵分成了

A = M - N

其中

 M = (D + L) and N = -U

D是对角线部分，L是下三角部分，U是上三角部分。然后

Pinv = scipy.linalg.inv(M)
x_k_1 = np.dot(Pinv,np.dot(N,x_k)) + np.dot(Pinv,b)

我还做了一些收敛性测试。结果是有效的。

Answer 1

这段内容是关于编程问题的讨论，涉及到一些技术细节和解决方案。具体来说，大家在讨论如何处理某个编程难题，分享各自的经验和看法。每个人都提出了自己的想法，有的人给出了代码示例，有的人则解释了为什么这样做有效。

总的来说，这里有很多实用的建议和技巧，适合那些刚开始学习编程的人。通过这些讨论，初学者可以更好地理解如何解决类似的问题，并且学到一些编程的基本知识。

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()

Answer 2

我知道这个问题有点老了，但我在Python里没有找到现成的高斯-赛德尔库。于是我自己写了一个小函数，利用一个排列矩阵，像我之前的回答中提到的排列矩阵，可以为任何方阵（包括对角线上有零的情况）计算出解（x向量）。

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

经过重新搜索，我发现有一个可以直接使用的包，但我自己还没用过，链接在这里numerica PyPI