Transcript Méthode du pivot de Gauss

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
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def pivot_gauss(A0,Y0):
’’’Algorithme de résolution du système matriciel A0.X = Y0.
* La matrice A0 doit être inversible pour que l’algorithme fonctionne
(système de Cramer) et est donnée sous forme de liste de listes où chaque
sous-liste représente une ligne de la matrice.
* Y0 est une simple liste.
L’algorithme renvoie la valeur attendue de X sous forme d’une liste.
’’’
# Initialisation: on copie les matrices/listes en profondeur pour ne pas
# modifier les données de l’utilisateur
A,Y = [copy.deepcopy(v) for v in [A0,Y0]]
n = len(A)
# Quelques vérifications de base:
assert n == len(A[0])
# A doit être une matrice carrée n*n,
assert n == len(Y)
# Y doit être un vecteur de taille n
for i in range(n):
# On itère sur chacune des n lignes
j = cherche_pivot(A,i)
# Recherche parmi les lignes j >= i
if j > i:
# Si le pivot n’est pas en i, on échange
echange_lignes(A,i,j)
# à la fois pour A
Y[i],Y[j] = Y[j],Y[i]
# et pour Y
for k in range(i+1,n):
# On itère sur les lignes restantes
# On effectue la transvection Lk = Lk - alpha*Li
alpha = A[k][i]/A[i][i]
# avec alpha bien choisi
transvection(A,k,i,alpha) # transvection sur A
Y[k] = Y[k] - alpha*Y[i] # et sur Y
# Il ne reste plus qu’à remonter en partant de la fin.
X = [0]*n
for i in range(n-1,-1,-1):
X[i] = (Y[i] - sum(A[i][j]*X[j] for j in range(i+1,n)))/A[i][i]
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
2
3
4
5
6
7
8
9
10
def cherche_pivot(A,i):
’’’Recherche du meilleur pivot pour la i-eme variable dans les lignes
suivant la ligne courante. Pour éviter d’avoir un pivot presque nul (pb de
comparaison à 0 pour les flottants...), on va prendre le plus grand en
valeur absolue.’’’
i_mem = i # Initialisation de la mémoire
for j in range(i+1,len(A)):
# On regarde les lignes ultérieures
if abs(A[j][i]) > abs(A[i_mem][i]):# Si le pivot de la j-eme ligne est meilleur,
i_mem = j
# on le garde.
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
2
3
4
5
def echange_lignes(A,i,j):
’’’Échange des lignes i et j dans la matrice A.
Ne renvoie rien mais *modifie* la matrice A directement.’’’
for k in range(len(A)):
# On itère sur les éléments de la ligne
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
2
3
4
5
def transvection(A,k,i,alpha):
’’’Transvection de la ligne Lk sous la forme Lk = Lk - alpha*Li.
Ne renvoie rien mais *modifie* la matrice A directement.’’’
for j in range(len(A)):
# On itère sur les éléments de la ligne
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
2
3
4
5
6
7
8
9
10
A = [[1, 1, 1, 1],
[1, 1, 2, 3],
[2, 3, 1,-1],
[3,-1,-1, 1]]
Y = [ 4,
13,
-5,
6]
X = pivot_gauss(A,Y)
print(X)
11
12
import numpy as np
13
14
15
X_np = np.matrix(A)**-1 * np.matrix(Y).transpose()
print(X_np)
16
17
18
19
20
21
22
A = np.random.rand(9,9)
Y = np.random.rand(9)
X2= pivot_gauss(A,Y)
X2_np = np.matrix(A)**-1 * np.matrix(Y).transpose()
print([round(v,4) for v in X2])
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.]]
[-85.9627, 28.4468, -39.9808, -24.3522, 48.3874, 45.2101, 2.1598, 16.0131, 3.416]
[[-85.9627]
[ 28.4468]
[-39.9808]
[-24.3522]
[ 48.3874]
[ 45.2101]
[ 2.1598]
[ 16.0131]
[ 3.416 ]]
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.016 |
0.025 |
0.023 |
0.003 |
4 | 0.008 |
0.013 |
0.017 |
0.018 |
8 | 0.009 |
0.021 |
0.032 |
0.087 |
16 | 0.016 |
0.036 |
0.059 |
0.390 |
32 | 0.031 |
0.067 |
0.111 |
1.659 |
64 | 0.062 |
0.134 |
0.221 |
6.738 |
128 | 0.121 |
0.294 |
0.447 |
27.481 |
256 | 0.243 |
0.518 |
0.855 |
107.868 |
512 | 0.498 |
1.047 |
1.721 |
442.391 |
1024 | 0.977 |
2.088 |
3.422 |
1792.431 |
2048 | 2.272 |
5.860 |
9.618 |
7188.753 |
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& -
INS3
Pivot de Gauss
Code INS3.9: Tests de vitesse globaux
Taille n | Triangulation | Remontee | Pivot complet | Pivot via NumPy |
----------------------------------------------------------------------2 |
0.130 |
0.024 |
0.073 |
0.249 |
4 |
0.178 |
0.039 |
0.208 |
0.148 |
8 |
0.951 |
0.099 |
1.033 |
0.140 |
16 |
6.281 |
0.288 |
6.653 |
0.155 |
32 |
46.058 |
0.963 |
48.535 |
0.310 |
64 |
358.827 |
3.349 |
372.968 |
0.783 |
128 |
2774.608 |
12.749 |
2723.760 |
2.435 |
5/8
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& -
INS3
Pivot de Gauss
Code INS3.10: Même NumPy marche tout de même en n3 ...
Taille n | Pivot via NumPy |
---------------------------2 |
0.212 |
4 |
0.130 |
8 |
0.134 |
16 |
0.142 |
32 |
0.176 |
64 |
0.488 |
128 |
1.840 |
256 |
14.125 |
512 |
65.049 |
1024 |
544.379 |
2048 |
3617.972 |
6/8
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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def inversion(A0):
’’’Algorithme d’inversion matricielle de la matrice A0 par pivot de
Gauss. L’idée est d’appliquer les transformations induite sur A0 à la
matrice identité au fur et à mesure. Quand A0 arrive à l’identité, alors
les transformations initiées sur l’identité donne l’inverse de A0.
* La matrice A0 doit être inversible pour que l’algorithme fonctionne et
est donnée sous forme de liste de listes où chaque sous-liste représente
une ligne de la matrice. Le format np.array() fonctionne aussi.
NB: la matrice doit être carrée.
L’algorithme renvoie L’inverse de la matrice sous forme de np.array()’’’
# Initialisation: on copie la matrice A0 en profondeur pour ne pas
# modifier les données de l’utilisateur
A = copy.deepcopy(A0)
n = len(A)
# Quelques vérifications de base:
assert n == len(A[0])
# A doit être une matrice carrée n*n,
# Initialisation de l’inverse sous forme de la matrice identité
A_inv = np.eye(n)
# renvoie automatiquement un np.array()
for i in range(n):
# On itère sur chacune des n lignes
j = cherche_pivot(A,i)
# Recherche parmi les lignes j >= i
if j > i:
# Si le pivot n’est pas en i, on échange
echange_lignes(A,i,j)
# à la fois pour A
echange_lignes(A_inv,i,j) # et pour son "inverse to be"
for k in range(i+1,n):
# On itère sur les lignes restantes
# On effectue la transvection Lk = Lk - alpha*Li
alpha = A[k][i]/A[i][i]
# avec alpha bien choisi
transvection(A,k,i,alpha) # transvection sur A
transvection(A_inv,k,i,alpha) # et sur l’inverse to be
# Il ne reste plus qu’à remonter en partant de la fin.
for i in range(n-1,-1,-1):
a = A[i][i]
ligne = A[i][:]
substitution(A,i,a,ligne)
# Substitution sur A
substitution(A_inv,i,a,ligne) # puis sur son "inverse to be"
return A_inv
35
36
37
38
39
40
41
42
43
44
def substitution(A,i,a,ligne):
’’’Substitution vue comme une multiplication de ligne puis une série de
transvections avec les lignes suivantes pour virer les termes non diagonaux’’’
n = len(A)
# La substitution se fait en deux phases:
for j in range(n):
A[i][j] /= a
# d’abord un passage Li = Li/a
for k in range(n-1,i,-1):
# ensuite une suite de transvections
transvection(A,i,k,ligne[k]/a) # Li = Li - Lik/a * Lk
45
46
47
48
A = [[1,2,3],[4,5,6],[1/2,3,4]]
print(inversion(A))
# Résultat de notre algorithme
print(np.matrix(A)**(-1)) # Résultat by NumPy
E Bougnol, JJ Fleck,
M Heckmann & M Kostyra,
Kléber, PCSI& -
INS3
Pivot de Gauss
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]]
8/8