Ce script 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 d'ingénieur de SFAX
Source / Exemple :
// mini projet ; le 20/06/2008
#include <stdio.h>
#include <math.h>
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;
}
}
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];
}
}
// 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);
}
// -------------------------------------------------------------
// 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");
}
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");
}
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");
}
//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");
}
//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");
}
//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");
}
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");
}
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");
}
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");
}
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");
}
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);
}
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);
switch(valider)
{
case 'A' : {
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);
}cramer(a,b,n); break;
case 'B' :{
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);
inversion(a,b,n); break;
}
case 'C' :{
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);
gauss(a,b,n); break;
}
case 'D' :{
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);
gauss_pivot_partiel(a,b,n); break;
}
case 'E' :{
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);
gauss_pivot_total(a,b,n); break; }
case 'F' :{
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);
jordan(a,b,n); break; }
case 'G' :{
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);
LU_doolittle(a,b,n); break; }
case 'H' :{
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);
LU_crout(a,b,n); break; }
case 'I' :{
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);
cholesky(a,b,n); break;
}
case 'J' :{
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);
jacobie(a,b,n); break;
}
case 'K' :{
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);
gauss_seidel(a,b,n); break; }
case 'X' :{
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);
printf("\nProgramme interompu par l'utilisateur.......\n\n");
break; }
default : goto reprendre_choix;
}
printf("\n");
}
// kolsi ahmed
// ahmedkolsi@gmail.com
Conclusion :
il contient des fonctions d'analyse numérique
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.