时间戳

首页 > 解决方案 > 高斯赛德尔在python中的实现

python - 高斯赛德尔在python中的实现

问题描述

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) 替换它）它可以工作。为什么三角求解在这里没有按预期运行？

标签： pythonmatlabnumpynumerical-methods

解决方案

较低的参数是第四个将您的行替换为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,
)

推荐阅读