# -*- coding: utf-8 -*-
"""
Created on Tue Dec 23 16:17:33 2014

@author: grivet
Listing 6.4 : Résolution d'un système d'équations 
linéaires par la méthode itérative de Jacobi

"""
import numpy as np

def affiche_mat(titre,M):
    print (titre)
    for row in M:
        for val in row:
            print (' %10.5f' %val, end='')
        print()

seuil = 10**(-6)
N = 5;
A = np.random.rand(N,N)-0.5
A = A + A.transpose()+2.0*np.eye(N,N)
affiche_mat('A = ',A)
print()

b = range(N)
x = np.ones(N); y = np.ones(N)
kmax = 20
for k in range(kmax):
    for i in range(N):
        s = b[i] - np.dot(A[i,:],x)
        y[i] = x[i] + s/A[i,i]
    if max(abs(y-x)) < seuil:
        break
    x = np.copy(y)
print ('nombre d iterations: ',k)
print ('solution: ', x)
print ('vecteur d erreurs')
print (np.dot(A,x) - b)