Transcript Le pivot de Gauss

Ingénierie numérique TP no 4
3 semaines
Le pivot de Gauss
MPSI
20142015
Objectifs : L'objectif de ce TP est de savoir programmer l'algorithme du pivot de Gauss pour
résoudre des systèmes linéaires, pour calculer l'inverse d'une matrice carré inversible et comparer le
temps de calcul de notre algorithme par rapport à ceux fournis par la bibliothèque numpy. Pour ceux
qui auront le temps, ils pourront, en plus, programmer la décomposition LU et le calcul de rang par
l'algorithme de Gauss.
Exercice No 1 : Dans un sous-dossier adapté, créer un sous-dossier consacré à ce TP.
1
Résolution d'un système linéaire
Rappelons que pour résoudre le système Ax = y à l'aide du pivot de Gauss, on eectue l'algorithme
suivant :
Pour i de 0 à n − 2 faire
Trouver j entre i et n − 1 tel que |aj,i | soit maximal. # Cela permet de trouver un terme
#non nul et de minimiser les erreurs d'arrondis.
Échanger Li et Lj . # Coecients de la matrice et membre de droite
Pour k de i + 1 à n − 1 faire
Lk ← Lk −
Fin Pour
Fin Pour
ak,i
ai,i Li
Arrivé ici, le système est sous forme triangulaire et il n'y a plus qu'à remonter , via des substitutions.
Le résultat est mis dans un vecteur x et il s'agit donc de calculer :
1
xi =
ai,i
yi −
n−1
X
!
ai,k xk
.
k=i+1
On procède alors de la manière suivante :
Pour i de n − 1 à 0 faire
Pour k de i + 1 à n − 1 faire
yi ← yi − ai,k xk
Fin Pour
xi ←
yi
ai,i
Fin Pour
Exercice No 2 : Créer un nouveau chier intitulé systeme_lineaire, modier l'encodage UTF-8 en
1
latin1 et importer la bibliothèque numpy. Ce chier prendra en entrée A et y sous forme de liste et
retournera sous forme de liste, l'élément x vériant Ax = y (sous réserve d'existence).
Exercice No 3 :
Dans le chier systeme_lineaire, créer la procédure
echange_ligne(i,j)
qui interverti les lignes i et j du système linéaire (on prendra garde à bien intervertir les lignes de A
et de y ).
Exercice No 4 :
Dans le chier systeme_lineaire, créer la procédure
transvection(i,j,mu)
qui eectue l'opération élémentaire Li ← Li + µ × Lj .
Exercice No 5 :
Dans le chier systeme_lineaire, créer la procédure
dilatation(i,mu)
qui eectue l'opération élémentaire Li ← µ × Li .
Exercice No 6 :
Dans le chier systeme_lineaire, créer la fonction
pivot_partiel(A,j0)
qui retourne en sortie le i > j0 tel que |ai,j0 | soit maximal.
Exercice No 7 :
Compléter le chier systeme_lineaire pour qu'il prenne entrée A et y sous forme de liste et retourne
sous forme de liste, l'élément x vériant Ax = y (on supposera que x existe et qu'il est unique).
Exercice No 8 :
Tester votre programme systeme_lineaire sur les exemples suivants :
2x + 3y = 5
5x − 2y = −16

 2x +2y −3z = 2
y
−6z = −3

z
=4

+2y −3z = 2
 2x
−2x −y −3z = −5

