Diagonalisation de matrices 3x3 symetriques

Soyez le premier à donner votre avis sur cette source.

Snippet vu 18 634 fois - Téléchargée 31 fois

Contenu du snippet

voici un algorithme qui permet de diagonaliser une matrice symetrique
pour l'exemple c'est une matrice 3x3
comme la matrice est symetrique, elle est diagonalisable et de plus
il est facile de demontrer que les vecteurs propres sont orthogonaux

Source / Exemple :


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

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

#include <conio.h>

//-------------------------------------------------
typedef struct tagEIGEN_SPACE
  {
  double  lambda; // valeur propre
  VECT3D  X;      // vecteur propre
  }EIGEN_SPACE,*P_EIGEN_SPACE,**PP_EIGEN_SPACE;

//-------------------------------------------------
void print(P_EIGEN_SPACE E)
{
printf("lambda = %lf   et   X = {%lf ; %lf ; %lf}\n",E->lambda,E->X.x,E->X.y,E->X.z);
} // print()

//-------------------------------------------------
double eigenValue(P_VECT3D Y,P_VECT3D X)
{
double normY,normX;
normY = NormVect3D(Y);
normX = NormVect3D(X);
if(CmpDouble0(normX))
  return 0.;
else
  return normY/normX;
} // eigenValue()

//-------------------------------------------------
BOOL diag(  P_MAT3X3      M,
            P_EIGEN_SPACE E1,
            P_EIGEN_SPACE E2,
            P_EIGEN_SPACE E3)
{
int     i;
VECT3D  X1,X2,X3;
VECT3D  Y1,Y2,Y3;

// initialement on est sur la base canonique
SetVect3D(&X1,1.,0.,0.);
SetVect3D(&X2,0.,1.,0.);
SetVect3D(&X3,0.,0.,1.);

for(i=0;i<10;i++)
  {
  // on calcule les images
  LinAppVect3D(&X1,M,&X1);
  LinAppVect3D(&X2,M,&X2);
  LinAppVect3D(&X3,M,&X3);

  // on normalise <X1>
  UnitVect3D(&X1,&X1);

  // on orthogonalise <X2> par rapport a <X1>
  // puis on normalise <X2>
  OrthoVect3D(&X1,&X2);
  UnitVect3D(&X2,&X2);

  // on orthogonalise <X3> par rapport a <X1> puis <X2>
  // puis on normalise <X3>
  OrthoVect3D(&X1,&X3);
  OrthoVect3D(&X2,&X3);
  UnitVect3D(&X3,&X3);
  }
E1->X = X1;
E2->X = X2;
E3->X = X3;

// on calcule maintenant les valeurs propres
LinAppVect3D(&Y1,M,&X1);
LinAppVect3D(&Y2,M,&X2);
LinAppVect3D(&Y3,M,&X3);
E1->lambda  = eigenValue(&Y1,&X1);
E2->lambda  = eigenValue(&Y2,&X2);
E3->lambda  = eigenValue(&Y3,&X3);

return TRUE;
} // diag()

//-------------------------------------------------
int main(int argc,char **argv)
{
InitLibUtil();
  {
  double        a,b,c,d,e,f;
  MAT3X3        mat;
  EIGEN_SPACE   E1,E2,E3;

  // (a b c)
  // (b d e)
  // (c e f)
  a = 3;
  b = 1;
  c = 2;
  d = 3;
  e = 2;
  f = 3;
  mat.a11 = a;mat.a12 = b;mat.a13 = c;
  mat.a21 = b;mat.a22 = d;mat.a23 = e;
  mat.a31 = c;mat.a32 = e;mat.a33 = f;

  diag(&mat,&E1,&E2,&E3);

  print(&E1);
  print(&E2);
  print(&E3);

  getch();
  }
CloseLibUtil();
return 0;
} // main()

Conclusion :


le resultat est :

lambda = 6.372281 et X = {0.541783 ; 0.541766 ; 0.642621}
lambda = 2.000000 et X = {-0.707100 ; 0.707113 ; 0.000008}
lambda = 0.627719 et X = {-0.454401 ; -0.454401 ; 0.766185}

A voir également

Ajouter un commentaire

Commentaires

frisou65
Messages postés
1
Date d'inscription
mardi 23 août 2011
Statut
Membre
Dernière intervention
25 août 2011

Bonjour,

est ce qu'il existe le même code mais en VBA car je voudrais diagonaliser une matrice 3*3 via un fichier Excel ?

Cordialement.
sebastien_lacaze01
Messages postés
1
Date d'inscription
mercredi 4 juillet 2007
Statut
Membre
Dernière intervention
3 septembre 2007

Salut,

je viens de creer un projet avec le main ( de matrice 3x3) et les sources du projet diagNxN (match et util). Est ce que c'est tjs compatible car j'ai un pbm dans l'execution?

d'avance merci de ta réponse
cs_JCDjcd
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
2
Stiko
Messages postés
38
Date d'inscription
jeudi 29 juin 2006
Statut
Membre
Dernière intervention
20 février 2008

est ce que tu peu mettre le lien stp??
cs_JCDjcd
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
2
j'ai fait une nouvelle source pour le cas général NxN
Afficher les 8 commentaires

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.