Calcul de valeurs propres

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...

Quelqu'un de vous est-il familier avec ceci ?

Merci

Reply

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

Reply

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

http://lumimath.univ-mrs.fr/~jlm/t [...] ode21.html

Message cité 1 fois

Reply

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

Joel F a écrit :

http://lumimath.univ-mrs.fr/~jlm/t [...] ode21.html

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

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

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

Message cité 1 fois

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

Message cité 1 fois

---------------
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
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

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 :

complexe *get_eig2D(...bla bla...){
double R, S;
complexe A[2][2], delta, Z1, Z2, *D,
zero={0, 0}, deux={2, 0}, quatre={4, 0};
//blabla... remplissage de la matrice A
// delta = (a+d)²-4(ad-bc)
delta = c_sub(c_mul(c_add(A[0][0], A[1][1]), c_add(A[0][0], A[1][1])),
c_mul(quatre, c_sub(c_mul(A[0][0], A[1][1]), c_mul(A[0][1], A[1][0])))
);
// R = +-sqrt( (re(delta)+|delta|)/2 )
// S = +-sqrt( -(re(delta)-|delta|)/2 )
R = sqrt((delta.re + sqrt(delta.re*delta.re + delta.im*delta.im))/2);
S = sqrt(fabs((delta.re - sqrt(delta.re*delta.re + delta.im*delta.im))/2));
if(delta.im > 0){
Z1.re = R; Z1.im = S;
Z2.re = -R; Z2.im = -S;
}else{
Z1.re = R; Z1.im = -S;
Z2.re = -R; Z2.im = S;
}
//Eigenvalues are ((a+d)+Z1)/2 and ((a+d)+Z2)/2
D=(complexe *)malloc(sizeof(complexe)*2);
D[0] = c_div(deux ,c_add(c_add(A[0][0], A[1][1]), Z1));
D[1] = c_div(deux ,c_add(c_add(A[0][0], A[1][1]), Z2));
return D;
}

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 ç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

Calcul de valeurs propres

Sujets relatifs:

Leave a Replay