INS3
Pivot de Gauss



Code INS3.1: Implémentation de la fonction principale pour le pivot de Gauss

1 import copy # pour la copie profonde
2

3 def pivot_gauss(A0,Y0):
4 ’’’Algorithme de résolution du système matriciel A0.X = Y0.
5 * La matrice A0 doit être inversible pour que l’algorithme fonctionne
6 (système de Cramer) et est donnée sous forme de liste de listes où chaque
7 sous-liste représente une ligne de la matrice.
8 * Y0 est une simple liste.
9 L’algorithme renvoie la valeur attendue de X sous forme d’une liste.
10 ’’’
11 # Initialisation: on copie les matrices/listes en profondeur pour ne pas
12 # modifier les données de l’utilisateur
13 A,Y = [copy.deepcopy(v) for v in [A0,Y0]]
14 n = len(A) # Quelques vérifications de base:
15 assert n == len(A[0]) # A doit être une matrice carrée n*n,
16 assert n == len(Y) # Y doit être un vecteur de taille n
17 for i in range(n): # On itère sur chacune des n lignes
18 j = cherche_pivot(A,i) # Recherche parmi les lignes j >= i
19 if j > i: # Si le pivot n’est pas en i, on échange
20 echange_lignes(A,i,j) # à la fois pour A
21 Y[i],Y[j] = Y[j],Y[i] # et pour Y
22 for k in range(i+1,n): # On itère sur les lignes restantes
23 # On effectue la transvection Lk = Lk - alpha*Li
24 alpha = A[k][i]/A[i][i] # avec alpha bien choisi
25 transvection(A,k,i,alpha) # transvection sur A
26 Y[k] = Y[k] - alpha*Y[i] # et sur Y
27 # Il ne reste plus qu’à remonter en partant de la fin.
28 X = [0]*n
29 for i in range(n-1,-1,-1):
30 X[i] = (Y[i] - sum(A[i][j]*X[j] for j in range(i+1,n)))/A[i][i]
31 return X
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 2/8



Code INS3.2: Algorithme général du pivot de Gauss

# Initialisation
# Quelques vérifications de base, puis
# On itère sur chacune des n lignes via le compteur i
# On cherche le meilleur pivot pour la i-eme variable
# parmi les lignes j >= i
# Si le pivot n’est pas en ligne i, on échange les deux lignes.
# On itère alors sur les lignes restantes via le compteur k
# On effectue la transvection Lk = Lk - alpha*Li
# avec alpha bien choisi
# Il ne reste plus qu’à remonter en partant de la fin.



Code INS3.3: Implémentation de la recherche du pivot à partir de la ie ligne

1 def cherche_pivot(A,i):
2 ’’’Recherche du meilleur pivot pour la i-eme variable dans les lignes
3 suivant la ligne courante. Pour éviter d’avoir un pivot presque nul (pb de
4 comparaison à 0 pour les flottants...), on va prendre le plus grand en
5 valeur absolue.’’’
6 i_mem = i # Initialisation de la mémoire
7 for j in range(i+1,len(A)): # On regarde les lignes ultérieures
8 if abs(A[j][i]) > abs(A[i_mem][i]):# Si le pivot de la j-eme ligne est meilleur,
9 i_mem = j # on le garde.
10 return i_mem # On renvoie le numéro de la ligne du pivot.



Code INS3.4: Implémentation de l’échange des lignes i et j d’une matrice.

1 def echange_lignes(A,i,j):
2 ’’’Échange des lignes i et j dans la matrice A.
3 Ne renvoie rien mais *modifie* la matrice A directement.’’’
4 for k in range(len(A)): # On itère sur les éléments de la ligne
5 A[i][k],A[j][k] = A[j][k],A[i][k] # et on fait l’échange



Code INS3.5: Implémentation de la procédure de transvection Lk ← Lk − αLi

1 def transvection(A,k,i,alpha):
2 ’’’Transvection de la ligne Lk sous la forme Lk = Lk - alpha*Li.
3 Ne renvoie rien mais *modifie* la matrice A directement.’’’
4 for j in range(len(A)): # On itère sur les éléments de la ligne
5 A[k][j] = A[k][j] - alpha*A[i][j] # et on effectue la transvection
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 3/8




Code INS3.6: Vérification sur l’exemple initial

1 A = [[1, 1, 1, 1],
2 [1, 1, 2, 3],
3 [2, 3, 1,-1],
4 [3,-1,-1, 1]]
5 Y = [ 4,
6 13,
7 -5,
8 6]
9 X = pivot_gauss(A,Y)
10 print(X)
11

12 import numpy as np
13

14 X_np = np.matrix(A)**-1 * np.matrix(Y).transpose()
15 print(X_np)
16

17 A = np.random.rand(9,9)
18 Y = np.random.rand(9)
19 X2= pivot_gauss(A,Y)
20 X2_np = np.matrix(A)**-1 * np.matrix(Y).transpose()
21 print([round(v,4) for v in X2])
22 print(np.around(X2_np,decimals=4))




Code INS3.7: Résultat du test précédent

[-1.0000000000000167, 3.000000000000041, -5.0000000000000595, 7.000000000000032]
[[-1.]
[ 3.]
[-5.]
[ 7.]]
[0.5642, -1.5298, -0.0893, -0.83, -0.41, 1.2085, 1.3018, 0.9811, 0.3951]
[[ 0.5642]
[-1.5298]
[-0.0893]
[-0.83 ]
[-0.41 ]
[ 1.2085]
[ 1.3018]
[ 0.9811]
[ 0.3951]]
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 4/8




