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 :
Les k premières équations ne seront plus modifiées. Il faut
maintenant éliminer xk. On n'utilise pas nécessairement l'équation
car il se peut que le coefficient akkksoit nul ; on va donc chercher un indice l, avec
,
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
et
.
Après cet échange, les k+1 premières équations ne seront plus
modifiées. On a donc désormais
.
On peut donc
éliminer xk de la ligne i (
)
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 >> :
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; } }