La méthode de Gauss transforme un système linéaire en un système
triangulaire équivalent. L'idée est de résoudre la première équation en
, puis d'éliminer
des équations suivantes, ensuite de
résoudre la seconde équation en
, puis d'éliminer
des
équations suivantes, et ainsi de suite. Initialement, on a le système
de
équations à
inconnues :
const int N = 10; double a[N][N]; // la matrice double b[N]; // le membre droit double x[N]; // la solution
Pour décrire l'algorithme, on posera
et
. Après
itérations, on obtient le système :
Les premières équations ne seront plus modifiées. Il faut
maintenant éliminer
. On n'utilise pas nécessairement l'équation
car il se peut que le coefficient
soit nul ; on va donc chercher un indice
, avec
, tel
que
ait la plus grande valeur absolue (pour minimiser les
problèmes d'arrondi) :
int pivot(int k) { // cherche l'indice d'un pivot dans k..N-1 int i; int l = k; double max = fabs(a[k][k]); for (i=k+1; i<N; i++) { if (fabs(a[i][k]) > max) { l = i; max = fabs(a[i][k]); } } cout << "pivot(" << k << ") =" << l << endl; return l; }
Si la valeur du pivot est non nulle, on échange les équations
et
:
void echanger(int k, int l) { // échange les lignes k et l du système // à partir de la colonne k double temp; int j; for (j=k; j<N; j++) { temp = a[k][j]; a[k][j] = a[l][j]; a[l][j] = temp; } temp = b[k]; b[k] = b[l]; b[l] = temp; }
Après cet échange, les premières équations ne seront plus
modifiées. On a donc désormais
. On peut donc
éliminer
de la ligne
(
) en ajoutant cette
ligne à
fois la ligne
: on obtient ainsi les
coefficients du système à l'issue de cette
-ème itération, par <<
pivotage >> :
void pivoter(int k) { // élimine x_k des lignes k+1..N-1 int i,j; for (i=k+1; i<N; i++) { double q = a[i][k]/a[k][k]; b[i] = b[i] - q * b[k]; for (j=k+1; j<N; j++) { a[i][j] = a[i][j] - q * a[k][j]; } } }
Si, à l'une des itérations, la valeur du pivot est nulle, c'est que le
système n'est pas régulier, et la résolution s'arrête (c'est le rôle du
break dans la fonction
gauss). Si le système est régulier (c'est-à-dire de rang ),
il est transformé en un système triangulaire en
itérations.
bool gauss() { // si le système "a x = b" est régulier, le transforme en un // système triangulaire en N-1 itérations et retourne true, // sinon retourne false int k, l; bool inversible; for (k=0; k<N-1; k++) { l = pivot(k); inversible = (fabs(a[l][k]) > DBL_EPSILON); if (inversible) { if (l > k) { echanger(k,l); } pivoter(k); } else break; } return inversible; }
Il reste à résoudre un système triangulaire supérieur :
bool systeme_triangulaire() { // calcule la solution x d'un système triangulaire supérieur // "a x = b" et retourne true s'il est régulier, // et retourne false sinon int i, j; double v; for (i=N-1; i>=0; i--) { if (fabs(a[i][i])<EPS) { return false; } else { v = b[i]; for (j=i+1; j<N; j++) { v = v - a[i][j]*x[j]; } x[i] = v/a[i][i]; } } return true; }
La complexité de cet algorithme est en .