import numpy as np
import matplotlib.pyplot as plt
n =30
eps =1E-5
A = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        if i == j:
            A[i, j] = 2
        elif abs(i - j) == 1:
            A[i, j] = -1
D = np.diag(np.diag(A))
E = -np.triu(A) + D
F = -np.tril(A) + D
uexact = np.ones(n)
b = np.dot(A,uexact)
uinit = np.random.rand(n)
errinit = np.max(np.abs(uinit - uexact))
# Jacobi
M1 = np.linalg.inv(D)
N1 = D - A
u = uinit
niter=1
errj = errinit
errj_list = [errinit]
while errj > eps:
    niter=niter+1
    u = np.dot(M1, np.dot(N1, u) + b)
    errj = np.max(np.abs(u - uexact))
    errj_list.append(errj)
print('nombre d''iterations Jacobi') 
print(niter)
plt.figure(1)
plt.plot(range(niter), np.log(errj_list), color='black')

# Gauss Seidel
M2 = np.linalg.inv(D - E)
N2 = F
u = uinit
niter=1
errgs = errinit
errgs_list = [errinit]
while errgs > eps:
    niter=niter+1
    u = np.dot(M2, np.dot(N2, u) + b)
    errgs = np.max(np.abs(u - uexact))
    errgs_list.append(errgs)
print('nombre d''iterations Gauss Seidel') 
print(niter)   
plt.figure(1)
plt.plot(range(niter), np.log(errgs_list), color='blue')
plt.show()