Gauss seidel在python中的实现

Question

import numpy as np
from scipy.linalg import solve_triangular as triSolve

#O(n) per iteration, so overall O(nN), good for large SPD/SDD matrices
def GS_iter(A, b, N):
    m = len(A)
    L = np.tril(A)
    P = L-A
    print(P)
    
    x = np.zeros(m)
    print(x)
    for k in range(N):
        x = triSolve(L,b+P@x, True)
        
    return x
        


#examples
A = np.array([[10,2,3,1],[1,10,0,1],[0.2,1,10,2],[0.1,3,3,10]])
b = np.array([1,2,1,0])

x = GS_iter(A,b,50000)

ans = A@x-b
print(ans)
print(np.linalg.norm(ans))

上面是我在Python中的Gauss-Seidel方法。由于某种原因，即使在矩阵A是严格对角占优的情况下，即使在50000次迭代后，它也不收敛。下面是在MATLAB中工作的相同实现：

function x = gSeidel(A,B,N)
    [n,~] = size(A);
    L = tril(A);
    P = L-A; %P = -U
    
    x = zeros(n,1); %x_0
    for k = 1:N
        x = L\(B+P*x);
    end
end 

我犯了什么错误？我认为它是在TriSolve方法中，因为如果我用常规的LU求解器（如np.linalg.solve）替换它，它就能工作。为什么三角解的行为不像这里想要的那样？

Answer 1

lower参数是第四个用x=triSolve(L,b+p@x,lower=true)替换您的行

Signature:
triSolve(
    a,
    b,
    trans=0,
    lower=False,
    unit_diagonal=False,
    overwrite_b=False,
    debug=None,
    check_finite=True,
)