#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Sun Nov 26 13:03:32 2017

@author: pascal
"""

import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt

N=100
#Calcul de la solution de référence
d0=np.ones(N)
d1=np.ones(N-1)
A=4*np.diag(d0,0)+np.diag(d1,1)+np.diag(d1,-1)
b=np.ones(N)
x_exact=npl.solve(A,b)

#Méthode de Jacobi
x0=np.zeros(N) #donnée initiale
x=x0
tol=1e-7 #niveau de résolution de Ax=b
test=1
compteur=0
X_jacobi=[npl.norm(x0-x_exact)] #vecteur des erreurs pour la méthode de Jacobi
while test>tol and compteur<200:
    x0p=np.hstack((x0[1:N],0))
    x0m=np.hstack((0,x0[0:N-1]))
    x=(b-x0p-x0m)/4
    test=max(npl.norm(x-x0), npl.norm(np.dot(A,x)-b))
    x0=x
    X_jacobi=np.hstack((X_jacobi,npl.norm(x0-x_exact)))
    compteur+=1
    
plt.plot(np.log(X_jacobi))   
##On obtient une droite: la convergence est linéaire.

#Méthode de Gauss Seidel
x0=np.zeros(N)
x=x0
tol=1e-7
test=1
compteur=0
X_gauss=[npl.norm(x0-x_exact)]
while test>tol and compteur<200:
    x[0]=(b[0]-x0[1])/4
    for j in np.arange(1,N-1,1):
        x[j]=(b[j]-x[j-1]-x0[j+1])/4
    x[N-1]=(b[N-1]-x[N-1])/4
    test=max(npl.norm(x-x0), npl.norm(np.dot(A,x)-b))
    x0=x
    X_gauss=np.hstack((X_gauss,npl.norm(x0-x_exact)))
    compteur+=1
    
plt.plot(np.log(X_gauss),'r')    
#On obtient une droite de pente supérieure à la pente de Jacobi: la vitesse de
#convergence est plus rapide.