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;
  5. // initialisation de l'inverse (= partie droite de la matrice etendue)
  6. for (i = 1 ; i <= nbp ; i++) {
  7. for (j = 1 ; j <= nbp ; j++) {
  8. if (i != j) derinv[i][j] = 0; else derinv[i][j] = 1;
  9. } }
  10. // boucle principale de triangulation de la partie gauche de la matrice etendue
  11. for (i = 1 ; i <= nbp ; i++) {
  12. // recherche du coef max
  13. max = 0;
  14. for (j = i ; j <= nbp ; j++) {
  15. temp = abs(deriv2[j][i]);
  16. if (temp > max) {lmax = j ; max = temp;}
  17. } // fin de recherche du coef max
  18. // echange des 2 lignes
  19. for (j = i ; j <= nbp ; j++) {
  20. temp = deriv2[i][j];
  21. deriv2[i][j] = deriv2[lmax][j];
  22. deriv2[lmax][j] = temp;
  23. }
  24. for (j = 1 ; j <= nbp ; j++) {
  25. temp = derinv[i][j];
  26. derinv[i][j] = derinv[lmax][j];
  27. derinv[lmax][j] = temp;
  28. } // fin d'echange des 2 lignes
  29. // normalisation
  30. temp = deriv2[i][i];
  31. for (j = i ; j <= nbp ; j++) {deriv2[i][j] /= temp;}
  32. for (j = 1 ; j <= nbp ; j++) {derinv[i][j] /= temp;}
  33. // fin de normalisation
  34. // reduction
  35. for (j = i+1 ; j <= nbp ; j++) {
  36. temp =  deriv2[j][i]; // combien de fois il faut retrancher la ligne normalisee
  37. for (k = i+1 ; k <= nbp ; k++) {deriv2[j][k] -= temp * deriv2[i][k];}
  38. for (k = 1 ; k <= nbp ; k++) {derinv[j][k] -= temp * derinv[i][k];}
  39. } // fin de reduction
  40. } // fin de boucle principale de triangulation
  41. // normalisation de la derniere ligne
  42. temp = deriv2[nbp][nbp];
  43. deriv2[nbp][nbp] = 1;
  44. for (i = 1 ; i <= nbp ; i++) {derinv[nbp][i] /= temp;}
  45. // boucle principale de reduction
  46. for (i = nbp-1 ; i >= 1 ; i--) {
  47. for (j = i+1 ; j <= nbp ; j++) { // boucle sur les colonnes d'une ligne donnee
  48. temp = deriv2[i][j]; // combien de fois il faut retrancher la ligne j
  49. for (k = 1 ; k <= nbp ; k++) {derinv[i][k] -= temp * derinv[j][k];}
  50. } // fin boucle sur une ligne donnee
  51. } // fin de boucle principale de reduction
  52. }
  53. // 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