Transcript TP 1 : Pivot de Gauss - IC Lycée Paul Constans

PT
IC – TP N O 1 : P IVOT DE G AUSS
1 octobre 2014
I. Quelques rappels sur la théorie
1) Introduction
Pour résoudre un système linéaire de n équations à p inconnues, on se ramène à un problème matriciel : AX = B avec
A ∈ Mn,p (R), X ∈ Mp,1 (R) et B ∈ Mn,1 (R).
Dans ce TP, on traite uniquement la résolution des systèmes de Cramer, c’est à dire tels que la matrice A est une matrice
carrée inversible. Dans ce cas, le système admet une unique solution. Nous nous intéresserons à l’algorithme correspondant à la méthode de Gauss avec recherche partielle du pivot. Cet algorithme permet d’exécuter des opérations
répétitives mais simples sur les lignes de A et B.
Nous étudierons la mise en œuvre de cet algorithme pour mettre en évidence les problèmes sous-jacents notamment
liés à la représentation machine des nombres réels. Enfin, nous détaillerons la complexité de cet algorithme après avoir
identifié la taille du problème.
2) Opérations élémentaires sur les lignes d’un système
On appelle opérations élémentaires sur les lignes d’un système les opérations suivantes :
• La permutation, codée Li ↔ Lk : échange des lignes i et k ;
• La dilatation, codée Li ← αLi : multiplication de la ligne i par un scalaire α non nul ;
• La transvection, codée Li ← Li + αLk : ajout d’un multiple de coefficient α quelconque d’une ligne k à une ligne i .
On dira que deux systèmes (S) et (S 0 ) sont équivalents lorsque l’on peut passer de l’un à l’autre par une suite finie d’opérations élémentaires sur les lignes. On a : (S) ⇔ (S 0 ). Ils ont alors le même ensemble de solutions.
3) Algorithme du pivot de Gauss-Jordan
Soit A ∈ GLn (R). On note a i j le coefficient situé ligne i , colonne j et C j la colonne j de A.
On appelle matrice augmentée du système la matrice à n lignes et n + 1 colonnes telle que la dernière colonne est le
second membre du système. On note A0 cette matrice et on lui applique l’algorithme suivant :
1. On procède au besoin à une permutation de lignes dans la matrice pour que a 11 soit non nul (c’est possible car, A
étant inversible, C1 est non nulle).
2. On annule alors les a i 1 pour 2 É i É n grâce à des transvections.
3. On recommence le procédé sur les colonnes suivantes : pour j ∈ ‚2; nƒ, permutation pour avoir un pivot a j j non nul,
puis transvection pour annuler les coefficients a i j pour j + 1 É i É n. On obtient ainsi une matrice de la forme :


∗ ... ... ∗ ∗

.. .. 
0 . . .
. .


.
.. .. 

 . .. ..
.
. . .
.
0 ... 0 ∗ ∗
4. On attaque ensuite la « remontée »pour obtenir une matrice pseudo diagonale : on annule les coefficients a 1n à a n−1,n
par transvection en utilisant a nn comme pivot, puis on recommence de même sur les colonnes n −1 à 2 pour obtenir
une matrice de la forme :


∗ 0 ... 0 ∗

.. .. 
0 . . .
. .


.
.. 
 . .. ..

.
. 0 .
.
0 ... 0 ∗ ∗
5. Il reste à obtenir des 1 sur la diagonales par dilatation : la dernière colonne donne alors la solution du système.
Remarque : En pratique sous Python, on ne « construira » pas la matrice A0 , mais on réalisera exactement les mêmes
opérations sur le vecteur B que sur la matrice A.
1/3
PT
IC – TP N O 1 : P IVOT DE G AUSS
1 octobre 2014
II. Mise en œuvre sous Python
1) Type de variable
Même s’il serait beaucoup plus pratique d’utiliser des array numpy, on n’utilisera ici que des listes Python. Une matrice
de Mn,p (R) sera donc stockée sous la forme d’une liste de n listes de longueur
 p.
 3x + y + z + t = 252


 3x − 3y + 2z + 3t = 210
Pour la suite du TP, on pourra par exemple chercher à résoudre le système
.
1

x + 34 y + 12
z + 13 t = 231

2

 3x − 3y + 3z + t = 168
