CPGE OUJDA SPE

Résolution d'un système linéaire inversible à l'aide de la méthode de Gauss

1. Introduction

Il existe deux types de méthodes de résolution d’un système linéaire :

Ø Résolution dite directe à l’aide du pivot de Gauss, que nous allons étudier

Ø Les méthodes itératives (ou indirectes) : on part d’un vecteur et on considère une suite récurrente du type :

On approximera la solution par une valeur de la suite qui converge vers la solution.

2. Description de l’algorithme du pivot de Gauss

Dans ce texte, on suppose que les systèmes linéaires sont de Cramer, c’est-à-dire admettent une unique solution. La matrice A est donc inversible. On suppose que A est de taille n de coefficients et que attentioncomme en Python, les indices commencent à 0.

L’algorithme du pivot de Gauss, consiste à trianguler la matrice A à l’aide d’opérations élémentaires.

a. Le cas triangulaire

Si la matrice A est triangulaire, on résoud facilement le système par "substitutions remontantes".

Le système s’écrit :

On a ainsi :

On en déduit le pseudo-code suivant pour résoudre un système triangulaire :

Algo: triangle
Données: A matrice triangulaire supérieure inversible et b un vecteur
Résultat: solution de Ax =b
n = longueur de b
x = un tableau de longueur n
x[n-1] = b[n-1]/ a[n-1][n-1]
Pour i de n-2 à 0 (pas de -1)
s = 0
Pour j de i+1 à n-1
s = s + a[i][j]*x[j]
Fin pour
x[i] = (b[i] -s)/a[i][i]
Fin Pour

b. Mise sous forme triangulaire, phase d’élimination

A chaque étape, on obtient un sous-système carré qui a une inconnue de moins.

Un exemple modèle tracé en Python.

Le pseudo-code :
Algo: triangulation
Données: A matrice inversible et b un vecteur
Résultat: mise sous forme triangulaire du système Ax =b
n = longueur de b
Pour i de 0 à n-2
chercher premier indice i0 entre i et n-1 tel que A[i0][i] non nul # pivot
échanger les lignes d’indice i0 et i pour A et b
Pour k de i+1 à n-1
L_k<- L_k - L_i0 *A[k][i0]/A[i0][i0] # pour A et b
Fin Pour
Fin Pour

Quelques justifications :

Ø Comme A est inversible, on est sûr de trouver un pivot non nul à chaque itération.

Ø À chaque étape k, pour tout i compris entre et , la variable a disparu de toutes les équations du système à partir de la . C’est l’invariant de boucle.

Remarque :

Ø Cette méthode d’élimination est extrêmement importante et s’adapte à des matrices rectangulaires, on peut par exemple alors récupérer le rang.

Ø Certaines matrices ne nécessitent pas de recherche de pivot et donc pas d’échanges de lignes, c’est le cas par exemple des matrices à diagonales strictement dominante ou des matrices définies positives.

c. La problématique supplémentaire des flottants

Quel type de données doit-on choisir pour représenter les coefficients du système ? On ne peut pas choisir le type entier puisque des divisions interviennent, on peut ainsi prendre le type flottant. Deux problématiques apparaissent alors :

Ø La comparaison à zéro, lorsque que l’on recherche un pivot non nul. On rappelle qu’on ne teste jamais l’égalité à zéro d’un flottant.

Ø Les erreurs arrondis, un exemple pour comprendre : Supposons que les nombres flottants soient codés en base 10 avec une mantisse de 3 chiffres et considérons le système suivant :

Son unique solution est le couple (x, y) = (1.0001, 0.9999).

En choisissant le nombre 10⁻⁴ comme pivot, en effectuant la transvection, le système est équivalent à :

Mais et comme la mantisse ne possède que 3 chiffres, ce nombre est arrondi à De même − est arrondi à

Ainsi est arrondi à 1 et par suite donc On a ici une perte de précision importante.

En divisant par un petit pivot, on a obtenu de grands nombres.

Résolvons maintenant d’une autre façon : on permute les lignes 1 et 2 et on choisit 1 comme pivot. On effectue ensuite la transvectionet le système est équivalent à:

Mais est arrondi à est aussi arrondi à Donc puis en remontant. Cette fois-ci le résultat est satisfaisant.

Morale : il vaut mieux perdre des chiffres significatifs pour un nombre petit que pour un nombre grand.

