Un exemple : la triangulation de systèmes linéaires

Les tableaux forment la structure de données de base pour de nombreux algorithmes numériques. En voici un exemple classique. 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 109#109, puis d'éliminer 109#109 des équations suivantes, ensuite de résoudre la seconde équation en 110#110, puis d'éliminer 110#110 des équations suivantes, et ainsi de suite. Initialement, on a le système de 111#111 équations à 111#111 inconnues :

112#112

Pour décrire l'algorithme, on posera 113#113 et 114#114. Après 115#115 itérations, on obtient le système :

116#116

Les 115#115 premières équations ne seront plus modifiées. Il faut maintenant éliminer 117#117. On n'utilise pas nécessairement l'équation 118#118 car il se peut que le coefficient 119#119 soit nul ; on va donc chercher un indice 120#120, avec 121#121, tel que 122#122 ait la plus grande valeur absolue (pour minimiser les problèmes d'arrondi). Si la valeur du pivot est non nulle, on échange les équations 118#118 et 123#123. Après cet échange, les 124#124 premières équations ne seront plus modifiées. On a donc désormais 125#125. On peut donc éliminer 117#117 de la ligne 107#107 ( 126#126) en ajoutant cette ligne à 127#127 fois la ligne 115#115 : on obtient ainsi les coefficients du système à l'issue de cette 124#124-ème itération, par « pivotage » :

128#128

Ceci se transcrit en la procédure pivoter. 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 111#111), il est transformé en un système triangulaire supérieur en 129#129 itérations. Il reste à résoudre ce système triangulaire. La complexité de cet algorithme est en 130#130. Voici les fonctions réalisant cet algorithme :

public class Gauss {

  static int pivot(double[][] a, int k) {
    // cherche l'indice d'un pivot dans k..dim-1
    int l = k;
    double max = Math.abs(a[k][k]);

    for (int i=k+1; i<dim; i++) {
      if (Math.abs(a[i][k]) > max) {
        l = i;
        max = Math.abs(a[i][k]);
      }
    }
    return l;
  }

  static void échanger(double[][] a, double[] b, 
                               int k, int l) {
    // échange les lignes k et l du système
    double[] ligne = a[k];
    a[k] = a[l];
    a[l] = ligne;
    double temp = b[k];
    b[k] = b[l];
    b[l] = temp;
  }

  static void pivoter(double[][] a, double[] b, int k) {
    // élimine x_k des lignes k+1..dim-1
    for (int i=k+1; i<dim; i++) {
      double q = a[i][k]/a[k][k];
      b[i] = b[i] - q * b[k];
      for (int j=k+1; j<dim; j++) {
        a[i][j] = a[i][j] - q * a[k][j];
      }
    }
  }

  public static boolean gauss(double[][] a, double[] b) {

    // si le système "a x = b" est régulier, le transforme
    // en un système triangulaire en dim-1 itérations et
    // retourne true, sinon retourne false

    boolean inversible = false;
    for (int k=0; k<dim-1; k++) {
      int l = pivot(k);
      inversible = (Math.abs(a[l][k]) > EPS);
      if (inversible) {
        if (l > k) {
          échanger(k,l);
        }
        pivoter(k);
      } else break;
    }
    return inversible;
  }

  public static double[] solutionTriangulaire(double[][] a, 
                                              double[] b) {

  // calcule la solution x d'un système triangulaire
  // supérieur "a x = b" et retourne x s'il est régulier,
  // et déclenche une exception sinon

    double[] x = new double[dim];
    for (int i=dim-1; i>=0; i--) {
      if (Math.abs(a[i][i])<EPS) {
        throw new ArithmeticException("système irrégulier");
      } else {
        double v = b[i];
        for (int j=i+1; j<dim; j++) {
          v = v - a[i][j]*x[j];
        }
        x[i] = v/a[i][i];
      }
    }
    return x;
  }
}

La fonction solutionTriangulaire ne retourne un tableau solution que si le système est régulier ; sinon, elle déclenche une exception (voir § 5.3). L'utilisation de cet algorithme se fait de la façon suivante :

class GaussTest {
  public static void main(String[] args) {
    double[][] membreGauche = new double[][] { ... };
    double[] membreDroit = new double[] { ... };
    if (Gauss.gauss(membreGauche, membreDroit)) {
      double[] solution = 
        Gauss.solutionTriangulaire(membreGauche, membreDroit);
      // ...
    }
  }
}