Résolution de systèmes linéaires

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 $x_0$, puis d'éliminer $x_0$ des équations suivantes, ensuite de résoudre la seconde équation en $x_1$, puis d'éliminer $x_1$ des équations suivantes, et ainsi de suite. Initialement, on a le système de $N$ équations à $N$ inconnues :

$$\sum_{j=0}^{N-1} a_{ij} x_j = b_i, \qquad 0 \le i\le N$$

Les fonctions implémentant cet algorithme vont opérer sur les tableaux globaux :

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 $a_{ij}^0 = a_{ij}$ et $b_i^0 = b_i$. Après $k$ itérations, on obtient le système :

$$\renewedcommand {arraycolsep}{1pt}\begin{array}{cccccccccccc... ...\ldots &+ &a_{N-1,N-1}^k x_{N-1} & = &b_{N-1}^k \\ \end{array}$$

Les $k$ premières équations ne seront plus modifiées. Il faut maintenant éliminer $x_k$. On n'utilise pas nécessairement l'équation $a_{kk}^k x_k + \ldots$ car il se peut que le coefficient $a_{kk}^k$ soit nul ; on va donc chercher un indice $l$, avec $k\le l\le N-1$, tel que $a_{lk}^k$ 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 $a_{kk}^k x_k + \ldots$ et $a_{lk}^k x_k + \ldots$ :

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 $k+1$ premières équations ne seront plus modifiées. On a donc désormais $a_{kk}^k \not= 0$. On peut donc éliminer $x_k$ de la ligne $i$ ( $k+1\le i \le N-1$) en ajoutant cette ligne à $-a_{ik}^k/a_{kk}^k$ fois la ligne $k$ : on obtient ainsi les coefficients du système à l'issue de cette $k+1$-ème itération, par << pivotage >> :

$$\begin{array}{rcll} a_{ij}^{k+1} &= &a_{ij}^k - a_{ik}^k a_{... ...{ik}^k b_{k}^k / a_{kk}^k, &\quad k+1\le i \le N-1 \end{array}$$

Ceci se transcrit en la procédure pivoter :

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 $N$), il est transformé en un système triangulaire en $N-1$ 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 $\Theta(N^3)$.