Calcul de valeurs propres

Calcul de valeurs propres - Algo - Programmation

Marsh Posté le 30-06-2006 à 23:41:20    

Hello à tous :)
 
Alors voilà, je dois écrire (en C, mais bon pour l'instant on s'en fiche du langage d'implémentation) un algo pour récupérer les valeurs propres d'une matrice 2x2 à coefficients complexes
J'ai un peu cherché sur le net, mais impossible de mettre la main sur un pseudo-algo sympa, je pense que je dois regarder du côté de l'algo QR mais je sais pas trop...
 
Dire qu'en matlab la fonction eig(X) est native... :d
 
Quelqu'un de vous est-il familier avec ceci ?
 
Merci :d

Reply

Marsh Posté le 30-06-2006 à 23:41:20   

Reply

Marsh Posté le 01-07-2006 à 10:11:03    

Reply

Marsh Posté le 01-07-2006 à 11:50:18    


 

Citation :

Soit A une matrice n*n symétrique réelle.


 
 
regarde du côté de la librairie LAPACK


---------------
What if I were smiling and running into your arms? Would you see then what I see now?  
Reply

Marsh Posté le 01-07-2006 à 18:43:52    

zut j'ai mal lu ...
sinon ouias LAPACK ou son interface C++ :
http://lpp.sourceforge.net


Message édité par Joel F le 01-07-2006 à 18:44:08
Reply

Marsh Posté le 01-07-2006 à 18:56:05    

mmh oui je viens de regarder LAPACK, ça serait exactement ça qu'il me faudrait :
 
http://www.nacse.org/demos/lapack/codes/eigen-c.html
http://www.nacse.org/demos/lapack/routines/zgeev.html
 
le problème c'est que le code C appelle la routine zgeev en fortran...? c'est un peu embêtant ça quand même :d
 
mais merci pour cette piste :jap:

Message cité 1 fois
Message édité par Zipo le 01-07-2006 à 18:56:19
Reply

Marsh Posté le 01-07-2006 à 19:31:18    

bon en fait je crois que je vais recoder la fonction à la main...
 
Comme c'est tout le temps des 2x2 c'est pas très compliqué

Reply

Marsh Posté le 02-07-2006 à 00:24:51    

Zipo a écrit :

mmh oui je viens de regarder LAPACK, ça serait exactement ça qu'il me faudrait :
 
http://www.nacse.org/demos/lapack/codes/eigen-c.html
http://www.nacse.org/demos/lapack/routines/zgeev.html
 
le problème c'est que le code C appelle la routine zgeev en fortran...? c'est un peu embêtant ça quand même :d
 
mais merci pour cette piste :jap:


 
qu'est-ce qu'il y a d'embétant, tu peux le faire... [:jagstang]


---------------
What if I were smiling and running into your arms? Would you see then what I see now?  
Reply

Marsh Posté le 02-07-2006 à 01:05:38    

pas besoin de lapack si tu as une matrice 2x2 car tu as une formule analytique :
par definition : r est vp si det(A-r*I)=0
 
si A = (a11 a12; a21 a22) alors : r² - (a11+a22)*r + a11*a22 - a12*a21 = 0
il suffit de resoudre cette equation du second degre et la encore tu as des formule analytique

Reply

Marsh Posté le 02-07-2006 à 01:13:43    

joneal a écrit :

pas besoin de lapack si tu as une matrice 2x2 car tu as une formule analytique :
par definition : r est vp si det(A-r*I)=0
 
si A = (a11 a12; a21 a22) alors : r² - (a11+a22)*r + a11*a22 - a12*a21 = 0
il suffit de resoudre cette equation du second degre et la encore tu as des formule analytique


lis bien le premier post, surtout la partie en gras


---------------
What if I were smiling and running into your arms? Would you see then what I see now?  
Reply

Marsh Posté le 03-07-2006 à 19:29:28    

jagstang a écrit :

lis bien le premier post, surtout la partie en gras


il n'y a pas de restriction les coefficient rij peuvent etre reel ou complexes, et tu as toujours une equation du second degrés à resoudre, le fait quelle soit complexe n'est pas nécessairemnt un souci
 
@pluche

Reply

Marsh Posté le 03-07-2006 à 19:29:28   

Reply

Marsh Posté le 05-07-2006 à 02:43:23    

c'est bon ça marche impec :d
yavait pas besoin de librairies finalement, c'est bidon à faire  
 
si ça intéresse du monde je met le code, mais sinon bof, vu qu'il y a pas mal de conventions à moi (notamment pour les fonctions de manipulation de complexes) alors ça fait un peu cafouilli :whistle:

Reply

Marsh Posté le 05-07-2006 à 08:52:55    

et pour des matrices d'ordre supérieures ça fontionne ?
 
donne toujours le code


---------------
What if I were smiling and running into your arms? Would you see then what I see now?  
Reply

Marsh Posté le 05-07-2006 à 09:36:40    

ah nan ça marche juste pour mes matrices 2x2 :d
 
Voilà en gros le bout de code, je vous copie pas les fonctions de manipulation de complexes mais en gros dans un autre module ya :
c_add(complexe, complexe);
c_sub(complexe, complexe);
c_mul(complexe, complexe);
c_div(complexe, complexe);
et le type complexe est composé d'un champ double re et d'un champ double im :)
 

Code :
  1. complexe *get_eig2D(...bla bla...){
  2. double R, S;
  3. complexe A[2][2], delta, Z1, Z2, *D,
  4.          zero={0, 0}, deux={2, 0}, quatre={4, 0};
  5. //blabla... remplissage de la matrice A
  6. // delta = (a+d)²-4(ad-bc)
  7. delta = c_sub(c_mul(c_add(A[0][0], A[1][1]), c_add(A[0][0], A[1][1])),
  8.         c_mul(quatre, c_sub(c_mul(A[0][0], A[1][1]), c_mul(A[0][1], A[1][0])))                           
  9.         );
  10. // R = +-sqrt( (re(delta)+|delta|)/2 )
  11. // S = +-sqrt( -(re(delta)-|delta|)/2 )
  12. R = sqrt((delta.re + sqrt(delta.re*delta.re + delta.im*delta.im))/2);
  13. S = sqrt(fabs((delta.re - sqrt(delta.re*delta.re + delta.im*delta.im))/2));
  14. if(delta.im > 0){
  15.             Z1.re = R; Z1.im = S;
  16.             Z2.re = -R; Z2.im = -S;           
  17. }else{
  18.             Z1.re = R; Z1.im = -S;
  19.             Z2.re = -R; Z2.im = S;     
  20. }
  21. //Eigenvalues are ((a+d)+Z1)/2  and  ((a+d)+Z2)/2
  22. D=(complexe *)malloc(sizeof(complexe)*2);
  23. D[0] = c_div(deux ,c_add(c_add(A[0][0], A[1][1]), Z1));
  24. D[1] = c_div(deux ,c_add(c_add(A[0][0], A[1][1]), Z2));
  25. return D;
  26. }

Reply

Marsh Posté le 05-07-2006 à 09:48:09    

just deux remarques :
 
1/ au lieu de gérer ta matrice comme etant [N][N], gére la en [N*N] et fait toi une chtite fonction d'accés.
 
2/ on cast pas devant un malloc [:pingouino]

Reply

Marsh Posté le 19-07-2006 à 06:52:40    

1/ quel interêt ?
 
2/ ouai je sais :d ça fait parti des LU dont il faudrait que je perde l'habitude :/

Reply

Marsh Posté le 19-07-2006 à 18:53:12    

c'est juste moins chiant

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed