[Matrices] Comment calculer le Déterminant d'une matrice 4x4

Version imprimable

Voir 40 message(s) de cette discussion en une page

25/07/2005, 17h40
cyber_N

[Matrices] Comment calculer le Déterminant d'une matrice 4x4

HELLO LE MONDE!! :D

Je souhaiterais savoir comment calculer le déterminant d'une matrice 4x4
(ET NON D'UNE MATRICE 3X3 ou 2x2)

Avez des explications et/ou des liens intéressants ??

Merci encore.
Séb.
:wink: :wink:
25/07/2005, 18h23
Wavyx

tu peux commencer par là:
http://mathworld.wolfram.com/Determinant.html

mais il existe pe des methodes spécifiques pour les 4x4 ??

sinon ben ya la méthode générale :roll:
25/07/2005, 18h37
alveric

Tu as besoin de performances optimisees, ou tu veux juste un truc qui marche ? C'est une matrice entiere, reelle, complexe ?
Si tu veux un truc qui ne marche que pour les 4x4, mais qui le fait bien... bah tape la formule complete ! :wink:
Ou alors tu prends la formule 7 du lien que Wavyz t'a donne, et tu retombes sur 4 determinants de degre 3, pour lesquels la formule est donnee quelques lignes plus bas.
Ou alors tu fais une triangulation de la matrice, et tu fais le produit de la diagonale grace a la formule 24.

Il y a le choix... Mais je ne me souviens pas de methode particuliere pour les matrice d'ordre 4...
27/07/2005, 13h52
Matthieu Brucher

La triangularisation à l'aidede matrices de Houseolder est sans doute une des meilleures méthodes que je connaisse pour calculer un déterminant - complexité en O(n^3) -
28/07/2005, 08h59
Invité

Bonjour, je fais remonter ce topic, car moi aussi j'aurais besoin d'un algo permettant de calculer le déterminant d'une matrice. Mais il me faudrait une formule générale, car l'ordre de mes matrices va varier de 2 à 6 en général.
J'ai bien entendu regarder le lien de Wavyx, mais il n'est pas évident d'en déduire la forumle générale.
J'ai bien compris qu'il fallait faire :
Somme(i=1 à k) des a(i,j)C(i,j) où C(i,j) = (-1)^(1+j)M(i,j)
mais il faut aussi calculer M(i,j)

Bref, ma date de fin de projet approche à grands pas, et j'ai pas énormément de temps (hélas) à passer à pondre de nouveaux algos, je voudrais donc savoir s'il n'y a pas un algo existant quelque part, voir même une implémentation dans un language connu (comme du C par exemple) qui me permettrait de déduire l'algo

Merci beaucoup :)
28/07/2005, 10h12
Wavyx

en gros il suffit de faire un appel récursif jusqu'à une dimension 3 par exemple.
sinon tu peux regarder dans matlab comment ils ont fait ou encore en cherchant un peu sur le web et sur source forge:
http://sourceforge.net/projects/jlinalg/
http://sourceforge.net/projects/slate/
http://sourceforge.net/projects/ezmatrixcalc/

mais je m'y connais pas en librairie pour les matrices ou l'algèbre linéaire donc pê que qqun aura de meilleurs liens/sources ;)
28/07/2005, 10h16
mathieu_t

Franchement la méthode donnée par Wavx est de mon point de vue la plus simple à implémenter, d'autant que ça prend la forme d'un algo récursif (on décompose le déterminant d'une matrice n*n par des déterminants de matrice (n-1)*(n-1) etc ...)

Sinon j'ai trouvé en tapant sous Google "class Matrix C++" le lien suivant : http://www.techsoftpl.com/matrix/download.htm mais il y en a bien d'autres...

Cela calcule le déterminant apparemment...

A+
28/07/2005, 10h49
Invité

en effet, j'avais pas pense a regarder du cote de sourceforge :)

je vais essayer d'eplucher ca rapidement ;)

mais c'est vrai qu'au final, je sens qu'il sera plus rapide de chercher l'algo puis de l'implementer, je n'ai pas non plus envie de passer par une librairie externe (je sais, je suis ch... difficile) :P

merci pour votre aide :)
28/07/2005, 10h55
alveric