C’est pourquoi, on préférera choisir comme pivot "le plus grand des pivots". On utilisera donc l’amélioration suivante où l’on va chercher un pivot maximal. On parle de pivot partiel (on parle de pivot total lorsque le pivot est maximal en ligne et en colonne, mais on ne le fera pas ici, car on ne raisonne que sur les lignes).

Remarque : on peut éviter les problématiques de nombres flottants si les coefficients des matrices sont des rationnels ou des éléments de corps finis. On peut alors représenter les coefficients de manière exacte en utiliser le type fraction de Python par exemple ou un logiciel de calcul formel. Mais alors, pourquoi ne choisit-on pas toujours les rationnels au lieu des flottants ? Et bien cela dépend des besoins. Avec des rationnels, on calcule de manière exacte mais lentement et avec des flottants, on calcule très vite mais de manière approchée.

3. Implémentation en Python

a. Découpage du travail

Codons tout d’abord les deux opérations élémentaires suivantes :

defechange_ligne(A,i,j):
#Applique l’opération élémentaire Li <->Lj à la matrice A
n = len(A[0]) # nbre de colonnes de A
for k inrange(n):
temp = A[i][k]
A[i][k] = A[j][k]
A[j][k] = temp

Attention cette version ne fonctionne pas pour une matrice colonne car n’est plus une liste mais un entier. Cela fonctionne si on écrit comme une liste de liste au lieu de

defechange_ligne_bis(A,i,j):
#Applique l’opération élémentaire Li <->Lj à la matrice A
temp = A[i]
A[i] = A[j]
A[j] = temp

def transvection(A,i,j,mu):
# Applique l’opération élémentaire Li <-Li + mu*Lj à la matrice A
# cas où A est une matrice colonne
if type(A[0]) != list:
A[i] = A[i] + mu*A[j]
else:
n = len(A[0]) # nbre de colonnes de A
for k in range(n):
A[i][k] = A[i][k] + mu*A[j][k]

Remarque : ces deux procédures suivantes ne renvoient rien, elles modifient la matrice passée en paramètre.
Codons maintenant la fonction qui recherche le pivot maximal.

defpivot_partiel(A, j0):
n = len(A) # nbre de lignes de A
imax = j0 # indice de ligne avec pivot max
for i inrange(j0+1, n):
if abs(A[i][j0]) > abs(A[imax][j0]):
imax = i
returnimax

Écrivons enfin une fonction triangle qui résoud les systèmes triangulaires supérieurs.

deftriangle(A,b):
""" Renvoie la solution du système Ax =b lorsque A est tringulaire supérieure inversible"""
n =len(b)
x = [0 for i inrange(n)]
x[n-1] = b[n-1]/A[n-1][n-1]
for i inrange(n-2,-1,-1):
s = 0
for j in range(i+1,n):
s = s + A[i][j]*x[j]
x[i] = (b[i] - s)/ A[i][i]
return x

b. Recoller les morceaux

Puisque nous allons appliquer des échanges et des transvections à la matrice de départ, celle-ci va être modifiée. Nous allons donc faire en premier lieu une copie de notre donnée. Nous allons pour cela utiliser la fonction copy.deepcopy. Rappelons que si on pense réaliser une copie de A0 en effectuant l’affectation A =A0, il n’en est rien puisque A et A0 sont deux noms correspondant à un même objet.

import copy
defsolution_systeme(A0,b0):
""" Renvoie la solution X du système de Cramer A0 X= b0"""
A,b = copy.deepcopy(A0), copy.deepcopy(b0)

# on ne détruit pas les données initiales
n = len(A)
# Mise sous forme triangulaire
for i inrange(n-1): # on itère n-1 fois
i0 = pivot_partiel(A,i)
echange_ligne_bis(A,i,i0)
echange_ligne_bis(b,i,i0)
for k inrange(i+1,n):
x = A[k][i]/A[i][i]
transvection(A,k,i,-x)
transvection(b,k,i,-x) # attention b matrice colonne
# phase de remontée
returntriangle(A,b)

Testons avec :

>>> A = [[2,2,-3], [-2,-1,-3], [6, 4, 4]]

>>> b = [2,-5,16]

>>>solution_systeme(A,b)
[-14.000000000000036, 21.000000000000046, 4.000000000000007]

c. Comparaison avec Numpy et un logiciel de calcul formel

La fonction numpy.linalg.solve permet de résoudre les systèmes linéaires. Elle utilise des tableaux de type array.
>>> x = numpy.linalg.solve(A,b)
array([-14., 21., 4.])
>>> x[0]
-14.000000000000036