Dans l’éditeur, créer un nouveau script Python et le sauver sous un nom bien choisi, dans un répertoire itou. Définir les
variables A et B correspondant au système précédent.
2) Les fonctions « élémentaires »
Implémenter trois fonctions correspondant aux opérations élémentaires vues au I.2)
On aura besoin de plus d’une quatrième fonction permettant de trouver l’indice de la ligne du pivot de la colonne k
(pour k ∈ ‚0; n − 2ƒ : attention aux indices en Python), c’est à dire un entier i ∈ ‚k, n − 1ƒ tel que le coefficient a i k soit non
nul.
3) L’algorithme lui-même
En utilisant les fonctions précédentes, écrire une fonction resolution(A0, B0) qui implémente l’algorithme décrit au
I.3). Les deux premières instructions de la fonction consisteront à faire une copie profonde des listes A0 et B0. Pour cela,
utiliser la fonction deepcopy du module copy.
4) La méthode du pivot partiel
On vu que que le travail avec des flottants entraîne souvent des erreurs d’arrondi (calculer par exemple (12*(1/3-1/4)-1).
Ici, le risque est de faire apparaître des coefficients très petits, mais non nuls, censés représenter le réel nul. Si on prend
ces coefficients comme pivot, on multipliera par leur inverse qui produira de très grands flottants n’ayant aucun sens,
ce qui amènera à un résultat n’ayant plus aucune signification.
 1
1
 4x + 3y +z = 3
0.3x + 0.4y + 2z = 4.4 en utilisant la fonction que vous venez d’implémenter. Que
Par exemple, résoudre le système

x + 2y + 3z = 13
pensez-vous du résultat ?
Une méthode consiste à choisir dans la colonne en cours le plus grand coefficient en valeur absolue. Prendre ce terme
comme pivot n’empêche pas les erreurs, mais est assez efficace dans la pratique. C’est la méthode du pivot partiel.
Écrire une fonction cherche_pivot_p(M, k) qui l’implémente, puis modifier la fonction resolution pour l’utiliser.
Reprendre l’exemple précédent. Le résultat est-il plus satisfaisant ?
III. Complexité
1) Calcul
On va maintenant évaluer la complexité de cet algorithme, c’est à dire le nombre d’opérations élémentaires (addition,
soustraction ,multiplication, division, comparaison) effectuées en fonction de la taille n de la matrice A.
Pour cela, on commence par évaluer la complexité des fonctions utilisées :
• La permutation :
• La dilatation :
• La transvection :
• La recherche du pivot :
2/3
PT
IC – TP N O 1 : P IVOT DE G AUSS
1 octobre 2014
Or, pour la mise sous forme triangulaire, le nombre d’appel de ces fonctions est au pire de :
• recherche de pivot :
• permutation :
• transvection :
• dilatation :
Le nombre total d’opérations élémentaires est ainsi de :
Remarque : on pouvait de façon beaucoup plus rapide dire « pour chaque i , la complexité de la recherche du pivot puis
de l’échange est linéaire. Chaque transvection a une complexité linéaire, et il y en a au plus n, donc on a un complexité
quadratique à i fixé. Comme i varie entre 0 et n − 2, on a une complexité en O(n 3 ) ».
De même, pour la remontée, la complexité est
et enfin le calcul des solutions a une complexité
Au final, la résolution par la méthode du pivot de Gauss a une complexité en O(n 3 ).
2) Vérification
En utilisant la fonction time() du module time, tester votre fonction pour des matrices de taille n ∈ {25; 50; 100; 200; 400}.
Pour créer un système de taille n, on pourra utiliser la fonction random() du module random :
from random import random
A = [[random() for col in range(n)] for ligne in range(n)]
B = [[random()] for ligne in range(n)]
L’évolution du temps d’exécution est-elle conforme aux prévisions ?
3) Comparaison avec numpy
Comparer avec les temps de résolution de la fonction solve() du sous module linalg du module numpy.
IV. Pour ceux qui ont fini
1) Inversion
Écrire une fonction qui calcule l’inverse d’une matrice (supposée inversible).
Indication : commencer par faire un copier-coller des fonctions que l’on vient de taper permettra de n’avoir quasiment
rien à faire.
2) Optimisation
Essayer d’optimiser vos fonctions. Quelques pistes :
• utiliser des matrix et des array du module numpy pour remplacer les listes ;
• utiliser le fait qu’au moment où on s’occupe de la ligne k, tous les indices a i j pour k + 1 É i É n et 1 É j É k − 1 sont
nuls. Cela peut permettre d’améliorer les fonctions permutation et transvection.
This constructive proof works for any m × n matrix !
3/3