Résolution des système ax=b

Soyez le premier à donner votre avis sur cette source.

Snippet vu 33 294 fois - Téléchargée 24 fois

Contenu du snippet

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; } }

A voir également

Ajouter un commentaire Commentaires
cresuscresus
Messages postés
2
Date d'inscription
lundi 25 décembre 2006
Statut
Membre
Dernière intervention
22 novembre 2008

14 déc. 2008 à 19:06
Pas de quoi, heureux d'avoir pu t'être utile
Uarel
Messages postés
1
Date d'inscription
lundi 15 septembre 2008
Statut
Membre
Dernière intervention
14 décembre 2008

14 déc. 2008 à 14:30
Merci grandement ça m'a énormément aidé !!
cs_zizou23
Messages postés
3
Date d'inscription
lundi 28 janvier 2008
Statut
Membre
Dernière intervention
14 mars 2008

14 mars 2008 à 02:57
merci pour cette idée pour m'aider a engager en Tp mais c tu veux ,je veux la solution avec c++builder,....j'attend ta reponse,MERCI
theg56
Messages postés
4
Date d'inscription
vendredi 10 décembre 2004
Statut
Membre
Dernière intervention
29 janvier 2008

29 janv. 2008 à 12:01
c'est déja du c
hissou
Messages postés
5
Date d'inscription
dimanche 19 juin 2005
Statut
Membre
Dernière intervention
28 janvier 2008

28 janv. 2008 à 23:35
est ce que c possible de me donner la solution d'un systeme d'eqautions lineaires AX=B a l'aide de la methode de Gauss seulement en C. c'est urgent
Afficher les 11 commentaires

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.