Mathc matrices/b0a4c

Un livre de Wikilivres.
Aller à la navigation Aller à la recherche


Sommaire


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


Crystal Clear mimetype source c.png b0a4c.c
'
/* ------------------------------------ */
/*  Save as :  b0a4c.c               */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
void Xgj_freevariable_eignvector_mR(
double **Ab,
double **eignV
)
{
double **Ab_T = c_mR(Ab,i_duplicateAbR0_Ab_c_A_c_mR(
                             Ab[R_SIZE][C0],
                             Ab[C_SIZE][C0],
                             Ab[C_SIZE_A][C0]));
double **new_Ab;
int r_not_0;
int r;

  gj4_T_mR(Ab_T);
  sort4_r_coef_r0_mR(Ab_T);
  r_not_0 =  is_system_consistent_mR(Ab,Ab_T);
  
  
      new_Ab = i_duplicateAbR0_Ab_c_A_c_mR(
                             Ab[C_SIZE_A][C0],
            Ab[C_SIZE][C0]+((Ab[C_SIZE_A][C0]-C1)-(r_not_0-C1)),
                             Ab[C_SIZE_A][C0]);

      c_n_r_mR(Ab,r_not_0,new_Ab);
      zero4_under_all_pivot_mR(new_Ab);
      put_freevariable_mR(new_Ab,r_not_0);

      r=Ab[C_SIZE_A][C0]-C1;
      while( r>R1 )
        zero_below_pivot_gj1Ab_mR(new_Ab,r--);

      sort_c_mR(new_Ab);
      sort_r_mR(new_Ab);

      c_Ab_b_mR(new_Ab,eignV);

             
      f_mR(Ab_T);
      f_mR(new_Ab);
}
/* ------------------------------------ */
/* ------------------------------------ */
void XeignVectoruv_mR(
double **A,
double **EignVector
)
{
int i;	
int r = rsize_mR(A);
	
double **EigsValue = i_mR(r,C1);
double **eignV     = i_mR(r,r);

double **Ab;
double **b         =   m0_mR(i_mR(r,C1));

double **Ide       =  eye_mR(i_mR(r,r));
double **sIde      =         i_mR(r,r);
double **sIde_A    =         i_mR(r,r);

	
  eigsuv_mR(A,EigsValue);

/* ********** DEBUG 1 ************* */  
  printf(" EigsValue:");
  p_mR(EigsValue,S10,P4,C6);
/* ******************************** */ 
  
  
  for(i = R1; i <= rsize_mR(EigsValue); i++)
  {  
   smul_mR(c_coef_R(EigsValue,i,C1),Ide,sIde);
   sub_mR(sIde,A,sIde_A);   
   Ab = c_A_b_Ab_mR(sIde_A,b,i_AbR0_mR(r,r,C1));
   Xgj_freevariable_eignvector_mR(Ab,eignV);
   c_c_mR(eignV,C2,EignVector,i);  
  }    
   
  Normalizeuv_mR(EignVector);  

  f_mR(EigsValue);
  f_mR(eignV);
  
  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(sIde_A);
}
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{
double ab[R4*C4]={
+150,-925,+32,-475,
-925,-253,+630,+302,
+32,+630,-183,-59,
-475,+302,-59,+668
};

int r = RC4;

double **A         = ca_A_mR(ab,i_mR(r,r));
double **EignVector=               i_mR(r,r);
double **T_EignVector=             i_mR(r,r);
double **EigsValue =               i_mR(r,C1);

double **T1=             i_mR(r,r);
double **T2=             i_mR(r,r);

  clrscrn();   
  XeignVectoruv_mR(A,EignVector); 
  printf(" EignVector:");
  p_mR(EignVector,S10,P4,C6);
  stop();
   
  f_mR(A);
  
  f_mR(EignVector);
  f_mR(T_EignVector);
  
  f_mR(EigsValue);  
     
  f_mR(T1);
  f_mR(T2);
}
/* ------------------------------------ */
int main(void)
{

  fun();

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


Vérifions dans un premier temps que les valeurs propres calculées par la fonction XeignVectoruv_mR(); correspondent bien à celles données par octave.


Exemple de sortie écran :
 ------------------------------------ 
 EigsValue:
+1394.5038 
-1256.5800 
 +405.3394 
 -161.2632 

 EignVector:
   +0.6008    +0.6008    +0.1824    +0.6000 
   -0.5106    -0.5106    -0.4330    +0.0296 
   -0.1696    -0.1696    -0.5249    +0.7018 
   -0.5913    -0.5913    +0.7098    +0.3829 

 Press return to continue.
Exemple de sortie Octave :
  >>  a=[
+150,-925,+32,-475;
-925,-253,+630,+302;
+32,+630,-183,-59;
-475,+302,-59,+668]
a =

   150  -925    32  -475
  -925  -253   630   302
    32   630  -183   -59
  -475   302   -59   668

>>
>>  [V, E] = eigs (a,10)
V =

  -0.6007990  -0.4957134   0.1823898   0.6000357
   0.5105608  -0.7422830  -0.4329758   0.0295902
   0.1695987   0.4508004  -0.5248908   0.7017870
   0.5912737   0.0079513   0.7097574   0.3828533

E =

Diagonal Matrix

   1394.50         0         0         0
         0  -1256.58         0         0
         0         0    405.34         0
         0         0         0   -161.26

>>