In the following code I have made an SOR iteration in python. I seem to be getting the wrong output. It seems to be just outputting the b array when it should be solving the system. Is there something wrong with my function or call of my function? Here is the code below! I don't think my values of x should be all 1's and I don't think it should just take one iteration! When I run this code for smaller x, b, xo it works! but I am trying to run it for x which is a 10x10, b which is a 10x1 and xo which is a 10x1
import numpy as np
import math
x = np.array([[3.0, 1.0, 0., 0., 0., 0., 0., 0., 0., 0.],[1.0, 3.0, 1.0, 0., 0., 0., 0., 0., 0., 0.], [0., 1.0, 3.0, 1.0, 0., 0., 0., 0., 0., 0.], [0., 0, 1.0, 3.0, 1.0, 0., 0., 0., 0., 0.], [0., 0., 0., 1.0, 3.0, 1.0, 0., 0., 0., 0.], [0., 0., 0., 0., 1.0, 3.0, 1.0, 0., 0., 0.], [0., 0., 0., 0., 0., 1.0, 3.0, 1.0, 0., 0.], [0., 0., 0., 0., 0., 0., 1.0, 3.0, 1.0, 0.], [0., 0., 0., 0., 0., 0., 0., 1.0, 3.0, 1.0], [0., 0., 0., 0., 0., 0., 0., 0., 1.0, 3.0]])
b = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
x0 = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
tol = 10 ** (-15)
max_iter = 20
w = 1.5
def SOR(A, b, x0, tol, max_iter, w):
if (w<=1 or w>2):
print('w should be inside [1, 2)');
step = -1;
x = float('nan')
return
n = b.shape
x = x0
for step in range (1, max_iter):
for i in range(n[0]):
new_values_sum = np.dot(A[i, 1 : (i - 1)], x[1 : (i - 1)])
for j in range(i + 1, n[0]):
old_values_sum = np.dot(A[i, j], x0[j])
x[i] = b[i] - (old_values_sum + new_values_sum) / A[i, i]
x[i] = np.dot(x[i], w) + np.dot(x0[i], (1 - w))
if (np.linalg.norm(x - x0) < tol):
print(step)
break
x0 = x
print("X = {}".format(x))
print("The number of iterations is: {}".format(step))
SOR(x, b, x0, tol, max_iter, w)
Which gives the following output, which is not what I am looking for
1
X = [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
The number of iterations is: 1
You set
x = x0. This means that both variables reference the same object, changes to one are automatically changes to the other. You need to force a copy, this happens in one of the variantswhere the last is obviously the preferred version for a descriptive programming style.
Later you assign
x0=xwhich has the same problem. Here you might not want to create a new copy and delete the previous object referenced byx0, so you can useto force copying of the data into the existing object.