Le probleme avec la formule utilisant les mineurs d'ordre n-1, c'est qu'il faut appliquer l'algo recursivement sur les mineurs. Or, ces mineurs sont des portions pour la plupart non-contigues de la matrice d'origine, ce qui signifie qu'il faut reconstruire, a chaque niveau d'appel, une matrice n-1 en recopiant les valeurs voulues...
Du coup, on se retrouve a recopier (n-1)^2 valeurs dans la matrice qui sert de mineur, on applique l'algo sur la nouvelle matrice pour connaitre le determinant, on multiplie si necessaire par -1, et on fait ca n fois !
Autrement dit, C(n) = n * C(n-1) + n*((n-1)^2 copies de valeur) + (n-1 additions) + (n/2 multiplications par -1)
A moins qu'il y ait une optimisation possible pour les copies de valeurs, ou un moyen efficace de travailler directement sur la matrice d'origine, ce dont je doute (mais je peux me tromper), ca devient vite lourd.
De plus, ca signifie qu'il faut avoir en memoire une matrice d'ordre n, une d'ordre n-1... jusqu'a 1.
Bon, pour n <= 6, ca peut etre acceptable, apres tout je n'ai pas fait les tests, mais la triangularisation serait plus efficace a priori (enfin, a mon priori).
28/07/2005, 10h59
Invité

je suis bien conscient du probleme en effet, mais je ne connais pas d'autre méthode :(

peux tu m'en dire plus sur la triangularisation ?
28/07/2005, 11h21
alveric
Cette methode utilise la propriete de multilinearite du determinant (16 du premier lien de Waxyz), et le fait que le determinant d'une matrice triangulaire est le produit des termes de la diagonale (24).
Le but est d'annuler toute la partie sous la diagonale, par des combinaison lineaire de lignes.
ex:
Code:

1 2 3 |a11 a12 a13| |a21 a22 a23| |a31 a32 a33|
A la premiere iteration, tu annules a21 en ajoutant -(a21/a11) fois la premiere ligne a la deuxieme. Tu obtiens
Code:

1 2 3 |a11 a12 a13| | 0 b22 b23| |a31 a32 a33|
Puis tu recommences pour annuler a31, puis a32. Tu obtiens alors
Code:

1 2 3 |a11 a12 a13| | 0 b22 b23| | 0 0 c33|
Et le determinant vaut

Code:

det(A) = a11 * b22 * c33

Evidemment, si a11 est nul, il faut permuter des lignes... (vu que le determinant est aussi invariant par permutation de lignes)
Je ne sais plus la complexite de l'algo, mais ca doit se retrouver.
28/07/2005, 11h31
Invité

ok, c'est peut etre en effet plus rapide :)

aurais tu un bout d'algo ou de code de cette methode ? :roll: :P
28/07/2005, 11h43
alveric

Je te l'ai deja presque donne...
Tu commences par la premiere ligne. Si a11 est nul, tu cherches sur la premiere colonne un terme non nul. Si tu en trouves un, tu permutes les deux lignes et passes a la suite. Sinon le determinant est nul et tu peux t'arreter.
Une fois que tu as un a11 non nul, tu cherches a annuler a21, puis a31...an1. (un truc: calcule 1/a11 une seule fois dans une variable temporaire. Les divisions sont toujours des calculs longs pour le CPU, alors autant en faire le moins possible)
Et tu recommences avec b22...
La condition d'arret, c'est soit arriver a un akk nul et tous les ajk , j>k nuls aussi, soit avoir triangule totalement la matrice.
28/07/2005, 12h01
Invité

oui, je sais, je suis deja parti sur l'algo, mais ca aurait bien pratique si quelqu'un l'avait deja fait, car bcp plus rapide :P

