Determinants (gauss & cofacteurs)

Télécharger le projet
Commenter

Description

Voici un exemple d'algorithmes pour calculer le determinant d'une matrice.
Le premier algorithme utilise le pivot de Gauss,
et la complexite est en O(n^3)
Le second utilise les cofacteurs (developpement suivant une colonne),
et sa complexite est en O(n!)

Les fonctions se nomment FastDet et StupidDet, je vous laisse deviner lequel correspond avec laquelle.
Je voudrais ici pousser un coup de geule, arreter de calculer un determinant par les cofacteurs
car cette methode est INUTILISABLE, il est impossible de l'utiliser pour une matrice 10x10, c'est
beaucoup trop LENT. Alors que l'autre methode est ultra-rapide.

Source / Exemple : 


//-------------------------------------------------
// calcul d'un determinant
//-------------------------------------------------

#ifndef _UTIL_H_
	#include "util.h"
#endif // _UTIL_H_

#ifndef _MATH_H_
	#include "math.h"
#endif // _MATH_H_

#include <conio.h>

//-------------------------------------------------
typedef double FLT;
//-------------------------------------------------
FLT **CreateMat(int n)
{
FLT **mat; // la matrice
int   y;   // <y> pour les lignes

mat = Malloc(FLT*,n);
for(y=0;y<n;y++)
  {
  mat[y] = Malloc(FLT,n);
  }

return mat;
} // CreateMat()

//-------------------------------------------------
FLT **RandomMat(int n,FLT coeffMin,FLT coeffMax)
{
FLT **mat;
int   x,y;  // <x> pour les colonnes
            // <y> pour les lignes

mat = CreateMat(n);

for(x=0;x<n;x++)
  {
  for(y=0;y<n;y++)
    {
    mat[y][x] = (FLT)RandomDouble(coeffMin,coeffMax);
    }
  }

return mat;
} // RandomMat()

//-------------------------------------------------
void DeleteMat(FLT **mat,int n)
{
int y; // <y> pour les lignes
for(y=0;y<n;y++)
  {
  Free(mat[y]);
  }
Free(mat);
} // DeleteMat()

//-------------------------------------------------
// calcul d'un determinant avec le pivot de Gauss
// la complexite est en n^3
FLT FastDet(FLT **_mat,int n)
{
FLT **mat;  // la matrice "tampon", que l'on pourra modifier
FLT   det;  // la valeur du determinant
int   x,y;  // <x> pour les colonnes
            // <y> pour les lignes
int   i,j;  // <i> pour les colonnes
            // <j> pour les lignes

// on copie <_mat> dans <mat> car on n'a pas le droit
// de modifier les valeurs de la matrice <_mat>
mat = CreateMat(n);
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x] = _mat[y][x];
    }
  }

// on calcule le determinant par la methode des pivots de Gauss
// on rend la matrice triangulaire superieur tout en conservant
// le determinant de la matrice a chaque iteration
det = 1.;

// on balaye les lignes
for(j=0;j<n-1;j++)
  {
  // Le principe est simple, montrons(le sur une exemple 6x6
  // La matrice est sous cette forme :
  // 
  // ( a1 *  *  *  *  * )  [L1]
  // ( 0  a2 *  *  *  * )  [L2]  <-- noms des lignes
  // ( 0  0  a3 *  *  * )  [L3]
  // ( 0  0  0  #1 *  * )  [L4]
  // ( 0  0  0  #2 *  * )  [L5]
  // ( 0  0  0  #3 *  * )  [L6]
  // 
  // Le determinant courant <det> vaut a1.a2.a3
  // Si la matrice est non-singuliere, (i.e. que sont determinant est non-nulle,
  // autrement dit elle est inversible) alors parmis les trois nombres #1 #2 et #3
  // il y en a un qui est non-nul, on prend le plus grand de ceux la en module pour
  // une raison de stabilite numerique.
  // Si #1=#2=#3=0 alors <det> est nul car la matrice est singuliere !
  // Par exemple si le plus grand est #2 alors on echange [L4] et [L5], ce qui va modifier
  // le determinant d'un signe.
  // 
  // La matrice sera donc sous cette forme :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4] (avec a4=#2)
  // ( 0  0  0  #1 *  *  )  [L5]
  // ( 0  0  0  #3 *  *  )  [L6]
  // 
  // Ensuite il suffit de retrancher #1/a4.[L4] a [L5] et #3/a4.[L4] a [L6] pour obtenir :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4]
  // ( 0  0  0  0  *  *  )  [L5]
  // ( 0  0  0  0  *  *  )  [L6]
  // 
  // Le determinant est inchange par cette derniere operation.
  // Le determinant courant vaut maintenant a1.a2.a3.a4
  // A la fin de l'operation on aura <mat> :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4]
  // ( 0  0  0  0  a5 *  )  [L5]
  // ( 0  0  0  0  0  a6 )  [L6]
  // 
  // Remarquons qu'il suffit de faire l'operation <n-1> fois
  // car a la fin on aura a6 directement
  // 
  // En resume voici les differentes etapes :
  //   ( etape 1 ) : on recherche du plus grand coefficient (en module)
  //                 si il est nul, alors le determinant est 0 et l'algorithme est fini !
  //   ( etape 2 ) : on echange des lignes avec eventuellement changement du signe de <det>
  //                 et le determinant courant est multiplie par ledit coefficient
  //   ( etape 3 ) : on retranche la lignes a toutes les lignes du dessous
  //                 pour annuler les coefficients
  int  rankMax,rank;
  FLT  coeffMax;

  // ( etape 1 )
  rankMax = j;
  for(rank=j+1;rank<n;rank++)
    {
    if(fabs(mat[rankMax][j]) < fabs(mat[rank][j]))
      {
      rankMax = rank;
      }
    }

  coeffMax = mat[rankMax][j];
  if(fabs(coeffMax) <= MATH_EPSILON)
    {
    det = 0.;
    goto label_end;
    }
  // ( etape 2 )
  if(rankMax != j)
    {
    for(i=j;i<n;i++)
      {
      FLT tmp;
      tmp             = mat[j][i];
      mat[j][i]       = mat[rankMax][i];
      mat[rankMax][i] = tmp;
      }
    det *= -1.;
    }

  det *= coeffMax;
  // ( etape 3 )
  for(rank=j+1;rank<n;rank++)
    {
    FLT coeff; // c'est dans l'exemple le rapport #3/a4
    coeff = mat[rank][j]/coeffMax;
    for(i=j;i<n;i++)
      {
      mat[rank][i] -= coeff*mat[j][i];
      }
    }
    
  }

