Метод Гаусса-Зайделя в Python3, в начале каждого цикла, почему я должен обнулять массив с помощью самых последних решений? - PullRequest
Новая
       42

Метод Гаусса-Зайделя в Python3, в начале каждого цикла, почему я должен обнулять массив с помощью самых последних решений?

0 голосов
/ 11 июня 2019

В следующем коде для метода Гаусса Зейделя я ввожу одну данную матрицу A.Результат кажется правильным, но когда я комментирую вектор x1 в начале while, я получаю нежелательный результат:

Например, перед присвоением x0=x1, когда k=1, x0 равно x1;вместо x0, когда k=1, будет равно x1, когда k=0.

Следовательно, norm(x1-x0) всегда равен 0 после первого while.Я не знаю, почему это происходит.

Этот код, вы можете увидеть желаемый и нежелательный вывод ниже: 

   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

Нежелательный вывод, если я комментирую x1 = np.zeros(n) в начале while: 

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

1 Ответ

1 голос
/ 12 июня 2019

Я не совсем уверен, чего вы пытаетесь достичь и почему вы хотите ИЗБЕЖАТЬ обнуления 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 обновляются одна за другой. Итак, если вы не обнуляете x1 до for -loop:

  • на первой итерации while записи x1 содержат мусор, поэтому первый вывод print будет иметь правильную запись x1[0] и мусор для x1[1:].
  • на других while итерациях x1 записи будут содержать предыдущие (еще не перезаписанные) значения.

Итак, если вы хотите увидеть «прогресс внутри» без путаницы с данными, которые будут перезаписаны, вы должны обнулить их. В противном случае, с этой конкретной реализацией алгоритма, можно оставить ее без внимания, если вас не волнует вывод.

...