je suis desole, je suis un peu pris par le temps :(
28/07/2005, 12h40
mathieu_t

Citation:

mais la triangularisation serait plus efficace a priori

Oui en effet, j'ai regardé rapido un algo que j'ai, ça passe par la triangularisation...

Je ne peux pas poster la classe Matrix que j'ai (code trop long), mais je peux l'envoyer par mails à ceux qui le souhaitent...

En revanche dans le lien que j'ai donné au dessus, ils le font par le pivot (méthode récursive donnée au dessus)...

Je peux quand même donner les bouts intéressants :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// determinant(u, w, v) -- Returns the determinant of a square matrix,
// which is the product of the diagonal weights in W, O(n).
// Determinant of a rectangular matrix is not well defined, 0 is returned.
// Elimination is done on the matrix and singularies are marked away
// if W is empty, otherwise the decomposition cached in (U, W, V) is used.
 
double Matrix::determinant(Matrix& u, Matrix& w, Matrix& v) const {
  if (this->num_rows != this->num_cols) {
    cerr << "Return 0 for determinant of a rectangular matrix." << endl;
    return 0;
  }
  if (w.data == NULL)
    this->singularities(u, w, v);               // cache SVD in U,W,V.
  double product = 1;
  for (unsigned long k = 0; k < w.num_cols; k++)		// product of diagonals in W.
    product *= w(k, k);				// small weights should be zeroed
  return product;				// out, so that determinant = 0.
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
// singularities -- Returns the number of singular weights which make
// the well-condition number less than given lower bound threshold. O(n).
// Default lower bound on the well-condition number is 1.0e-6.
// These singular weights are marked by setting them to zero.
// The inverse of these zero weights will be set to zero instead of
// infinite, which corresponds to finding least-norm solution.
 
int Matrix::singularities(Matrix& u, Matrix& w, Matrix& v, double lbound) const {
  if (w.data == NULL)
    this->singular_value_decomposition(u, w, v); // cache SVD in U,W,V
  double max_w = 0;
  for (unsigned long k=0; k < w.num_cols; k++) {// find max absolute weight
    double weight = (double) fabs(w(k, k));		// in diagonal
    if (max_w < weight) max_w = weight;
  }
  double lbound_w = (double) (max_w * lbound);	// double must support "* double"
  unsigned long count = 0;
  for (unsigned long l=0; l < w.num_cols; l++) {// find near-singular weights
    double weight = (double) fabs(w(l, l));		// and zero them out
    if (weight < lbound_w) {
      w(l, l) = 0;
      count++;					// count the number of singular
    }						// diagonals.
  }
  return count;
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// svd(u, w, v) -- Computes the singular value decomposition of THIS matrix, and
// returns the 3 matrices U, W, V, such that: *this = U * W * V.transpose(). O(n^3).
// Small diagonals in W should be zeroed out, to clearly mark singularities
// in the matrix. The range space are columns of U such that the
// corresponding diagonals in W are non zero. The null space are columns
// of V such that the corresponding diagonals in W are zero.
// The vectors in U and V are orthonormal. If the matrix is symmetric, then
// U and V are equal and they contain the eigenvectors,
// W contains the eigenvalues.
// The eigenvalues are positive if the matrix is positive definite.
// The matrices U, W, V are returned for the user to mark out singularities
// based on some lower bound on the well-condition number, such as 1.0e-6.
// These matrices are cached so that subsequent computations such
// as rank, determinant, condition number, inverse, back substitution, will
// not need to compute the expensive SVD again.
// See Numerical Recipees for reference.
 
bool Matrix::singular_value_decomposition(Matrix& u, Matrix& w, Matrix& v) const {
    unsigned long n = this->num_cols;			 // dimension of vectors in V.
    unsigned long m = (n>this->num_rows)? n : this->num_rows; // dimension of vectors in U.
    u.resize(m, n);				 // fill *this into u.
    u.update(*this);
    for (unsigned long i = this->num_rows; i < m; i++)	// clear the extra rows
      for (unsigned long j = 0; j < n; j++)
	u.data[i][j] = 0;
    w.resize(n, n);
    w.fill(0);					// clear w
    v.resize(n, n);
    if (Matrix::svdcmp(u.data, m, n,	// NR does all work.
				 w.data[0], // and stores results
				 v.data)) { // in place
      for (unsigned long k = 1; k < n; k++) {
	w(k, k) = w(0, k);			// fill diagonal of w
	w(0, k) = 0;
      }
      return true;
    } else
      return false;
}

Voilà, en espérant que ça aidera...

A+

28/07/2005, 15h21
Invité

très intéressant merci :)

je t'ai envoyé mon mail en MP, pour avoir la suite si possible :roll:

voila ce que j'ai fait a present :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
    //----- parcours de la matrice pour verrifier que son determinant n'est pas nul -----
    Indice = 1
     Répéter
        //on cherche d'abord sur la ligne
        LigneNulle = Vrai
        //on se place sur la ligne voulu
        Ligne = Indice
        //et on parcourt chaque colonne
        Colonne = 1
         Répéter
            //des qu'on tombe sur un coefficient different de zero
            Si (TabCoeffs(Ligne, Colonne) <> 0) Alors
                //on dit que la ligne n'est pas nulle
                LigneNulle = Faux
             FinSi
            //on passe a la colonne suivante
            Colonne = Colonne + 1
        //on sort s'il n'y a plus de colonnes
        //ou si on sait que la ligne n'est pas nulle
         Jusqu 'à ((Colonne > Dimension) Ou (LigneNulle = Faux))
 
        //si la ligne est bonne, on regarde la colonne
        //par le meme principe
        Si (LigneNulle = Faux) Alors
            ColonneNulle = Vrai
            Colonne = Indice
            Ligne = 1
             Répéter
                Si (TabCoeffs(Ligne, Colonne) <> 0) Alors
                    ColonneNulle = Faux
                 FinSi
                Ligne = Ligne + 1
             Jusqu 'à ((Ligne > Dimension) Ou (ColonneNulle = Faux))
         FinSi
 
        //une fois qu'on a regarde ligne et colonne,
        //on passe a l'indice suivant
        Indice = Indice + 1
    //on s'arrete si on a parcouru toute la matrice
    //ou si on a trouve une ligne ou une colonne nulle
     Jusqu 'à ((Indice > Dimension) Ou (LigneNulle = Vrai) Ou (ColonneNulle = Vrai))
 
    //si une ligne ou une colonne est nulle, le determinant est nul
    Si (LigneNulle = Vrai) Ou (ColonneNulle = Vrai) Alors
        DeterminantCalcule = 0
    //sinon, on se lance dans le calcul
    Else
        //methode par la triangularisation (pivot de gauss)
        //le determinant d'une matrice triangulaire est le produit de sa diagonale
        //on parcourt la matrice ligne par ligne
        Pour LignePivot <- 1 à Dimension
            //si le coefficient de la diagonale est nul
            Si (TabCoeffs(LignePivot, LignePivot) = 0) Alors
                //on cherche a permutter la ligne avec une autre, dont le coeff n'est pas nul
                LignePermutation = LignePivot + 1
                //on parcourt les lignes en dessous de la ligne actuelle
                Tant que ((LignePivot + 1 <= Dimension) et (PermutationEffectuee = Faux))
                    //des qu'un coeff n'est pas nul
                    Si (TabCoeffs(LignePermutation, LignePivot) <> 0) Alors
                        //on inverse les deux lignes
                         PermutterLigne(LignePivot, LignePermutation)
                        PermutationEffectuee = Vrai
                     FinSi
                //on s'arrete si on a parcouru toutes les lignes
                //ou si on a effectue une permutation
                 FinTq
             FinSi
 
            //on met les coefficients dans la colonne du pivot a zero
            Pour Ligne <- LignePivot + 1 à Dimension - 1
                 //Lj <- Lj - (aji/aii) Li
                 OperationSurLigne(Ligne, LignePivot, -TabCoeffs(LignePivot, Ligne) / TabCoeffs(Ligne, Ligne))
             FinPour
         FinPour
 
        //on calcule le determinant en faisant le produit de la diagonale
        DeterminantCalcule = 1
        Pour Ligne <- 1 à Dimension
            DeterminantCalcule = DeterminantCalcule * TabCoeffs(Ligne, Ligne)
        FinPour
     FinSi
 
    Retourner(DeterminantCalcule)

il n'y a pas les declarations de variables ni les fonctions annexes, mais c'est deja assez long comme ca ;)

j'ai pas encore pu vérifier si ca marchait bien, je vais m'y pencher a present :)

