des zéros dans gauss

des zéros dans gauss - Algo - Programmation

Marsh Posté le 16-01-2006 à 23:13:17    

Salut à tous,
 
j'ai écrit un petit algorithme de Gauss-Jordan pour l'inversion matricielle.
le pb est qu'il y a des risques de divisions par zéros
par exemple, lorsque l'on veut normaliser à cette étape :
 
1 2 3 4 | 3 2 7 5
0 1 5 3 | 4 9 5 1
0 0 0 3 | 2 5 7 4
0 0 5 2 | 2 4 9 4
 
j'ai donc insérer une recherche sur le maximum (en valeur absolue)... pour qu'ici, par exemple, les 2 dernières lignes soient inversées
je n'ai trouvé ça nulle part sur internet, ce qui m'étonne un peu. Est-ce qu'on peut éventuellement faire autrement, et est-ce que c'est utile (d'autant que ça peut prendre un certain temps pour les grandes matrices) ?
 
code correspondant, sorti de tout contexte (mais testé) :
nbp est le nombre de paramètres
deriv2 est la matrice à inverser (matrice essienne nbp x nbp)
derinv est la matrice que deviendra, à la fin de l'algorithme, la matrice inverse (nbp x nbp également)
 

Code :
  1. // recherche du coef max
  2. max = 0;
  3. for (j = i ; j <= nbp ; j++) {
  4. temp = abs(deriv2[j][i]);
  5. if (temp > max) {lmax = j ; max = temp;}
  6. } // fin de recherche du coef max
  7. // echange des 2 lignes
  8. for (j = i ; j <= nbp ; j++) {
  9. temp = deriv2[i][j];
  10. deriv2[i][j] = deriv2[lmax][j];
  11. deriv2[lmax][j] = temp;
  12. }
  13. for (j = 1 ; j <= nbp ; j++) {
  14. temp = derinv[i][j];
  15. derinv[i][j] = derinv[lmax][j];
  16. derinv[lmax][j] = temp;
  17. } // fin d'echange des 2 lignes


 
P.S. : je viens de me rendre compte que je n'ai mis aucune condition sur l'échange de ligne, qui est particulièrement inutile si il n'y a rien à échanger... [:albator] ... j'ai en plus bien conscience que je prend de la mémoire inutilement en n'utilisant pas l'indice 0 des tableaux :D


Message édité par jimipage le 18-01-2006 à 22:47:36

---------------
un perlien qui programme salement
Reply

Marsh Posté le 16-01-2006 à 23:13:17   

Reply

Marsh Posté le 17-01-2006 à 15:48:56    

pour ceux que ça intéresse ou qui auraient des remarques, voici l'algorithme complet :
 

Code :
  1. // DEBUT DE LA SUBROUTINE GAUSS
  2. int Gauss(unsigned int nbp , double deriv2[][nbp+1] , double derinv[][nbp+1]) {
  3. unsigned int i , j , k , lmax;
  4. double max , temp;
  6. // initialisation de l'inverse (= partie droite de la matrice etendue)
  7. for (i = 1 ; i <= nbp ; i++) {
  8. for (j = 1 ; j <= nbp ; j++) {
  9. if (i != j) derinv[i][j] = 0; else derinv[i][j] = 1;
  10. } }
  12. // boucle principale de triangulation de la partie gauche de la matrice etendue
  13. for (i = 1 ; i <= nbp ; i++) {
  14. // recherche du coef max
  15. max = 0;
  16. for (j = i ; j <= nbp ; j++) {
  17. temp = abs(deriv2[j][i]);
  18. if (temp > max) {lmax = j ; max = temp;}
  19. } // fin de recherche du coef max
  20. // echange des 2 lignes
  21. for (j = i ; j <= nbp ; j++) {
  22. temp = deriv2[i][j];
  23. deriv2[i][j] = deriv2[lmax][j];
  24. deriv2[lmax][j] = temp;
  25. }
  26. for (j = 1 ; j <= nbp ; j++) {
  27. temp = derinv[i][j];
  28. derinv[i][j] = derinv[lmax][j];
  29. derinv[lmax][j] = temp;
  30. } // fin d'echange des 2 lignes
  31. // normalisation
  32. temp = deriv2[i][i];
  33. for (j = i ; j <= nbp ; j++) {deriv2[i][j] /= temp;}
  34. for (j = 1 ; j <= nbp ; j++) {derinv[i][j] /= temp;}
  35. // fin de normalisation
  36. // reduction
  37. for (j = i+1 ; j <= nbp ; j++) {
  38. temp =  deriv2[j][i]; // combien de fois il faut retrancher la ligne normalisee
  39. for (k = i+1 ; k <= nbp ; k++) {deriv2[j][k] -= temp * deriv2[i][k];}
  40. for (k = 1 ; k <= nbp ; k++) {derinv[j][k] -= temp * derinv[i][k];}
  41. } // fin de reduction
  42. } // fin de boucle principale de triangulation
  45. // normalisation de la derniere ligne
  46. temp = deriv2[nbp][nbp];
  47. deriv2[nbp][nbp] = 1;
  48. for (i = 1 ; i <= nbp ; i++) {derinv[nbp][i] /= temp;}
  50. // boucle principale de reduction
  51. for (i = nbp-1 ; i >= 1 ; i--) {
  52. for (j = i+1 ; j <= nbp ; j++) { // boucle sur les colonnes d'une ligne donnee
  53. temp = deriv2[i][j]; // combien de fois il faut retrancher la ligne j
  54. for (k = 1 ; k <= nbp ; k++) {derinv[i][k] -= temp * derinv[j][k];}
  55. } // fin boucle sur une ligne donnee
  56. } // fin de boucle principale de reduction
  58. }
  59. // FIN DE LA SUBROUTINE GAUSS


Message édité par jimipage le 17-01-2006 à 17:10:36

---------------
un perlien qui programme salement
Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed