In the following code for the Gauss Seidel method, I enter one given matrix A. The result seems to be correct, but when I comment the vector x1 at the beginning of the while, I obtain an unwanted result:
For example, before the assignment x0=x1, when k=1, x0 is equal to x1; instead x0 when k=1, would be equal to x1 when k=0.
Consequently, the norm(x1-x0) is always 0, after the first while. I don't know why it happens.
This is the code, you can see the wanted and unwanted output below:
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))
Wanted output: - As you can see x0 contains the x1 of of the previous while iteration
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
Unwanted output if I comment x1 = np.zeros(n) at the beginning of 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
Even if I don't assing zeros to x1 in each cycle, the solution are computed correctly. Can you help me?
You can execute online at: https://www.jdoodle.com/a/1h6N
I am not entirely sure what are you trying to achieve and why do you want to AVOID zeroing out x1 if you care about what is being printed.
However, I think I see the source of confusion. Notice, that in the following piece of code:
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)
the print expression is inside the for-loop, while the entries of x1 are updated one by one.
So, if you don't zero out x1 before the for-loop:
while iteration, x1 entries contain garbage, so the first print output will have a correct entry of x1[0] and garbage for x1[1:].while iterations x1 entries would contain previous (not yet overwritten) values.So, if you want to see the "progress inside" without confusion on data that will be overwritten, you should zero it out. Otherwise, with this particular algorithm implementation it is safe to leave it out if you don't care about the output.