Citation:
Envoyé par Peltchag
voila ce que j'ai fait a present :
Code:

1 2 //des qu'on tombe sur un coefficient different de zero Si (TabCoeffs(Ligne, Colonne) <> 0) Alors
Avec des flottants, il ne faut jamais comparer directement avec zero, mais tester si la valeur absolue est inferieure a une constante (le langage peut en fournir une adaptee a la precision du systeme).
Citation:
Envoyé par Peltchag
Code:

1 2 3 4 5 6 7 8 9 10 Pour LignePivot <- 1 à Dimension Si (TabCoeffs(LignePivot, LignePivot) = 0) Alors LignePermutation = LignePivot + 1 Tant que ((LignePivot + 1 <= Dimension) et (PermutationEffectuee = Faux)) Si (TabCoeffs(LignePermutation, LignePivot) <> 0) Alors PermutterLigne(LignePivot, LignePermutation) PermutationEffectuee = Vrai FinSi FinTq FinSi
Meme remarque pour le test de nullite d'un flottant.
Apres cette boucle, il faut verifier si PermutationEffectuee est vraie ou non. Si elle est fausse, alors c'est un cas d'arret et le determinant est nul. (cf. ce que j'ai dit plus haut)
Citation:
Envoyé par Peltchag
Code:

1 2 3 //on met les coefficients dans la colonne du pivot a zero Pour Ligne <- LignePivot + 1 à Dimension - 1 //Lj <- Lj - (aji/aii) Li
Le 1/aii devrait etre calcule avant (ok, pour une matrice 4x4, le gain est faible. Mais ici l'algo est ecrit pour etre independant de l'ordre de la matrice, donc ca compte).

Citation:

Envoyé par Peltchag

Code:

OperationSurLigne(Ligne, LignePivot, -TabCoeffs(LignePivot, Ligne) / TabCoeffs(Ligne, Ligne))

"OperationSurLigne" ? Pas tres explicite, comme nom...

29/07/2005, 08h04
Invité

merci pour les remarques :)