6x
+4y +4z = 16
Comparer les résultats avec ceux obtenus à l'aide de la fonction linsolve du package numpy.linalg.
On rappelle que pour utiliser numpy.linalg.solve(A,y) il faut que A et y soient de type array.
2
Inversion de matrice
L'algorithme du pivot de Gauss nous a amené à faire des opérations sur les lignes d'une matrice,
ce qui revient à multiplier à gauche cette matrice par des matrices élémentaires, disons L1 , · · · , LN
(après N opérations). Le système Ax = y est alors équivalent à L1 Ax = L1 y puis L2 L1 Ax = L2 L1 y ,
puis LN · · · L1 Ax = LN · · · L1 y et si on a correctement appliqué l'algorithme du pivot de Gauss, le
2
membre de droite vaut exactement x. Ainsi A−1 = LN · · · L1 . Cette matrice est en fait exactement
celle obtenue en appliquant à la matrice identité (et dans le même ordre) les opérations réalisées sur
A. On peut donc appliquer l'algorithme suivant pour calculer l'inverse d'une matrice A :
X ← In
Pour i de 0 à n − 2 faire
Trouver j entre i et n − 1 tel que |aj,i | soit maximale.
Échanger Li et Lj . # Pour la matrice A et la matrice X .
Pour k de i + 1 à n − 1 faire
a
Li # Pour la matrice A et la matrice X .
Lk ← Lk − ak,i
i,i
Fin Pour
Fin Pour
Pour i de n − 1 à 0 faire
Pour j de i − 1 à 0 faire
a
Lj ← Lj − aj,i
Li # Pour la matrice A et la matrice X .
i,i
Fin Pour
Fin Pour
Pour i de 0 à n − 1 faire
Li ← a1i,i Li # Pour la matrice A et la matrice X .
Fin Pour
Exercice No 9 :
Créer un nouveau chier intitulé inversion qui prend en entrée A et retourne en sortie A−1 sous
forme de matrice. Il s'agit donc d'écrire en langage Python, l'algorithme précédent. Contrairement
au programme précédent, on prendra garde à ne pas modier la matrice A donnée en entrée. On
s'inspirera bien évidement de la méthode eectuée pour le programme systeme_lineaire.
Exercice No 10 :
Tester votre programme inversion sur les exemples suivants :


2 2 −3
0 1 −6
0 0 1
1 2
3 5


2
2 −3
−2 −1 −3
6
4
4
Comparer les résultats avec ceux obtenus à l'aide de la fonction inv du package numpy.linalg. On
rappelle que pour utiliser numpy.linalg.inv(A) il faut que A soit de type array.
3
Comparaison des temps de calcul
La matrice de Virginie d'ordre n est la matrice V ∈ Mn (R) suivante :


2 −1
(0)
−1 2 −1





.
.
.
.
.
.


.
.
.



−1 2 −1
(0)
−1 2
3
Exercice No 11 : Créer un nouveau chier intitulé temps_calcul. Ce programme prendra en entrée
un entier n et retournera en sortie le temps de calcul de V (n)−1 par notre algorithme et par la fonction
inv.
Exercice No 12 :
Écrire le programme temps_calcul. On utilisera les mots-clefs et la syntaxe suivante :
import timeit
t0=timeit.timeit()
· · · [Mon calcul] · · ·
t1=timeit.timeit() # t1-t0 est le temps (en secondes) du calcul
Exercice No 13 :
Faire tourner le programme précédent pour n ∈ {50, 100, 200, 400}. Conclure.
4
La décomposition
LU
Exercice No 14 :
Écrire une fonction decomposition(A) prenant en entrée une matrice carré A inversible admettant
une décomposition LU et retournant en sortie les deux matrices L et U de cette décomposition. On
utilisera l'algorithme donné dans le cours.
Exercice No 15 :
En utilisant la fonction précédente, écrire un programme decomposition_systeme_lineraire qui
prend en entrée A et Y et retourne en sortie la solution X de AX = Y .
5
Calcul de rang
Un théorème de mathématiques dit la chose suivante :
Théorème 1
Soient A ∈ Mn,p (K) et r ∈ N tel que 0 6 r 6 min(n, p). On a équivalence entre :
(i) rg(A) = r,
(ii) ∃P ∈ GLp (K), ∃Q ∈ GLn (K) telles que A = QJr P .
Exercice No 16 :
Écrire un programme intitulé rang qui prend en entrée une matrice A et retourne en sortie les données
r, Q et P du théorème précédent 1 .
1. Si vous arrivez à faire cet exercice, c'est que je n'ai plus rien à vous apprendre sur le pivot de Gauss.
4