RÉSOLUTION DES SYSTÈME AX=B

/* ========================================================================== */
/*   TP informatique SF2A G7 Khelifi Walid 2006/2007                                                 */
/*   Ecole nationale polytechnique, El Harrach Alger                          */
/*   Description : Resolution des systèmes                                    */
/* ========================================================================== */

/*
Résolution de A x = b par :
           
    Cramer , Inversion , Gauss , Gauss pivot partiel , 
    Gauss pivot total , Jordan , Décomposition LU croot , 
    LU doolittle , Cholseky , 
    et les methodes itératives de Jacobie et Gauss-Seidel.
*/

#include <stdio.h>
#include <math.h>

void cramer(float [][],float [],int);
void inversion(float [][],float [],int);
float det(float [][],int);
void comatrices(float [][],float [][],int,int,int);
void gauss(float [][],float [],int);

void gauss_pivot_partiel(float [][],float [],int);
void gauss_pivot_total(float [][],float [],int);
void jordan(float [][],float [],int);
void LU_doolittle(float [][],float [],int);

void LU_crout(float [][],float [],int);
void cholesky(float [][],float [],int);
void jacobie(float [][],float [],int);

void gauss_seidel(float [][],float [],int);
float norme(float [],int);

void remplissage(float [][],float [],int);
void aff_syst(float [][],float [],int);
void aff_mat(float [][],int);
void zero(float [][],float [],int);

main()
{
float a[19][19],b[19];int n,i,j; char valider;

// choix de la methode 
printf("\n *******************************************************\n");    
printf(" *                                                     *\n");      
printf(" *           Resolution de A x = B                     *\n");    
printf(" *                                                     *\n");      
printf(" *******************************************************\n");    
printf(" *                                                     *\n");      
printf(" *   Choix de la methode  :                            *\n"); 
printf(" *                                                     *\n");      
printf(" *      A => Cramer                                    *\n");
printf(" *      B => Inversion                                 *\n");
printf(" *      C => Gauss                                     *\n");
printf(" *      D => Gauss pivot partiel                       *\n");
printf(" *      E => Gauss pivot total                         *\n");
printf(" *      F => Jordan                                    *\n");
printf(" *      G => LU Doolittle                              *\n");
printf(" *      H => LU Crout                                  *\n");
printf(" *      I => Cholesky                                  *\n");
printf(" *      J => Jacobie                                   *\n");
printf(" *      K => Gauss-Seidel                              *\n");
printf(" *                                                     *\n");    
printf(" *******************************************************\n");                  
printf(" *                                                     *\n");    
printf(" *    X => E X I T                                     *\n");
printf(" *                                                     *\n");    
printf(" *******************************************************\n\n");                  

reprendre_choix:

printf("------ => donner votre choix ------ : ");
scanf("%s",&valider);

void remplir()
{
printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
scanf("%d",&n);  
remplissage(a,b,n);
printf("\n\n * Le systeme est : ");
aff_syst(a,b,n);
}

switch(valider)
{
case 'A'  : remplir();cramer(a,b,n);                break;
case 'B'  : remplir();inversion(a,b,n);             break;
case 'C'  : remplir();gauss(a,b,n);                 break;
case 'D'  : remplir();gauss_pivot_partiel(a,b,n);   break;
case 'E'  : remplir();gauss_pivot_total(a,b,n);     break;
case 'F'  : remplir();jordan(a,b,n);                break;
case 'G'  : remplir();LU_doolittle(a,b,n);          break; 
case 'H'  : remplir();LU_crout(a,b,n);              break; 
case 'I'  : remplir();cholesky(a,b,n);              break;
case 'J'  : remplir();jacobie(a,b,n);               break;
case 'K'  : remplir();gauss_seidel(a,b,n);          break;
case 'X'  : printf("\nProgramme interompu par l'utilisateur.......\n\n"); 
            system("PAUSE");exit(0);                break; 
default   : goto reprendre_choix;
}
printf("\n");
system("PAUSE");
main();
}

// -------------------------------------------------------------
//  Fonctions utiles pour la résolution par Cramer et inversion
// -------------------------------------------------------------

// calcul du determinant