mais je ne vais pas travailler avec des flottants moi, juste des decimaux (2 ou 5 chiffres apres la virgule, jamais plus). de plus, mon systeme d'equations concerne une production/consommation, et il y aura forcement des valeurs a zero dans la matrice, c'est pourquoi je suis parti sur de tels tests en fait.

pour la procedure, je sais que le nom est un peu sordide, mais ne me souvenant plus du nom mathematique de l'operation, j'ai mis un nom vite fait (flemme oblige), et ai ajoute un commentaire juste au dessus pour dire ce que cetait ;)

je vais tenir compte de tout ce que tu m'as dit, et y apporte une petite retouche, merci bien !

Retouche effectuee :)

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    //----- parcours de la matrice pour verrifier que son determinant n'est pas nul ----- 
    Indice = 1 
     Répéter 
        //on cherche d'abord sur la ligne 
        LigneNulle = Vrai 
        //on se place sur la ligne voulu 
        Ligne = Indice 
        //et on parcourt chaque colonne 
        Colonne = 1 
         Répéter 
            //des qu'on tombe sur un coefficient different de zero 
            Si (TabCoeffs(Ligne, Colonne) <> 0) Alors 
                //on dit que la ligne n'est pas nulle 
                LigneNulle = Faux 
             FinSi 
            //on passe a la colonne suivante 
            Colonne = Colonne + 1 
        //on sort s'il n'y a plus de colonnes 
        //ou si on sait que la ligne n'est pas nulle 
         Jusqu 'à ((Colonne > Dimension) Ou (LigneNulle = Faux)) 
 
        //si la ligne est bonne, on regarde la colonne 
        //par le meme principe 
        Si (LigneNulle = Faux) Alors 
            ColonneNulle = Vrai 
            Colonne = Indice 
            Ligne = 1 
             Répéter 
                Si (TabCoeffs(Ligne, Colonne) <> 0) Alors 
                    ColonneNulle = Faux 
                 FinSi 
                Ligne = Ligne + 1 
             Jusqu 'à ((Ligne > Dimension) Ou (ColonneNulle = Faux)) 
         FinSi 
 
        //une fois qu'on a regarde ligne et colonne, 
        //on passe a l'indice suivant 
        Indice = Indice + 1 
    //on s'arrete si on a parcouru toute la matrice 
    //ou si on a trouve une ligne ou une colonne nulle 
     Jusqu 'à ((Indice > Dimension) Ou (LigneNulle = Vrai) Ou (ColonneNulle = Vrai)) 
 
    //si une ligne ou une colonne est nulle, le determinant est nul 
    Si (LigneNulle = Vrai) Ou (ColonneNulle = Vrai) Alors 
        DeterminantCalcule = 0 
    //sinon, on se lance dans le calcul 
    Else 
        //----- methode par la triangularisation (pivot de gauss) ----------------------------
        //----- le determinant d'une matrice triangulaire est le produit de sa diagonale -----
        //on parcourt la matrice ligne par ligne 
        LignePivot = 1
        Répéter
            //initialisation du booleen d'arret de la fonction
            ArretFonction = Faux
 
            //si le coefficient de la diagonale est nul
            Si (TabCoeffs(LignePivot, LignePivot) = 0) Alors 
                //on cherche a permutter la ligne avec une autre, dont le coeff n'est pas nul 
                LignePermutation = LignePivot + 1 
                //on parcourt les lignes en dessous de la ligne actuelle 
                Tant que ((LignePivot + 1 <= Dimension) et (PermutationEffectuee = Faux)) 
                    //des qu'un coeff n'est pas nul 
                    Si (TabCoeffs(LignePermutation, LignePivot) <> 0) Alors 
                        //on inverse les deux lignes 
                        PermutterLigne(LignePivot, LignePermutation) 
                        PermutationEffectuee = Vrai 
                    FinSi 
                //on s'arrete si on a parcouru toutes les lignes 
                //ou si on a effectue une permutation 
                FinTq
 
                //si pas de permutation effectuee, on s'arrete car un coefficient de la diagonale est nul
                Si (PermutationEffectuee = Faux) Alors
                    ArretFonction = Vrai
                FinSi
             FinSi
 
             Si (ArretFonction = Faux) Alors            
                 //on met les coefficients dans la colonne du pivot a zero 
                 Pour Ligne <- LignePivot + 1 à Dimension - 1 
                     //calcul du coefficient pour l'opération
                     CoefficientMultiplication = -TabCoeffs(Ligne, Ligne) / TabCoeffs(LignePivot, Ligne)
                     //Ly <- Ly - u Lx 
                     OperationSurLigne(Ligne, LignePivot, CoefficientMultiplication) 
                 FinPour 
            FinSi 
        Jusqu'à ((LignePivot > Dimension) Ou (ArretFonction = Vrai))
 
        //----- Calcul du determinant -----
        //Si on est sorti a cause d'un coeff sur la diagonale nul, le determinant est nul
        Si (ArretFonction = Vrai) Alors
            DeterminantCalcule = Faux
        Sinon 
            //on calcule le determinant en faisant le produit de la diagonale 
            DeterminantCalcule = 1 
            Pour Ligne <- 1 à Dimension 
                DeterminantCalcule = DeterminantCalcule * TabCoeffs(Ligne, Ligne) 
            FinPour 
        FinSi
    FinSi 
 
    Retourner(DeterminantCalcule)