Code INS3.8: Tests de rapidité

Taille n | pivot | echange | transvection simple | transvections a partir de n/2 |
----------------------------------------------------------------------------------
2 | 0.006 | 0.010 | 0.009 | 0.002 |
4 | 0.006 | 0.005 | 0.006 | 0.007 |
8 | 0.004 | 0.009 | 0.009 | 0.027 |
16 | 0.008 | 0.017 | 0.017 | 0.113 |
32 | 0.013 | 0.032 | 0.032 | 0.439 |
64 | 0.022 | 0.056 | 0.056 | 1.702 |
128 | 0.043 | 0.110 | 0.111 | 6.060 |
256 | 0.085 | 0.183 | 0.198 | 23.405 |
512 | 0.152 | 0.357 | 0.359 | 110.292 |
1024 | 0.350 | 0.911 | 1.389 | 414.461 |
2048 | 0.867 | 2.061 | 1.947 | 1693.602 |
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 5/8




Code INS3.9: Tests de vitesse globaux

Taille n | Triangulation | Remontee | Pivot complet | Pivot via NumPy |
-----------------------------------------------------------------------
2 | 0.074 | 0.013 | 0.037 | 0.129 |
4 | 0.056 | 0.012 | 0.059 | 0.063 |
8 | 0.251 | 0.025 | 0.272 | 0.070 |
16 | 1.669 | 0.065 | 1.620 | 0.075 |
32 | 11.870 | 0.242 | 11.742 | 0.119 |
64 | 88.498 | 0.777 | 90.760 | 0.674 |
128 | 700.960 | 2.921 | 699.308 | 2.113 |
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 6/8




Code INS3.10: Même NumPy marche tout de même en n3 ...

Taille n | Pivot via NumPy |
----------------------------
2 | 0.102 |
4 | 0.056 |
8 | 0.066 |
16 | 0.061 |
32 | 0.091 |
64 | 0.377 |
128 | 1.648 |
256 | 9.230 |
512 | 39.846 |
1024 | 232.507 |
2048 | 1556.124 |
4096 | 10732.557 |
8192 | 95080.232 |
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 7/8


Code INS3.11: Aménagement du pivot de Gauss pour inverser la matrice A

1 def inversion(A0):
2 ’’’Algorithme d’inversion matricielle de la matrice A0 par pivot de
3 Gauss. L’idée est d’appliquer les transformations induite sur A0 à la
4 matrice identité au fur et à mesure. Quand A0 arrive à l’identité, alors
5 les transformations initiées sur l’identité donne l’inverse de A0.
6 * La matrice A0 doit être inversible pour que l’algorithme fonctionne et
7 est donnée sous forme de liste de listes où chaque sous-liste représente
8 une ligne de la matrice. Le format np.array() fonctionne aussi.
9 NB: la matrice doit être carrée.
10 L’algorithme renvoie L’inverse de la matrice sous forme de np.array()’’’
11 # Initialisation: on copie la matrice A0 en profondeur pour ne pas
12 # modifier les données de l’utilisateur
13 A = copy.deepcopy(A0)
14 n = len(A) # Quelques vérifications de base:
15 assert n == len(A[0]) # A doit être une matrice carrée n*n,
16 # Initialisation de l’inverse sous forme de la matrice identité
17 A_inv = np.eye(n) # renvoie automatiquement un np.array()
18 for i in range(n): # On itère sur chacune des n lignes
19 j = cherche_pivot(A,i) # Recherche parmi les lignes j >= i
20 if j > i: # Si le pivot n’est pas en i, on échange
21 echange_lignes(A,i,j) # à la fois pour A
22 echange_lignes(A_inv,i,j) # et pour son "inverse to be"
23 for k in range(i+1,n): # On itère sur les lignes restantes
24 # On effectue la transvection Lk = Lk - alpha*Li
25 alpha = A[k][i]/A[i][i] # avec alpha bien choisi
26 transvection(A,k,i,alpha) # transvection sur A
27 transvection(A_inv,k,i,alpha) # et sur l’inverse to be
28 # Il ne reste plus qu’à remonter en partant de la fin.
29 for i in range(n-1,-1,-1):
30 a = A[i][i]
31 ligne = A[i][:]
32 substitution(A,i,a,ligne) # Substitution sur A
33 substitution(A_inv,i,a,ligne) # puis sur son "inverse to be"
34 return A_inv
35

36 def substitution(A,i,a,ligne):
37 ’’’Substitution vue comme une multiplication de ligne puis une série de
38 transvections avec les lignes suivantes pour virer les termes non diagonaux’’’
39 n = len(A)
40 # La substitution se fait en deux phases:
41 for j in range(n):
42 A[i][j] /= a # d’abord un passage Li = Li/a
43 for k in range(n-1,i,-1): # ensuite une suite de transvections
44 transvection(A,i,k,ligne[k]/a) # Li = Li - Lik/a * Lk
45

46 A = [[1,2,3],[4,5,6],[1/2,3,4]]
47 print(inversion(A)) # Résultat de notre algorithme
48 print(np.matrix(A)**(-1)) # Résultat by NumPy
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& - INS3 Pivot de Gauss 8/8




Code INS3.12: Résultat du test d’inversion de matrice: ça marche !

[[ 0.44444444 0.22222222 -0.66666667]
[-2.88888889 0.55555556 1.33333333]
[ 2.11111111 -0.44444444 -0.66666667]]
[[ 0.44444444 0.22222222 -0.66666667]
[-2.88888889 0.55555556 1.33333333]
[ 2.11111111 -0.44444444 -0.66666667]]