Python3 中的 Gauss-Seidel 方法，为什么我必须在每个循环开始时将最新解的数组设置为零？

Question

在 Gauss Seidel method 的以下代码中，我输入了一个给定的矩阵 A。结果似乎是正确的，但是当我在 while 的开头注释向量 x1 时，我得到了一个不需要的结果：

比如赋值前x0=x1，当k=1时，x0等于x1；相反 x0 当 k=1 时等于 x1 当 k=0.

因此，第一个 while 之后的 norm(x1-x0) 始终为 0。我不知道为什么会这样。

这是代码，您可以在下面看到需要和不需要的输出：

   def GaussSeidel(A,b):
       # dimension of the non-singular matrix
       n = len(A)

       # def. max iteration and criterions
       Kmax = 100;
       tol  = 1.0e-4;
       btol = la.norm(b)*tol


       x0   = np.zeros(n)
       k    = 0 ;
       stop = False
       x1   = np.empty(n)

       while not(stop) and k < Kmax:
           print ("begin while with k =", k)
           x1 = np.zeros(n)
           for i in range(n):          # rows of A
               x1[i] = ( b[i] - np.dot(A[i,0:i], x1[0:i]) - np.dot(A[i,i+1:n], x0[i+1:n]) ) / A[i,i]
               print("x1 =", x1)

           r    = b - np.dot(A,x1)
           stop = (la.norm(r) < btol) and (la.norm(x1-x0) < tol)
           print("end of for i \n")
           print("x0 =", x0)
           print("btol = %e \t; la.norm(r) = %e \t; tol = %e \t; la.norm(x1-x0) = %e; stop = %s " % (btol, la.norm(r), tol, la.norm(x1-x0), stop))
           x0   = x1
           print("x0 =", x0, end='\n')
           print("end of current while \n\n")
           k    = k + 1

       if not(stop): # or if k >= Kmax
           print('Not converges in %d iterations' % Kmax)

       return x1, k

   import numpy        as np
   import numpy.linalg as la
   import time

   A = np.array( [
          [  3, -0.1, -0.2],
          [0.1,    7, -0.3],
          [0.3, -0.2,   10]
       ], dtype='f')

   b = np.array( [7.85, -19.3, 71.4] )

   xsol = la.solve(A,b)

   start    = time.time()
   x, k     = GaussSeidel(A,b)
   ending   = time.time()
   duration = ending-start
   err      = la.norm(xsol-x)
   print('Iter.=%d  duration=%f  err=%e' % (k,duration,err))

想要的输出：- 如您所见，x0 包含前一次迭代

的 x1

begin while with k = 0
x1 = [2.61666667 0.         0.        ]
x1 = [ 2.61666667 -2.79452381  0.        ]
x1 = [ 2.61666667 -2.79452381  7.00560952]
end of for i 

x0 = [0. 0. 0.]
la.norm(r) = 2.382271e+00   ; la.norm(x1-x0) = 7.983412e+00; stop = False 
x0 = [ 2.61666667 -2.79452381  7.00560952]
end of current while 


begin while with k = 1
x1 = [2.99055651 0.         0.        ]
x1 = [ 2.99055651 -2.49962467  0.        ]
x1 = [ 2.99055651 -2.49962467  7.00029081]
end of for i 

x0 = [ 2.61666667 -2.79452381  7.00560952]
la.norm(r) = 2.847092e-02   ; la.norm(x1-x0) = 4.762220e-01; stop = False 
x0 = [ 2.99055651 -2.49962467  7.00029081]
end of current while 


begin while with k = 2
x1 = [3.0000319 0.        0.       ]
x1 = [ 3.0000319  -2.49998798  0.        ]
x1 = [ 3.0000319  -2.49998798  6.99999928]
end of for i 

x0 = [ 2.99055651 -2.49962467  7.00029081]
la.norm(r) = 1.288604e-04   ; la.norm(x1-x0) = 9.486833e-03; stop = False 
x0 = [ 3.0000319  -2.49998798  6.99999928]
end of current while 


begin while with k = 3
x1 = [3.00000036 0.         0.        ]
x1 = [ 3.00000036 -2.50000002  0.        ]
x1 = [ 3.00000036 -2.50000002  6.99999998]
end of for i 

x0 = [ 3.0000319  -2.49998798  6.99999928]
la.norm(r) = 1.084102e-06   ; la.norm(x1-x0) = 3.377360e-05; stop = True 
x0 = [ 3.00000036 -2.50000002  6.99999998]
end of current while 


Iter.=4  duration=0.234001  err=3.544580e-07

如果我在 while 的开头注释 x1 = np.zeros(n) 则不需要的输出：

begin while with k = 0
x1 = [2.61666667e+000 1.94626056e+227 2.04746603e+161]
x1 = [ 2.61666667e+000 -2.79452381e+000  2.04746603e+161]
x1 = [ 2.61666667 -2.79452381  7.00560952]
end of for i 

x0 = [0. 0. 0.]
la.norm(r) = 2.382271e+00   ; la.norm(x1-x0) = 7.983412e+00; stop = False 
x0 = [ 2.61666667 -2.79452381  7.00560952]
end of current while 


begin while with k = 1
x1 = [ 2.99055651 -2.79452381  7.00560952]
x1 = [ 2.99055651 -2.49962467  7.00560952]
x1 = [ 2.99055651 -2.49962467  7.00029081]
end of for i 

x0 = [ 2.99055651 -2.49962467  7.00029081]
la.norm(r) = 2.847092e-02   ; la.norm(x1-x0) = 0.000000e+00; stop = False 
x0 = [ 2.99055651 -2.49962467  7.00029081]
end of current while 


begin while with k = 2
x1 = [ 3.0000319  -2.49962467  7.00029081]
x1 = [ 3.0000319  -2.49998798  7.00029081]
x1 = [ 3.0000319  -2.49998798  6.99999928]
end of for i 

x0 = [ 3.0000319  -2.49998798  6.99999928]
la.norm(r) = 1.288604e-04   ; la.norm(x1-x0) = 0.000000e+00; stop = True 
x0 = [ 3.0000319  -2.49998798  6.99999928]
end of current while 


Iter.=3  duration=0.156000  err=3.409068e-05

即使我没有在每个循环中将零分配给 x1，解决方案也会正确计算。你能帮帮我吗？

您可以在线执行： https://www.jdoodle.com/a/1h6N

Answer 1

我不完全确定你想要实现什么，如果你关心正在打印的内容，为什么要避免归零 x1。

但是，我想我看到了混乱的根源。请注意，在以下代码中：

x1 = np.zeros(n)
    for i in range(n):          # rows of A
         x1[i] = ( b[i] - np.dot(A[i,0:i], x1[0:i]) - np.dot(A[i,i+1:n], x0[i+1:n]) ) / A[i,i]
         print("x1 =", x1)

print 表达式在 for 循环内，而 x1 的条目被逐一更新。因此，如果您不在 for 循环之前将 x1 归零：

在第一次 while 迭代中，x1 条目包含垃圾，因此第一个 print 输出将具有正确的 x1[0] 条目和 [=21 的垃圾=].
在其他 while 次迭代中 x1 条目将包含以前（尚未覆盖）的值。

因此，如果您想查看 "progress inside" 而不会混淆将被覆盖的数据，则应将其清零。否则，对于这个特定的算法实现，如果您不关心输出，可以安全地将其省略。

Python3 中的 Gauss-Seidel 方法，为什么我必须在每个循环开始时将最新解的数组设置为零？

Gauss-Seidel method in Python3, at the beginning of every cycle why I have to set to zero the array with most recent solutions?

python-3.6

numerical-methods