Citation:

Envoyé par Peltchag

mais je ne vais pas travailler avec des flottants moi, juste des decimaux (2 ou 5 chiffres apres la virgule, jamais plus).

En machine, ca s'appelle des flottants; et les flottants sont rarement representables exactement sur un ordinateur...
Quand tu calcules les nouveaux coefficients, par exemple

Code:

a'21 = a21 - (a21/a11)*a11

, les arrondis materiels feront que tu ne retomberas pas exactement sur zero, mais sur 10^-5, 10^-15... et le test (a'21 == 0) echouera. Tu me diras que pour a'21, ce n'est pas grave... Mais si, par exemple, a11 = a12 et a21 = a22, alors a'22 aura la meme valeur "presque nulle". Et la deuxieme iteration du pivot fera des trucs bizarres... Tu peux essayer, pour voir les resultats numeriques qui sortent, ca pourra montrer des ordres de grandeur des erreurs commises (ce qui est toujours interessant a connaitre).
Citation:
Code:

1 2 3 4 5 6 7 Si (ArretFonction = Faux) Alors Pour Ligne <- LignePivot + 1 à Dimension - 1 CoefficientMultiplication = -TabCoeffs(Ligne, Ligne) / TabCoeffs(LignePivot, Ligne) //Ly <- Ly - u Lx OperationSurLigne(Ligne, LignePivot, CoefficientMultiplication) FinPour FinSi
1) pourquoi jusqu'a (Dimension -1) ?
2) Si tu mets la formule en commentaire, mais que tu utilises d'autres notations, c'est moins comprehensible... Ly, c'est Ligne ou LignePivot ? Si on connait l'algo, on comprend, mais quelqu'un relit apres... il aura un peu de mal a suivre (de oute facon, ca doit etre precise dans la doc de la fonction que tu appelles).
3) L'operation, c'est une combinaison lineaire... (enfin, une combinaison particuliere).

Citation:

Envoyé par Peltchag

mais je ne vais pas travailler avec des flottants moi, juste des decimaux (2 ou 5 chiffres apres la virgule, jamais plus).

En machine, ca s'appelle des flottants; et les flottants sont rarement representables exactement sur un ordinateur...
Quand tu calcules les nouveaux coefficients, par exemple

Code:

a'21 = a21 - (a21/a11)*a11

, les arrondis materiels feront que tu ne retomberas pas exactement sur zero, mais sur 10^-5, 10^-15... et le test (a'21 == 0) echouera. Tu me diras que pour a'21, ce n'est pas grave... Mais si, par exemple, a11 = a12 et a21 = a22, alors a'22 aura la meme valeur "presque nulle". Et la deuxieme iteration du pivot fera des trucs bizarres... Tu peux essayer, pour voir les resultats numeriques qui sortent, ca pourra montrer des ordres de grandeur des erreurs commises (ce qui est toujours interessant a connaitre).
Citation:
Code:

1 2 3 4 5 6 7 Si (ArretFonction = Faux) Alors Pour Ligne <- LignePivot + 1 à Dimension - 1 CoefficientMultiplication = -TabCoeffs(Ligne, Ligne) / TabCoeffs(LignePivot, Ligne) //Ly <- Ly - u Lx OperationSurLigne(Ligne, LignePivot, CoefficientMultiplication) FinPour FinSi
1) pourquoi jusqu'a (Dimension -1) ?
2) Si tu mets la formule en commentaire, mais que tu utilises d'autres notations, c'est moins comprehensible... Ly, c'est Ligne ou LignePivot ? Si on connait l'algo, on comprend, mais si quelqu'un relit apres... il aura un peu de mal a suivre (de toute facon, ca doit etre precise dans la doc de la fonction que tu appelles).
3) L'operation, c'est une combinaison lineaire... (enfin, une combinaison particuliere).

