Resolution de systeme d equation methode iteratives jacobi gauss seidel en langage c

resolution de systeme d equation methode iteratives jacobi gauss seidel en langage c    
 Gauss–Seidel method

Le terme Gauss–Seidel method est cité dans le Wikipedia de langue anglaise. Il est défini comme suit:

In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and is similar to the Jacobi method. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either diagonally dominant, orsymmetric and positive definite. It was only mentioned in a private letter from Gauss to his student Gerling in 1823. A publication was not delivered before 1874 by Seidel.

Ceci est un extrait de l'article Gauss–Seidel method de l'encyclopédie libre Wikipedia. La liste des auteurs est disponible sur Wikipedia.
           
Méthode de Gauss-Seidel
Le terme Méthode de Gauss-Seidel du Wikipedia en langue allemande correspond au terme Gauss–Seidel method issu du Wikipedia en langue anglaise. Il est défini comme suit :

La méthode de Gauss-Seidel est une méthode itérative de résolution d'un système linéaire (de dimension finie) de la forme Ax = b, ce qui signifie qu'elle génère une suite qui converge vers une solution de cette équation, lorsque celle-ci en a une et lorsque des conditions de convergence sont satisfaites (par exemple lorsque A est symétrique définie positive). L'algorithme suppose que la diagonale de A est formée d'éléments non nuls.
La méthode se décline en une version « par blocs ».


