Mathc complexes/02c

Sommaire ◀ Utilise la commande "Retour en Arrière" de ton navigateur.

Installer et compiler ces fichiers dans votre répertoire de travail.

	c00a.c

/* ------------------------------------ */
/*  Save as :   c00a.c                  */
/* ------------------------------------ */
#include "w_a.h"
/* ------------------------------------ */
void fun(int r)
{
double **A          = r_mZ(i_mZ(r,r),99.);
double **EigsValue  =      i_mZ(r,C1);

  clrscrn();
  printf(" Copy/Past into the octave windows \n\n");
  p_Octave_mZ(A,"a",P0,P0);  
  printf(" EigenValues  = eigs (a,%d)\n\n\n",r);
         
  printf(" EigsValue :");
  eigs_mZ(A,EigsValue);  
  p_mZ(EigsValue, S10,P4, S10,P4, C5); 
       
  f_mZ(A);
  f_mZ(EigsValue);
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

do
{
  fun(rp_I(R2)+R3);
} while(stop_w());

  return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */

Je veux calculer les valeurs propres de la matrice M.

J'utilise la QR décomposition de M

  - Puis je calcule RQ 
  - RQ devient la nouvelle valeur de M 
  - Et je recommence un cycle.

Je copie M dans T pour ne pas modifier M. Puis j'applique 1000 fois l'algorithme vue plus haut.

 c_mZ(M,T);
 
 for(i=0;i<1000;i++)
     { 
	    QR_mZ(T,Q,R);  
       mul_mZ(R,Q,T);
     }

À la fin je trouve sur la diagonale de T les valeurs propres.

     
  for(re=R1,ce=C1;re<=r;re++,ce+=C2)
     { 
	   EigsValue[re][C1] = T[re][ce];
	   EigsValue[re][C2] = T[re][ce+C1];
     }

Exemple de sortie écran :

 Copy/Past into the octave windows 

 a=[
-53+82*i,-11-51*i,+88-49*i,+1-16*i,+82+37*i;
-60-6*i,-69+36*i,+51+27*i,-96+4*i,-67-61*i;
-79-78*i,-8-13*i,+93-43*i,-87+13*i,+87+27*i;
+29+38*i,+61-28*i,+90+80*i,+94-47*i,-48-87*i;
-3-92*i,-73-44*i,+49+96*i,+21-98*i,-72-5*i]

 EigenValues  = eigs (a,5)


 EigsValue :
  +74.6084 -164.7657i 
  +59.6316 +153.6893i 
  -78.3801 +129.8447i 
 -131.1392  -66.3198i 
  +68.2793  -29.4485i 


 Press   return to continue
 Press X return to stop

Le code de la fonction eigs_mZ(); :

/* ------------------------------------ */
/* ------------------------------------ */
double **eigs_mZ(
double **A,
double **EigsValue
)
{
int r = rsize_Z(A);

double **T      =      i_mZ(r,r);
double **Q      =      i_mZ(r,r);
double **R      =      i_mZ(r,r);
int i;
int re, ce;

 c_mZ(A,T);
 
 for(i=0;i<1000;i++)
     { 
	   QR_mZ(T,Q,R);  
        mul_mZ(R,Q,T);
     }
     
  for(re=R1,ce=C1;re<=r;re++,ce+=C2)
     { 
	   EigsValue[re][C1] = T[re][ce];
	   EigsValue[re][C2] = T[re][ce+C1];
     }
     
  f_mZ(T);
  f_mZ(Q);
  f_mZ(R); 
  
  return(EigsValue);    
}
/* ------------------------------------ */
/* ------------------------------------ */