29/07/2005, 10h22
Invité

je ne savais pas pour les flottants, j'apprends quelque chose ! donc il faut que je remplace le test par un test :

Code:

Si (abs(coeff) < 10^-x) Alors

mais quelle valeur doit avoir x, si je travaille avec des decimaux qui ont 5 chiffres apres la virgule ?

pour repondre aux remarques :
1. euh... parce que j'ai pas tres bien dormi :oops:
2. c'est vrai, mais bon, celui qui lit ce code, il a interet a connaitre le pivot de gauss de toute facon !!! :wink:
3. merci, j'avais oublie! c'est loin tout ca :P
29/07/2005, 16h54
alveric

Ce n'es pas forcement un nombre de la forme 10^(-x). Suivant le langage et le systeme, tu peux avoir des constantes predefinies, qui donnent par exemple la valeur strictement positive la plus petite representable sur le systeme (avec un nom comme EPSILON), et qui permet de dire "A = B" quand tu as en fait abs(A-B) < EPSILON. Ada permet meme de la configurer, il me semble.
Pour ce qui est de te donner une valeur acceptable, bah surement quelques ordres de grandeur inferieur a 10^(-5). Apres, une valeur precise...
31/07/2005, 14h11
Matthieu Brucher

Heureusement que j'ai dit que la manière la plus simple était par les matrices de Houseolder... Une SVD complète est inutile, il faut juste la partie capable de bidiagonaliser une matrice.
Si ton code a une complexité supérieure à O(n^3), prends cette solution. En plus la GSL propose ça "facilement".
Il y a pleins de librairies OS de calcul matriciel, si vous pouvez les utiliser, utilisez-les ! - je n'ai pas pu le faire dans mon cas et ma librairie se retrouve à être complexe... -
31/07/2005, 22h02
alveric

Citation:

Envoyé par Miles

Si ton code a une complexité supérieure à O(n^3), prends cette solution.

Justement,je ne me souviens plus de la complexité du pivot... Et j'ai un peu la flemme de la recalculer :roll:
01/08/2005, 09h34
Matthieu Brucher

Tout dépend comment on calcule ! La méthode Houseolder fait intervenir une multiplication par un vecteur puis le vecteur transposé. C'est la matrice de Housolder. Si on joue à l'idiot et qu'on calcule d'abord la matrice, on se retrouve avec une multplication matricielle pour chaque colonne du tableau, donc du O(n^4), tandis que si on fait d'abord le vecteur puis le vecteur transposé, ce ne sera plus que du O(n^3).
14/08/2005, 12h17
DrTopos

Bonjour.

Citation:

Envoyé par alveric

(vu que le determinant est aussi invariant par permutation de lignes)

Citation:

Envoyé par alveric

Si tu en trouves un, tu permutes les deux lignes et passes a la suite.

Citation:

Envoyé par Peltchag

PermutterLigne(LignePivot, LignePermutation)

Juste une petite remarque, qui évitera peut-être un bug: s'il est vrai que le déterminant est invariant par ajout à une ligne d'une combinaison linéaire des autres lignes, il n'est pas invariant par permutation de deux lignes. Dans ce cas, il change de signe (sauf si on permute une ligne avec elle-même, bien sûr). La même chose vaut pour les colonnes. En d'autres termes, le déterminant n'est pas seulement multilinéaire. Il est multilinéaire alterné.

Cordialement.
15/08/2005, 10h40
alveric

Bonjour,

Citation:

Envoyé par DrTopos

En d'autres termes, le déterminant n'est pas seulement multilinéaire. Il est multilinéaire alterné.

...et meeeeeeerde........ :boulet:
Voila ce que c'est de dire des choses sans verifier dans son cours de math... :oops:
16/08/2005, 10h11
Matthieu Brucher

J'avais pas vu cette ânerie non plus... :D
16/08/2005, 13h11
Rafy
La methode du pivot de gauss pour triangulariser la matrice est efficace, mais elle peut être employée de deux manières...
Tu peut prendre le plus grand terme pour pivot, alors ta précision sera plus grande mais les calculs plus lourd. Avec le terme le plus petit en pivot, c'est l'inverse... Tu peux aussi essayer de trier les colonnes pour voir s'il n'y a pas de zéro bien placé, et aussi si la matrice est symétrique.
Dans le cas symétrique, c'est super !!!!! Car tu peux mettre ta matrice sous la forme :
Code:

