Diagonalisation de matrices 3x3 symetriques

Commenter

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

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)