Ceci est un extrait de l'article Méthode de Gauss-Seidel de l'encyclopédie libre Wikipedia. La liste des auteurs est disponible sur Wikipedia
#include <math.h>
#include <conio.h>
#include <stdio.h>
#include <stdlib.h>
//declaration de contantes -------------------------------------------------------------------------------------------------
#define NMAX 20 //Taille maximum d'une matrice
#define N2MAX 21 //Taille max +1
#define ITERMAX 21
//declaration de variables globales (je sais, c'est pas beau, mais bon!) ---------------------------------------------------
float ident[NMAX][NMAX]; //Matrice identité
float a[NMAX][NMAX];
float diag[NMAX][NMAX]; //Matrice diagonale de a
float triInf[NMAX][NMAX]; //Triangle inférieur + diagonale de a
float t[NMAX][ITERMAX];
float x[NMAX];
float inversee[NMAX][NMAX]; //Matrice resultat de l'inversion
float * bi;
float solu[NMAX];
float * u; //Vecteur unitaire colonne de la resolution gaussienne
float eps;
float w=1; // precision souhaitee pour la resolution
float erreur[NMAX];
int iter; // nombre d'iterations max demandee pour la resolution
int n,k,i;
int P=1;
int temps=0; //N=nombre de ligne de la matrice inversee, M=la dim commune aus deux matrice, P=le nbre de colonne de u
int choix;
float matresmult[NMAX];
float resmult[NMAX];
float matiter[NMAX][ITERMAX];
float error = 1.0;
//declaration des fonctions ------------------------------------------------------------------------------------------------
//Fonction d'initialisation ------------------------------------------------------------------------------------------------
void init(){
printf("\t============================================\n\n");
printf("\t RESOLUTION DE SYSTEME LINEAIRE\n\n");
printf("\t PAR METHODES ITERATIVES\n\n");
printf("\tCHAFFARDON P, GERVAIS B, MAILLAND F, SOL P\n\n");
printf("\t============================================\n\n\n\n");
printf("Quelle methode desirez-vous utiliser?\n");
printf("1- Methode de Jacobi\n");
printf("2- Methode de Gauss-Seidel\n\n");
do{
printf("Entrez votre choix (1, 2)");
scanf("%d", &choix);
getchar();
printf("choix: %d\n\n\n", choix);
}while(choix != 1 && choix != 2);
}
//Fonction recuperant la diagonale de la matrice A -------------------------------------------------------------------------
void recupDiagMatrice(){
for(int i = 1; i <= n; i++){
diag[i][i] = a[i][i];
}
}
//Fonction calculant la triangulaire inférieure de la matrice passée en paramètre ------------------------------------------
void recupTriInf(){
for(int i = 1; i <= n; i++){
for(int j = 1; j <= i; j++){
printf("a[i](j] = %f\n", a[i][j]);
triInf[i][j] = a[i][j];
}
}
}
//Initialisation de la matrice résultante de la multiplication ------------------------------------------------------------
void iniresmult()
{
for (int i=1;i<=n;i++) //resmult aura bien N lignes
{ //et P colonnes
matresmult[i]=0;
}
}
//Initialisation de la matrice résultante de la multiplication resmult------------------------------------------------------------
void iniresmult2()
{
for (int i=1;i<=n;i++) //resmult aura bien N lignes
{ //et P colonnes
resmult[i]=0;
}
}
//Calcul de la matrice résultante -------------------------------------------------------------------------------------------
void calcresmult()
{
for (int i=1;i<=n;i++)
{
for (int k=1;k<=n;k++)
{
matresmult[i]=inversee[i][k]*u[k]+matresmult[i];
}
}
}
//Affichage de la matrice résultante ---------------------------------------------------------------------------------------
void affresmult()
{
printf ("\nProduit matriciel :\n");
for (int i=1;i<=n;i++)
{
for (int j=1;j<=P;j++)
{
printf ("ligne %d,colonne %d : ",i,j);
printf ("%f\n",matresmult[i]);
}
}
}
//Multiplication d'une matrice n*n par une matrice n*1 -------------------------------------------------------------------------------------------
void calcmult(float bob[NMAX][NMAX],float bib[NMAX])
{
for (int i=1;i<=n;i++)
{
for (int k=1;k<=n;k++)
{
resmult[i]=bob[i][k]*bib[k]+resmult[i];
}
}
}
//Fonction d'inversion d'une matrice par la methode du pivot de gauss ------------------------------------------------------
int pivot(float g[NMAX][NMAX],float b[NMAX])
{
int i,j,k,l,err;
float pivot,coef,s;
float t[NMAX][N2MAX];
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
t[i][j]=g[i][j];
t[i][n+1]=b[i];
}
err=1;
k=1;
while (err==1 && k<n)
{
pivot=fabs(t[k][k]);
l=k;
while(pivot==0 && l<n)
{
l++;
pivot=fabs(t[l][k]);
}
if(pivot!=0)
{
if(l!=k)
for(j=k;j<=n+1;j++)
{
pivot=t[k][j];
t[k][j]=t[l][j];
t[l][j]=pivot;
}
pivot=t[k][k];
for(i=k+1;i<=n;i++)
{
coef=t[i][k]/pivot;
for(j=1;j<=n+1;j++)
t[i][j] -= coef*t[k][j];
}
//printf(" A(%d) b(%d)\n",k+1,k+1);
/*for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
printf("%f ",t[i][j]);
printf(" I %f \n",t[i][n+1]);
}*/
}
else err=0;
k++;
}
if(t[n][n]==0) err=0;
if(err==1)
{
b[n]=t[n][n+1]/t[n][n];
for(i=n-1;i>=1;i--)
{
s=t[i][n+1];
for(j=i+1;j<=n;j++)
s-=t[i][j]*b[j];
b[i]=s/t[i][i];
}
}
return(err);
}
//Fonction d'inversion de la matrice ------------------------------------------------------------------------------------------
void invMatrice(float g[NMAX][NMAX]){
int err;
bi = (float *) malloc (sizeof(float)*n); //allocation dyn de notre vecteur unitaire
/*printf(" A(1) b(1)\n");
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
{
printf("%f ",g[i][j]);
printf(" I %f \n",bi[i]); //Affichage du premier resultat de la methode
}
}*/
for (int z = 1;z <= n; z++) //indice de la colonne en cours d'etude
{
for (int y = 1;y <= n; y++)
{
bi[y] = ident[y][z]; //initialiser le vecteur identitaire correspondant à l'equation y
}
err=pivot(g,bi); //Appel de la fonction pivot
if(err==1)
{
//printf("Solution:\n"); //Si matrice non singuliere, affichage des resultats pas a pas
for(int i=1;i<=n;i++)
{
//printf("x%d =%f \n",i,bi[i]);
inversee[i][z] = bi[i]; //Formation de notre matrice inversee par addition des resultats intermediaires
}
}
else
{
printf("Erreur, matrice singulière\n"); //Si matrice singuliere, affichage d'erreur et sortie
getchar();
exit(0);
}
}
}
//Fonction de resolution de l'aproximation de x ----------------------------------------------------------------------------
void switching(){
switch(choix){
case 1:
printf("\n******************** Methode de Jacobi *************************\n");
//Calcul de l'aproximation de x par jacobi
recupDiagMatrice(); //recuperation de la matrice diagonale de A
invMatrice(diag); //inversion de la matrice
iniresmult(); //multiplication de la matrice inversee * bi
calcresmult();
//affresmult();
break;
case 2:
printf("\n********************* Methode de Gauss-Seidel ***************************\n");
//Calcul de l'aproximation de x par gauss-seidel
recupTriInf(); //recuperation de la triangulaire inférieure de A
invMatrice(triInf); //inversion de la matrice
iniresmult(); //multiplication de la matrice inversee * bi
calcresmult();
//affresmult();
break;
default:
break;
}
}
//Fonction d'affichage de la matrice inversee ------------------------------------------------------------------------------
void aff(){
printf("Affichage de la matrice Inversee\n");
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++){
if(j == n){
printf("%f\n", inversee[i][j]);
}else{
printf("%f ", inversee[i][j]);
}
}
}
}
//Fonction de creation d'une matrice identite de taille n -------------------------------------------------------------
void creatematriceident(){
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n; j++){
if(i == j){
ident[i][j] = 1;
}else{
ident[i][j] = 0;
}
}
}
}
// Fonction demandant les parametres du systeme a resoudre a l'utilisateur ---------------------------------------------
void param()
{
printf("Veuillez indiquer l'ordre de la matrice A:");
scanf("%d", &n);
getchar();
u = (float *) malloc (sizeof (float)*n);
printf("Veuillez saisir les elements de la matrice A (la saisie se fera dans l'ordre des lignes)\n");
for(int i = 1; i <= n; i++)
{
for(int j = 1; j <= n; j++)
{
printf("element[%d][%d]= ", i, j);
scanf("%f", &a[i][j]);
getchar();
}
}
printf("Veuillez saisir le vecteur U:");
for(int k = 1; k <= n; k ++)
{
printf("element[%d]= ", k);
scanf("%f", &u[k]);
getchar();
}
printf("Veuillez saisir le nombre maximum d'iteration:");
scanf("%d", &iter);
getchar();
printf("Veuillez saisir la precision:");
scanf("%f", &eps);
getchar();
printf("*********** Vous avez finit les saisies, verifions les ensemble ***********\n\n");
printf("\nAffichage de la matrice A\n");
for(i = 1; i <= n; i++)
{
for(int j = 1; j <= n; j++)
{
if(j == n)
{
printf("%f\n", a[i][j]);
}else
{
printf("%f ", a[i][j]);
}
}
}
printf("\nAffichage du vecteur U \n");
for(int l = 1; l <= n; l++)
{
printf("%f\n", u[l]);
}
printf("\nAffichage de iter = %d\n", iter);
printf("Affichage de eps = %f\n\n", eps);
printf("****** Une fois ces verifications terminees, passons a la resolution! ******\n\n");
}
/** Fonction calculant la norme vectoriel euclidienne ------------------------------------------------------*/
float al_norme_vect(float x[NMAX],int n)
{
int i;
w=0;
for(i=1;i<=n;i++)
w += x[i]*x[i];
w=sqrt(w);
return(w);
}
/** Fonction effectuant la résolution du systeme par la méthode de jacobi --------------------------------------------------*/
void resolutionj(float p[NMAX][NMAX],float d[NMAX],float solu[NMAX],int n,int iter,float eps)
{
int i,j;
float alfa,s;
double x2[NMAX];
for(i=1;i<=n;i++)
for(j=1;j<=ITERMAX-1;j++)
matiter[i][j]=0,0;
alfa=1.0;
k=1;
while(k<=iter)
{
for(i=1;i<=n;i++)
{
s=d[i];
for(j=1;j<=n;j++)
if(i!=j)
s -= p[i][j]*solu[j];
x2[i]=s/p[i][i];
}
alfa=0.0;
for(i=1;i<=n;i++)
{
alfa+=pow(x2[i]-solu[i],2);
solu[i]=x2[i];
}
for(i=1;i<=n;i++)
{
t[i][k]=solu[i];
}
k++;
}
}
/** Fonction effectuant la résolution du systeme par la méthode de gauss-seidel --------------------------------------------------*/
void resolutiong(float p[NMAX][NMAX],float d[NMAX],float solu[NMAX],int n,int iter,float eps)
{
int i,j;
float alfa,s,ddd;
for(i=1;i<=n;i++)
for(j=1;j<=ITERMAX-1;j++)
matiter[i][j]=0,0;
alfa=1.0;
k=1;
while(k<=iter)
{
alfa=0.0;
for(i=1;i<=n;i++)
{
ddd=solu[i];
s=d[i];
for(j=1;j<=n;j++)
if(i!=j)
s -= p[i][j]*solu[j];
solu[i]=s/p[i][i];
ddd-=solu[i];
alfa+=ddd*ddd;
}
for(i=1;i<=n;i++)
{
t[i][k]=solu[i];
}
k++;
}
}
/* Fonction d'affichage de la solution et des differents parametres de la resolution --------------------------------------*/
void affichsolgenerale()
{
float temp[NMAX];
float temp2[NMAX];
float fin[NMAX];
int cc,aaa;
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
//printf("%f ",a[i][j]);
//printf(" I %f\n",u[i]);
x[i]=0;
}
//printf("Solution exacte :\n");
printf("\nVoici les n iterations demandees, a chaque fois sont proposees: \n");
printf("Le vecteur X(n) / Err, la norme de la difference entre l'itere X(n) et la solution exacte \n");
printf("et L'Erreur commise sur la solution du systeme... \n\n");
for(i=1;i<=n;i++)
//printf("%f\n",matresmult[i]);
for(int o=1; o<=n;o++)
{
temp[o] = matresmult[o];
}
if (choix == 1)
{
resolutionj(a,u,temp,n,iter,eps);
}
else
{
resolutiong(a,u,temp,n,iter,eps);
}
for(int j=1;j<=(k-1);j++)
{
printf("x(%2d)= (",j);
for(i=1;i<=n;i++)
{
x[i]=t[i][j]-matresmult[i];
printf("%f ",t[i][j]);
temp2[i] = t[i][j];
}
iniresmult2();
calcmult(a,temp2);
for(int dd=1;dd<=n;dd++)
{
erreur[dd]=u[dd]-resmult[dd];
}
printf(" err= %f",al_norme_vect(x,n));
printf(")erreur= %f\n",al_norme_vect(erreur,n));
error = al_norme_vect(erreur,n);
if (error<eps)
{
printf("La precision est atteinte !!\n Appuyer sur une touche pour avoir les resultats attendus\n");
getchar();
break;
}
temps++;
}
printf("\n\nLa solution estimee pour X est donc le dernier itere:\n");
for(i=1;i<=n;i++)
printf("%f\n",t[i][(k-1)]);
printf("\nLe resutat de l'equation avec cette solution etant de :\n");
for (aaa=1;aaa<=n;aaa++)
{
fin[aaa] = t[aaa][(k-1)];
}
iniresmult2();
calcmult(a,fin);
for (aaa=1;aaa<=n;aaa++)
{
printf("%f\n",resmult[aaa]);
}
printf("\n\nLe nombre d'iteration effective ayant ete de :%i\n",(temps+1));
printf("\n******* Merci d'avoir utilise le resolveur magique... !!!*******");
printf("\n******* Taper sur une touche pour quitter l'application *******");
getchar();
}
//fonction principale -------------------------------------------------------------------------------------------------------
main()
{
int i,j,err;
init(); //Choix de la methode de resolution
param();
creatematriceident(); //Creation de la matrice identite
switching(); //aproximation de x suivant la methode desiree
affichsolgenerale();
getchar();
}
view raw Resolution de systeme d equation methode iteratives jacobi gauss seidel en langage c hosted with ❤ by GitHub