1 2 A = (B)*(^tB)
(Comprendre tansposée)
avec A ta matrice, et B une matrice colonne....
16/08/2005, 13h25
Matthieu Brucher

...
Si A est symétrique, tu ne peux pas la mettre sous cette forme à moins qu'elle soit de rang 1 !
Ton écriture est sans doute celle de Cholesky où B est triangulaire ?

Au fait, pourquoi les calculs sont plus lourds en prenant le pivot l'élément le plus important ?
16/08/2005, 13h30
Rafy

Oui je suis allé trop vite effectivement.. Mais B n'est pas obligé d'être triangulaire...
Et oui pour les matrices de rang 1 c'est nickel !!! Pardon trop rapide.
:sm: Allez-y
Les calculs sont plus lourd car il faut diviser les coefficients, dans l'autre cas il faut les multiplier... C'est pas pareil
16/08/2005, 13h31
alveric
Citation:
Envoyé par Rafy

Dans le cas symétrique, c'est super !!!!! Car tu peux mettre ta matrice sous la forme :

Code:

1 2 A = (B)*(^tB)

(Comprendre tansposée)
avec A ta matrice, et B une matrice colonne....
Autrement dit, tu passes d'une matrice symetrique d'ordre n, donc un element d'un espace de dimension n*(n+1)/2, a un vecteur element d'un espace de dimension n. Il n'y aurait pas une erreur ? Je suis plutot de l'avis de Miles.
Et idem pour la "lourdeur" des calculs avec l'element le plus "important" (en valeur absolue, je suppose ?).

EDIT: mes reponses m'ont devance de peu...
16/08/2005, 13h34
Rafy

Ben ouais en valeur absolue... Sinon ça n'a pas de sens.
Il faut savoir que le calcul d'un déterminant c'est d'ordre !n (par la téchnique standard de développement par rapport à une ligne ou colonne), donc presque toutes les methodes sont mieux ! :twisted:
16/08/2005, 13h40
Rafy

Si les matrices dont il faut déterminer le déterminant sont des matrices de type directX (Transformation spaciale) il y a déjà des methodes pour le faire. (En plus si les matrices sont orthogonales ! Leur déterminant vaut 1 en valeur absolue... iark).
16/08/2005, 14h33
Matthieu Brucher

Citation:

Envoyé par Rafy

Oui je suis allé trop vite effectivement.. Mais B n'est pas obligé d'être triangulaire...
Et oui pour les matrices de rang 1 c'est nickel !!! Pardon trop rapide.
:sm: Allez-y
Les calculs sont plus lourd car il faut diviser les coefficients, dans l'autre cas il faut les multiplier... C'est pas pareil

Ah ? pour moi, il y a toujours une division quand même...

Oui, B n'est pas forcément triangulaire, mais c'est parce que ça rappelle beaucoup la décomposition de cholesky.
16/08/2005, 14h37
alveric

Citation:

Envoyé par Miles

Citation:

Envoyé par Rafy

Oui je suis allé trop vite effectivement.. Mais B n'est pas obligé d'être triangulaire...
Et oui pour les matrices de rang 1 c'est nickel !!! Pardon trop rapide.
:sm: Allez-y
Les calculs sont plus lourd car il faut diviser les coefficients, dans l'autre cas il faut les multiplier... C'est pas pareil

Ah ? pour moi, il y a toujours une division quand même...

Pour moi aussi. Je vois comment reduire le nombre total de divisions a (n-1) pour un determinant d'ordre n, mais de la a s'en passer...
16/08/2005, 16h36
Rafy

Heu dans la page tout le monde me fouette !!!!! Ca fait mal.
Attendez, je réfléchi....
hum....
C'est vrai, que le nombre de division n'est pas réduit à zéro.
Le truc ça vient de la précision, voila c'est ca.
Si on fait la division, puis la multiplication, on a pas la même précision que si on fait l'inverse... Voila. Ce coup ci c'est bon... Pas de fouet.
De toute manière, ça coute cher de calculer un déterminant...
Es ce que l'appel est fréquent ou pas ?
C'est pour un jeu, dans chaque frame... Es ce que des lignes ou colonnes de la matrice restent les mêmes entre deux calculs ?

Voir 40 message(s) de cette discussion en une page