next up previous contents index
Next: La fonction main Up: Tableaux d'objets Previous: Sous-typage

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

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

Pour décrire l'algorithme, on posera aij0 = aij et bi0 = bi. Après k itérations, on obtient le système :


\begin{displaymath}\renewcommand{\arraycolsep}{1pt}
\begin{array}{ccccccccccccc}...
...\ldots &+
&a_{N-1,N-1}^k x_{N-1} & = &b_{N-1}^k \\
\end{array}\end{displaymath}

Les k premières équations ne seront plus modifiées. Il faut maintenant éliminer xk. On n'utilise pas nécessairement l'équation $a_{kk}^k x_k + \ldots$ car il se peut que le coefficient akkksoit nul ; on va donc chercher un indice l, avec $k\le l\le N-1$, tel que alkk 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 $a_{kk}^k x_k + \ldots$ et $a_{lk}^k x_k + \ldots$. 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 xk de la ligne i ( $k+1\le i \le N-1$) en ajoutant cette ligne à -aikk/akkk fois la ligne k : on obtient ainsi les coefficients du système à l'issue de cette k+1-ème itération, par << pivotage >> :

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

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 N), il est transformé en un système triangulaire supérieur en N-1 itérations. Il reste à résoudre ce système triangulaire. La complexité de cet algorithme est en $\Theta(N^3)$. Voici les fonctions implémentant cet algorithme :

class SystemeLineaire {

  private 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;
  }

  private 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;
  }

  private 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];
      }
    }
  }

  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;
  }

  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;
  }
}


next up previous contents index
Next: La fonction main Up: Tableaux d'objets Previous: Sous-typage
R. Lalement
2000-10-23