from __future__ import division, absolute_import, print_function
import numpy as np
def p_print(x):
return str(x.tolist()).replace(']','').replace('[', '').replace(',','')
def seidel(A, b, x0, eps, maxiter=100):
x = np.zeros(b.shape)
B = np.zeros(A.shape)
c = np.zeros(b.shape)
n = 0
for i in range(A.shape[0]):
c[i] = b[i] / A[i, i]
for j in range(A.shape[1]):
if i != j:
B[i, j] = -A[i, j] / A[i, i]
B1 = np.tril(B, k=0)
B2 = np.triu(B, k=0)
eps2 = (1 - np.linalg.norm(B, ord=np.inf)) / np.linalg.norm(B2, ord=np.inf) * eps
err = 2 * eps2
while (n < maxiter) and (err > eps2):
for i in range(B.shape[0]):
x[i] = np.dot(x0, B2[i]) + np.dot(x, B1[i]) + c[i]
err = np.linalg.norm(x - x0, ord=np.inf)
x0 = x.copy()
n += 1
return p_print(x)
A = np.array([[5.94, 0.8, 0.6, -3.96, 0.2, 0.3],
[2.97, 6.4, 0, -2.97, 0.2, 0.2],
[2.97, 3.5, 8.7, 1.98, 0.2, 0],
[4.95, 1.6, 1.2, -8.91, 0.8, 0.3],
[-0.99, 2.5, 1.1, -3.96, 9, 0.4],
[5.94, 1.4, 2.4, 0, 3.2, 13]])
b = np.array([11.44, -54.75, -4.64, 20.47, -95.86, 26.92])
print(seidel(A, b, np.zeros(b.shape), 0.00001))