float det(float a[19][19],int n)
{
int k,j;float c[19][19],s;

k=n-1;

if(n==0) return(1);

s=0;
for(j=0;j<n;j++)
{
comatrices(a,c,k,j,n);
s=s+pow(-1,k+j)*a[k][j]*det(c,k);
}
return(s);
}

// calcul des comatrices

void comatrices(float a[19][19],float c[19][19],int i,int j,int n)
{
int l,k;
for(l=0;l<n;l++) for(k=0;k<n;k++)
{
if ((l<i)&&(k<j)) c[l][k]=a[l][k];
if ((l>i)&&(k<j)) c[l-1][k]=a[l][k];
if ((l<i)&&(k>j)) c[l][k-1]=a[l][k];
if ((l>i)&&(k>j)) c[l-1][k-1]=a[l][k];
}
}

// -------------------------------------------------------------
//                  Résolution Par Cramer
// -------------------------------------------------------------

void cramer(float a[19][19],float b[19],int n)
{
float A[19][19],x[19],deter;int i,j,k;

deter=det(a,n);

if (deter==0) 
{
printf("\n => Determinant nul, pas de solutions \n\n");
system("PAUSE");main(); 
}

for(j=0;j<n;j++) 
{
for(k=0;k<n;k++) 
{
if (k==j) for(i=0;i<n;i++) A[i][k]=b[i];
else for(i=0;i<n;i++) A[i][k]=a[i][k];
}                  
x[j]=det(A,n)/deter;
}
printf("\n-------------- Cramer -------------\n");
printf("\n * Le determiant du systeme : %f \n",deter);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//                  Résolution Par Inversion
// -------------------------------------------------------------

void inversion(float a[19][19],float b[19],int n)
{
float c[19][19],comat[19][19],a_1[19][19],x[19],deter,s;int i,j;

deter=det(a,n);

if (deter==0) 
{
printf("\n => Matrice non inversible, pas de solutions \n\n");
system("PAUSE");main(); 
}

for(i=0;i<n;i++) for(j=0;j<n;j++)
{
comatrices(a,c,i,j,n);
comat[i][j]=(pow(-1,i+j)/deter)*(det(c,n-1));
}

// transposée
for(i=0;i<n;i++) for(j=0;j<n;j++)
a_1[i][j]=comat[j][i];

//résolution
for(i=0;i<n;i++) 
{s=0;
for(j=0;j<n;j++) s=s+a_1[i][j]*b[j];
x[i]=s;
}

printf("\n-------------- Inversion -------------\n");
printf("\n * La matrice inverse A^-1 :");
aff_mat(a_1,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//                  Résolution Par Gauss
// -------------------------------------------------------------

void gauss(float a[19][19],float b[19],int n)

{
float x[19],p,s;int i,j,k;

for(k=0;k<n-1;k++)
{
if (a[k][k]==0) 
{
printf("\n\n * Un pivot nul ! => methode de Gauss non applicable\n\n");
system("PAUSE");main();
}                 

//réduction
for(i=k+1;i<n;i++)
{
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
}
}

//Résolution
for(i=n-1;i>=0;i--)
{
s=0;   
for(j=i+1;j<n;j++)s=s+a[i][j]*x[j];                    
x[i]=(b[i]-s)/a[i][i];
}
zero(a,b,n);
printf("\n-------------- Gauss -------------\n");
printf("\n * La matrice reduite :");
aff_syst(a,b,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//              Résolution Par Gauss pivot pariel
// -------------------------------------------------------------

void gauss_pivot_partiel(float a[19][19],float b[19],int n)

{
float x[19],p,s,ref,temp;int i,j,k,ligne;     

for(k=0;k<n-1;k++)
{
// max pour le pivot partiel                
ref=0;
for(i=k;i<n;i++) if(fabs(a[i][k])>ref) {ref=fabs(a[i][k]);ligne=i;}

// pivotations
for(j=k;j<n;j++) 
{temp=a[k][j]; a[k][j]=a[ligne][j] ;a[ligne][j]=temp;}

temp=b[k]; b[k]=b[ligne]; b[ligne]=temp;

if (a[k][k]==0) 
{
printf("\n\n* Un pivot nul ! => methode de Gauss pivot partiel non applicable\n\n");
system("PAUSE");main();
}                 


//réduction                
for(i=k+1;i<n;i++)
{
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
}
}

//Résolution
for(i=n-1;i>=0;i--)
{
s=0;   
for(j=i+1;j<n;j++) s=s+a[i][j]*x[j];                    
x[i]=(b[i]-s)/a[i][i];
}
zero(a,b,n);
printf("\n-------- Gauss avec pivot partiel --------\n");
printf("\n * La matrice reduite :");
aff_syst(a,b,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
printf("\n");
}

// -------------------------------------------------------------
//              Résolution Par Gauss pivot total
// -------------------------------------------------------------

void gauss_pivot_total(float a[19][19],float b[19],int n)

{
float x[19],p,s,ref,temp;int i,j,k,ligne,colonne,pivot_sol[19],temps;

// vecteur de pivotation des solutions
for(i=0;i<n;i++) pivot_sol[i]=i;

for(k=0;k<n-1;k++)
{                  
// max pour le pivot total                
ref=0;
for(i=k;i<n;i++) for (j=k;j<n;j++)
if(fabs(a[i][j])>ref) {ref=fabs(a[i][j]);ligne=i;colonne=j;}

// pivotations
for(j=k;j<n;j++) {temp=a[k][j]; a[k][j]=a[ligne][j] ;a[ligne][j]=temp;}

temp=b[k]; b[k]=b[ligne]; b[ligne]=temp;

for(i=0;i<n;i++) {temp=a[i][k]; a[i][k]=a[i][colonne] ;a[i][colonne]=temp;}

// remplissage du vecteur accordé aux pivotations 
temps=pivot_sol[k];
pivot_sol[k]=pivot_sol[colonne];
pivot_sol[colonne]=temps;

if (a[k][k]==0) 
{
printf("\n\n * Un pivot nul ! => methode de Gauss pivot total non applicable\n\n");
system("PAUSE");main();
}                 

//réduction                
for(i=k+1;i<n;i++)
{
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
}
}

// Résolution
for(i=n-1;i>=0;i--)
{   
s=0;    
for(j=i+1;j<n;j++)s=s+a[i][j]*b[j];                    
b[i]=(b[i]-s)/a[i][i];
}

// pivotation des solutions
for(i=0;i<n;i++) x[pivot_sol[i]]=b[i];

zero(a,b,n);
printf("\n--------- Gauss avec pivot total ---------\n");
printf("\n * La matrice reduite :");
aff_syst(a,b,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
printf("\n");
}

// -------------------------------------------------------------
//                  Résolution Par Jordan
// -------------------------------------------------------------

void jordan(float a[19][19],float b[19],int n)

{
float p;int i,j,k;

for(k=0;k<n;k++)
{
if (a[k][k]==0)
{
printf("\n\n * Un pivot nul ! => methode de Jordan non applicable\n\n");
system("PAUSE");main();
}                 

p=a[k][k]; 

//normalisation               
for (j=k;j<n;j++) a[k][j]=a[k][j]/p;
b[k]=b[k]/p;

//réduction                
for(i=0;i<n;i++)
{
 if (i!=k)
 {
  p=a[i][k];                
  for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
  b[i]=b[i]-p*b[k];
 }
}
}
zero(a,b,n);
printf("\n-------------- Jordan -------------\n");
printf("\n * La matrice reduite :");
aff_syst(a,b,n);
printf("\n * La resolution donne :\n\n");
for(i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,b[i]);
printf("\n");
}

// -------------------------------------------------------------
//       Résolution Par décomposition LU Doolittle
// -------------------------------------------------------------

void LU_doolittle(float a[19][19],float b[19],int n)

{
float L[19][19],U[19][19],x[19],y[19],s;int i,j,k,m;

for (i=0;i<n;i++) for (j=0;j<n;j++) 
{
if(i==j) L[i][j]=1; else L[i][j]=0;
U[i][j]=0;
}

for (m=0;m<n;m++)
{       
for (j=m;j<n;j++)   
{
s=0;     
for (k=0;k<m;k++) s=s+L[m][k]*U[k][j];     
U[m][j]=a[m][j]-s;
}
if (U[k][k]==0) 
{
printf("\n\n * Un mineur nul ! => methode de LU non applicable\n\n");
system("PAUSE");main();
}                 

for (i=m+1;i<n;i++) 
{
s=0;   
for (k=0;k<m;k++) s=s+L[i][k]*U[k][m];    
L[i][m]=(a[i][m]-s)/U[m][m];
}
}

// resolution
for(i=0;i<n;i++)
{
s=0;                      
for(j=0;j<i;j++) s=s+L[i][j]*y[j];                    
y[i]=(b[i]-s)/L[i][i];
}

for(i=n-1;i>=0;i--)
{
s=0;                      
for(j=i+1;j<n;j++)s=s+U[i][j]*x[j];                    
x[i]=(y[i]-s)/U[i][i];
}

printf("\n------------ LU Doolittle -----------\n");
printf("\n * A = L * U \n");
printf("\n * La matrice L :");
aff_mat(L,n);
printf("\n * La matrice U :");
aff_mat(U,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//       Résolution Par décomposition LU Crout
// -------------------------------------------------------------

void LU_crout(float a[19][19],float b[19],int n)

{
float L[19][19],U[19][19],x[19],y[19],s;int i,j,k,m;

for (i=0;i<n;i++) for (j=0;j<n;j++) 
{
if(i==j) U[i][j]=1; else U[i][j]=0;
L[i][j]=0;
}

for (m=0;m<n;m++)
{    
for (i=m;i<n;i++) 
{
s=0;    
for (k=0;k<m;k++) s=s+L[i][k]*U[k][m];    
L[i][m]=a[i][m]-s;
}

if (L[k][k]==0) 
{
printf("\n\n* Un mineur nul ! => methode de LU non applicable\n\n");
system("PAUSE");main();
}                 

for (j=m+1;j<n;j++)   
{
s=0;    
for (k=0;k<m;k++) s=s+L[m][k]*U[k][j];     
U[m][j]=(a[m][j]-s)/L[m][m];
}
}

// resolution
for(i=0;i<n;i++)
{
s=0;                      
for(j=0;j<i;j++) s=s+L[i][j]*y[j];                    
y[i]=(b[i]-s)/L[i][i];
}

for(i=n-1;i>=0;i--)
{
s=0;                      
for(j=i+1;j<n;j++)s=s+U[i][j]*x[j];                    
x[i]=(y[i]-s)/U[i][i];
}

printf("\n-------------- LU Crout -------------\n");
printf("\n * A = L * U \n");
printf("\n * La matrice L :");
aff_mat(L,n);
printf("\n * La matrice U :");
aff_mat(U,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//       Résolution Par décomposition Cholesky
// -------------------------------------------------------------

void cholesky(float a[19][19],float b[19],int n)
{
float L[19][19],Lt[19][19],x[19],y[19],s,p;int i,j,k;

// vérification de le symétrie
for (i=0;i<n;i++) for (j=0;j<n;j++) 
if (a[i][j]!=a[j][i])
{
printf("\n\n * Non symetrique => methode de Cholesky non applicable\n\n");
system("PAUSE");main();
}

for (i=0;i<n;i++) for (j=0;j<n;j++) L[i][j]=0;

for (i=0;i<n;i++) 
{
s=0;
for (k=0;k<i;k++) s=s+pow(L[i][k],2);
p=a[i][i]-s;

if (p<=0)
{
printf("\n\n * Non definie positive => methode de Cholesky non applicable\n\n");
system("PAUSE");main();
}

L[i][i]=sqrt(p);

for(j=i+1;j<n;j++)
{
s=0;
for (k=0;k<i;k++) s=s+L[i][k]*L[j][k];
L[j][i]=(a[j][i]-s)/L[i][i];
}
}

for (i=0;i<n;i++) for (j=0;j<n;j++) Lt[i][j]=L[j][i];

// resolution
for(i=0;i<n;i++)
{
s=0;                      
for(j=0;j<i;j++) s=s+L[i][j]*y[j];                    
y[i]=(b[i]-s)/L[i][i];
}

for(i=n-1;i>=0;i--)
{
s=0;                      
for(j=i+1;j<n;j++) s=s+Lt[i][j]*x[j];                    
x[i]=(y[i]-s)/Lt[i][i];
}

printf("\n--------------- Cholesky --------------\n");
printf("\n * A = L * Lt \n");
printf("\n * La matrice L :");
aff_mat(L,n);
printf("\n * La matrice Lt :");
aff_mat(Lt,n);
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x[i]);
}

// -------------------------------------------------------------
//         Résolution Par Jacobie
// -------------------------------------------------------------

void jacobie(float a[19][19],float b[19],int n)
{
float x[19],x1[19],x2[19],s,eps=1e-4;int i,j,k,iter=0;     

//initialisation du vecteur
printf("\n * Veuillez initialisez le vecteur solution : \n\n");
for (i=0;i<n;i++) {printf(" X(0)[%d]= ",i+1);scanf("%f",&x1[i]);} 

do
{      
for(i=0;i<n;i++)
{
s=0;
for (j=0;j<n;j++) if (i!=j) s=s+a[i][j]*x1[j];
x2[i]=(b[i]-s)/a[i][i];
}
for (k=0;k<n;k++) {x[k]=fabs(x1[k]-x2[k]);x1[k]=x2[k];}  
 
iter++;
}
while (norme(x,n)>eps) ;

printf("\n-------------- Jacobie -------------\n");
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x2[i]);
printf("\n * Apres %d iteration, precision 10^-4.  \n",iter);
}

// -------------------------------------------------------------
//        Résolution Par Gauss-Seidel
// -------------------------------------------------------------

void gauss_seidel(float a[19][19],float b[19],int n)
{
float x[19],x1[19],x2[19],s,p,eps=1e-4;int i,j,k,iter=0;     

//initialisation du vecteur
printf("\n * Veuillez initialisez le vecteur solution : \n\n");
for (i=0;i<n;i++) {printf(" X(0)[%d]= ",i+1);scanf("%f",&x1[i]);} 

do
{  
for(i=0;i<n;i++)
{
s=0;
p=0;
for (j=i+1;j<n;j++) s=s+a[i][j]*x1[j];
for (j=0;j<i;j++) p=p+a[i][j]*x2[j];
x2[i]=(b[i]-s-p)/a[i][i];
}
for (k=0;k<n;k++) {x[k]=fabs(x1[k]-x2[k]);x1[k]=x2[k];}  

iter++;
}
while (norme(x,n)>eps) ;

printf("\n-------------- Gauss-Saidel -------------\n");
printf("\n * La resolution donne :\n\n");
for (i=0;i<n;i++) printf(" X_%d =  %f  ;\n",i+1,x2[i]);
printf("\n * Apres %d iteration, precision 10^-4.  \n",iter);
}

// Calcul norme vecteur

float norme(float x[19],int n)
{
float ref;int i;
ref=0;
for(i=0;i<n;i++) if (x[i]>ref) ref=x[i];
return(ref);
}

// fonction de remplissage

void remplissage(float a[19][19],float b[19],int n)
{
int i,j;     
printf("\n * La matrice A:\n\n");
for (i=0;i<n;i++) 
{
printf(" | ");
for (j=0;j<n;j++) scanf("%f",&a[i][j]);
}
printf("\n * Le vecteur B:\n\n");
for(i=0;i<n;i++) 
{
printf(" | ");
scanf("%f",&b[i]);
}
}

// fonction d'affichage système

void aff_syst(float a[19][19],float b[19],int n)
{
int i,j;     
printf("\n\n");
for (i=0;i<n;i++)
{
printf(" [");
for (j=0;j<n;j++)
{
printf("  %.4f  ",a[i][j]);
}
printf("] [  %.4f  ]",b[i]);
printf("\n");
}
}

// fonction d'affichage matrice

void aff_mat(float a[19][19],int n)
{
int i,j;     
printf("\n\n");
for (i=0;i<n;i++)
{
printf(" [");
for (j=0;j<n;j++)
{
printf("  %.4f  ",a[i][j]);
}
printf("]\n");
}
}

// Mettre à Zero les elements qui doivent etre des zéro

void zero(float a[19][19],float b[19],int n)
{
int i,j;float eps=1e-4;
for(i=0;i<n;i++)
{
for (j=0;j<n;j++) if (fabs(a[i][j])<eps) a[i][j]=0;      
if (fabs(b[i])<eps) b[i]=0;
}
}