det *= mat[n-1][n-1];

// on libere <mat>
label_end:
DeleteMat(mat,n);
return det;
} // FastDet()

//-------------------------------------------------
// calcul d'un determinant avec les cofacteurs
// la complexite est en n!
FLT StupidDet(FLT **mat,int n)
{
// cas trivial d'une matrice 1x1
if(1 == n)
  {
  return mat[0][0];
  }
// cas trivial d'une matrice 2x2
// 
// | a b |
// | c d |  = a.d - c.d
// 
else if(2 == n)
  {
  return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
  }
else
  {
  FLT    **cof; // matrice qui represente un cofacteur
  FLT      det; // le determinant
  int      x,y; // <x> pour les colonnes
                // <y> pour les lignes
  int      i,j; // <i> pour les colonnes
                // <j> pour les lignes
  FLT      sgn;

  // on cree une matrice temporaire pour le calcul des cofacteurs
  cof = CreateMat(n-1);
  for(y=0;y<n-1;y++)
    {
    for(x=0;x<n-1;x++)
      {
      cof[y][x] = mat[y][x];
      }
    }

  // quelques initialisations
  det =  0.0;
  sgn = +1.0;
  // on balaye la derniere colonne (de bas en haut)
  for(j=n-1;j>=0;j--)
    {
    // Le principe est tres simple, il s'agit de developper le determinant
    // suivant la derniere colonne. Cet algorithme est evidemment recursif
    // puisqu'il faut les sous-determinants
    // 
    // Exemple sur une matrice 4x4 :
    // 
    // | *  *  *  a1 |
    // | *  *  *  a2 |
    // | *  *  *  a3 |
    // | *  *  *  a4 |
    // 
    // A chaque coefficient de la derniere colonne, on associe une sous-dertermiant ("sd") :
    // (notation : '#')
    // 
    // ds1 : | *  *  *  a1 |      ds2 : | #  #  #  a1 |
    //       | #  #  #  a2 |            | *  *  *  a2 |
    //       | #  #  #  a3 |            | #  #  #  a3 |
    //       | #  #  #  a4 |            | #  #  #  a4 |
    // 
    // ds3 : | #  #  #  a1 |      ds4 : | #  #  #  a1 |
    //       | #  #  #  a2 |            | #  #  #  a2 |
    //       | *  *  *  a3 |            | #  #  #  a3 |
    //       | #  #  #  a4 |            | *  *  *  a4 |
    // 
    // La formule du calcul du determinant par les cofacteurs est donnee par :
    // 
    //     det = a4.ds4 - a3.ds3 + a2.ds2 - a1.ds1
    //
    //     (notons l'alternance des signes)
    // 
    // On doit calculer les "ds", donc creer les matrices correspondants, et calculer
    // par recurrence leurs determinants. On remarque facilement que les matrices se
    // ressemblent à (toujours) 2 lignes pres. Donc pour eviter de toujours les recalculer,
    // on les initilise de facon "incrementale".
    // 
    
    // on modifie la bonne ligne pour la matrice des cofacteurs
    if(j < n-1)
      {
      for(i=0;i<n-1;i++)
        {
        cof[j][i] = mat[j+1][i];
        }
      }
    // la recurrence en personne !
    det += sgn * mat[j][n-1] * StupidDet(cof,n-1);
    // changement de signe (l'alternance remarquee)
    sgn *= -1.0;
    }

  DeleteMat(cof,n-1);
  return det;
  }
} // StupidDet()
Télécharger le projet

Codes Sources

A voir également

S'inscrire maintenant

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.

Du même auteur (cs_JCDjcd)