我有一个python代码，使用Numpy和Scipy，用Gauss-Seidel方法求解线性系统。我正在实现代码和书中的一个例子：* *’数值分析：《负担与仙子》***问题是我得到了精确的解，但需要更多的迭代：10次迭代，公差为0.0000001，但本书仅用6次迭代和0.001公差就得到了解决方案。***我认为问题是因为使用scipy计算误差的无穷范数。当我不使用代码中的错误（只有迭代）时，我获得了与书中相同的结果。下面是我的Python代码：
import numpy as np
import scipy as scp
def gauss_seidel(A, b, x_0, max_iterations=15, tolerance=0.0000001):
L = -np.tril(A, -1)
U = -np.triu(A, 1)
v = np.diagonal(A)
D = np.diag(v)
DL = D - L
Hg = np.linalg.inv(DL)
Tg = Hg @ U
Cg = Hg @ b
n = A.shape[0]
x = np.zeros(n)
diff = np.zeros(n)
error = 0.0
k = 1
while k <= max_iterations:
x = Tg @ x_0 + Cg
diff = x - x_0
error = scp.linalg.norm(diff, ord=np.inf, axis=None) / \
scp.linalg.norm(x, ord=np.inf)
x_0 = x
k += 1
if(error < tolerance):
break
return x, k
A = np.matrix([
[10, -1, 2, 0],
[-1, 11, -1, 3],
[2, -1, 10, -1],
[0, 3, -1, 8]
])
b = np.array([6, 25, -11, 15])
x_0 = np.array([0, 0, 0, 0])
solution = gauss_seidel(A, b, x_0, tolerance=0.001)
print('WITH TOLERANCE = 0.001')
print(
f'Solution = {solution[0]} with {solution[1]} iterations')
solution = gauss_seidel(A, b, x_0)
print('WITH TOLERANCE = 0.0000001')
print(
f'Solution = {solution[0]} with {solution[1]} iterations')这是我的终端输出：
WITH TOLERANCE = 0.001 Solution = [ 1.00009128 2.00002134 -1.00003115 0.9999881 ] with 6 iterations WITH TOLERANCE = 0.0000001 Solution = [ 1. 2. -1.1.]，10次迭代
谢谢
1条答案
按热度按时间gudnpqoy1#
您的GS求解器实现得很好！但是：这很重要吗？你自己跑GS？如果不是，那么不同的scipy求解器通常比GS快得多，即使具有更严格的收敛标准。
对于这个小方程组，直接求解器
np.linalg.solve()和scp.linalg.solve()比GS快10到20倍。对于更大的系统[25,25]和更多，迭代求解器更适合。bicgstab将比GS快得多。（但是这个例子太小了，不能有效）。