Ce document contient :
résolution des systèmes Ax=B par : inversion, cramer ,gauss, gauss pivot partiel, gauss pivot total, jordan et LU (croot et doolittle), Cholesky, jacobie, Gauss-Seidel.
Ecole nationale polytechnique d'alger, ENP
compilé sous Dev-C++
Source / Exemple :
/* ========================================================================== */
/* 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;
}
}
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.