# -*- coding: utf-8 -*-
import numpy as np
import numpy.linalg as npl

#Exercice 1

#question 2: Définition de la matrice H

H=np.array([[1, 1/2. ,1/3. ,1/4. ,1/5.],[1/2., 1/3., 1/4., 1/5., 1/6.], [1/3., 1/4., 1/5., 1/6., 1/7.], [1/4., 1/5., 1/6., 1/7., 1/8.],[1/5., 1/6., 1/7., 1/8., 1/9.]])
 
#QUESTION 3: Calcul des normes subordonnées de H (utilisation de la librairie numpy.linalg)

#Calcul de la norme subordonnée de H à la norme 1
print 'La norme de H subordonnée à la norme 1 vaut', npl.norm(H, ord=1)

#Calcul de la norme subordonnée de H à la norme 1
print 'La norme de H subordonnée à la norme 2 vaut', npl.norm(H, ord=2)

# Calcul de la norme subordonnée de H à la norme infinie
print 'La norme de H subordonnée à la norme infinie vaut', npl.norm(H, ord=np.inf)

#QUESTION 4:  Calcul des conditionnements de la matrice H

#Calcul de la norme subordonnée de H à la norme 1
print 'Le conditionnement de H subordonné à la norme 1 vaut', npl.cond(H, 1)

#Calcul de la norme subordonnée de H à la norme 1
print 'Le conditionnement de H subordonné à la norme 2 vaut', npl.cond(H, 2)

# Calcul de la norme subordonnée de H à la norme infinie
print 'La conditionnement de H subordonné à la norme infinie vaut', npl.cond(H,np.inf)

#QUESTION 5: Calcul du déterminant

print 'Le déterminant de H vaut', npl.det(H)

#QUESTION 6: Calcul des valeurs propres et des vecteurs propres de H

vap,vep= npl.eig(H)
print 'Les valeurs propres de H sont', vap
print 'Les vecteurs propres de H sont:'
print vep
#Remarque: pour vérifier que le premier vecteur colonne de la matrice vep est la vecteur propre associé à vap[0], on peut taper la commande:
# np.dot(H,vep[0:5,0])-vap[0]*vep([0:5,0])

#QUESTION 7: résolution de Hx_0=b

b=np.array([1,1,1,1,1])
x0=npl.solve(H,b)
print x0

#QUESTION 8: résolution de Hx1=b+delta b

delta_b=np.array([0,0,0,0, 0.01])
x1=npl.solve(H,b+delta_b)
print x1

#calcul de l'erreur relative (pour la norme 1)

err1=npl.norm(x1-x0,1)/npl.norm(x0,1)

print err1

#Explication du resultat obtenu: on part d'une erreur sur le second membre de 1% et on fait une erreur de 80% sur le résultat: 
#ceci est dû au mauvais conditionnement de la matrice H (le conditionnement de la matrice peut être vu comme un facteur d'amplification
#de l'erreur commise sur les données que sont la matrice H et le vecteur b.


# Exercice 2: résolution (directe) d'un système linéaire

# QUESTION 2: création de la matrice A
def A(n):
    return 4*np.diag(np.ones(n),0)+np.diag(np.ones(n-1),1)+np.diag(np.ones(n-1),-1)

# Affichage de la matrice A pour n=10        
print A(10)    

#QUESTION 3: factorisation LU de A(10)

#La décomposition LU n'est pas implémentée dans Numpy: il faut passer à SciPy

#IMPORT des librairies SCIPY et SCIPY.LINALG
import scipy.linalg as spl
  
L,U=spl.lu(A(10),True) #Décomposition de la 

print L
print U   

#Remarque: vérification que A(10)=LU 

print 'La norme de A-LU est égale à', npl.norm(np.dot(L,U)-A(10))  

# QUESTION 4
# Dans le calcul de la décomposition LU, on n'a besoin (comme pour le pivot de Gauss) que d'annuler un seul terme à chaque étape du calcul.
# Comme la matrice est de taille n, on fait donc au plus n opérations.

#QUESTION 5
# Calcul de X0
X0=npl.solve(A(100),np.ones(100))  
                                                              
#QUESTION 6
#définition du vecteur b+delta b
z=np.hstack((1.01,np.ones(99)))
X1=npl.solve(A(100),z)

# Calcul de l'erreur relative

err2=npl.norm(X1-X0,1)/npl.norm(X0,1)

print 'L erreur relative vaut', err2

#QUESTION 7

errb=0.01/npl.norm(np.ones(100),1)

est_cond=err2/errb

print 'Le conditionnement de la matrice A(100) est minoré par', est_cond

#QUESTION 8

print 'Le condition de la matrice A(100) donné par Python est', npl.cond(A(100),1)


# QUESTION 9: factorisation de Cholesky de la matrice A pour n=4.

B=npl.cholesky(A(4))

print B

#Remarque: vérification que A(4)=B*B^T: on calcule la norme de la différence

#print npl.norm(np.dot(L,L.T)-A(4))