Résolution de systèmes d'équations

Soyez le premier à donner votre avis sur cette source.

Snippet vu 13 506 fois - Téléchargée 29 fois

Contenu du snippet

Petite application permettant de résoudre des sytèmes d'équations de grande taille. Plusieurs méthodes de résolution sont proposées

Source / Exemple :

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

const int valmax=10;   // variable globale (constante) taille maximale des matrices et vecteurs

/*_________________________________________________________________________________________________________________________________________

                              PROCÉDURES ÉLÉMENTAIRES DE CALCUL SUR LES MATRICES ET LES VECTEURS
     ________________________________________________________________________________________________________________________________________*/
 
 
 

void lectureVect(int dim,double vect[]) //Permet à l'utilisateur de remplir un vecteur
{
  int i;
  for(i=0; i<dim;i++)
    {
      scanf("%lf",&vect[i]);
    }
}

void afficheVect(int dim,double vect[]) //Permet d'afficher à l'écran un vecteur
{
  int i;
  for(i=0;i<dim;i++)
    {
      printf("%lf\n",vect[i]);
    }
}

void lectureMat(int dim, double mat[][valmax]) //permet à l'utilisateur de remplir une matrice ligne par ligne
{
  int i,j; //i->ligne, j->colonne;
  for(i=0;i<dim;i++)
    {
      for(j=0;j<dim;j++)
 {
   scanf("%lf",&mat[i][j]);
 }
    }
}

void afficheMat(int dim, double mat[][valmax]) //Permet d'afficher à l'écran une matrice
{
  int i,j;
  for(i=0;i<dim;i++)
    {
      for(j=0;j<dim;j++)
 {
   printf("%lf ",mat[i][j]);
 }
      printf("\n");
    }
}

void sommeVect(int dim, double vect1[], double vect2[], double svect[])  //Calcule svect le vecteur somme des vecteurs vect1 et vect2
{
  int i;
  for(i=0;i<dim;i++)
    {
      svect[i]=vect1[i]+vect2[i];
    }
}
 
 

void diffVect(int dim, double vect1[], double vect2[], double svect[])  //Calcule svect le vecteur différence des vecteurs vect1 et vect2
{
  int i;
  for(i=0;i<dim;i++)
    {
      svect[i]=vect1[i]-vect2[i];
    }
}
 

void produitMatVect(int dimvect, int lmat, double vect[], double mat[][valmax],double pvect[]) //Calcule le vecteur pvect produit de la
{                                                        // matrice mat et du vecteur vect
  int i,j;
  for(i=0;i<lmat;i++)   //Change les lignes
    {
      pvect[i]=0;  //initialise l'élément avant d'effectuer la somme
      for(j=0;j<dimvect;j++)
 {
   pvect[i]=mat[i][j]*vect[j]+pvect[i];    // Somme les éléments d'une ligne de la matrice avec les éléments du vecteur.
 }
    }
}

double normeVect(int dim, double vect[]) // Renvoie la norme d'un vecteur
{
  int i;
  double norme=0;
  for(i=0;i<dim;i++)
    {
      norme=vect[i]*vect[i]+norme; //Somme des carrés des composantes du vecteur dans la variable norme
    }
  norme=sqrt(norme);  //calcul de la racine carré de la variable norme, la valeur dans norme est maintenant la norme du vecteur
  return norme;
}
 
 

void produitMat(int dim,double mat1[][valmax],double mat2[][valmax], double pmat[valmax][valmax])//Calcule le produit de deux matrices carrées
{                           // de même dimension (mat1*mat2) et met le résultat dans la matrice pmat
  int i,j,k=0;
  for (k=0;k<dim;k++)  // change la colonne de remplissage de pmat
    {
      for(i=0;i<dim;i++) //change la ligne de remplissage de pmat
 {
   pmat[i][k]=0;
   for(j=0;j<dim;j++)pmat[i][k]=mat1[i][j]*mat2[j][k]+pmat[i][k]; //somme des produits de chaque élément de la ligne de la mat1
 }                //avec les éléments de la colonne de mat2
    }
}

/*_________________________________________________________________________________________________________________________________________*/
/*
                                                    PROCÉDURES DE RÉSOLUTION

_________________________________________________________________________________________________________________________________________*/
 
 
 
 

void resolsys(int dim,double A[][valmax],double b[],double x[])//Résolution d'un système triangulaire supérieur
{
  int i,k,n=dim-1;
  double S;
  for(i=0; i<dim; i++)x[i]=0; //initialisation du vecteur x (solution)
 
  for(i=0;i<dim;i++)// fixe la ligne pour laquelle on cherche la solution (on commence par le bas)
    {
      S=0; // initialise S=0 pour effectuer la somme
 
      for(k=0;k<i;k++)
 {
   S=A[n-i][n-k]*x[n-k]+S;  // Calcul la somme des produits des éléments de x par leur coéfficient (élément de A)
 }
      x[n-i]=(b[n-i]-S)/A[n-i][n-i];// résolution de A[n-i][n-i]*x[n-i] + S = b[n-i]
    }
}

void resolsys2(int dim, double A[][valmax],double b[], double x[])//Résoloution d'un système triangulaire inférieur
{
  int i,k;// i=ligne;k=colonne
  double S;
  for(i=0; i<dim; i++)x[i]=0;//initialisation du vecteur x
  for(i=0;i<dim;i++)
    {
      S=0;
      for(k=0;k<i;k++)
 {
   S=A[i][k]*x[k]+S;//Calcul de la somme des produits des éléments de x par leur coefficient
 }
      x[i]=(b[i]-S)/A[i][i];// résolution de A[i][i]*x[i] + S = b [i]
    }
}
 
 

double maximum(int dim,int *i,double A[][valmax], int k)  // Retourne l'élément maximum de la colonne k d'une matrice à partir de la ligne k
{                                                       //et indique dans *i l'indice de cet élément
  int j;
  double max=abs(A[k][k]); //initialisation de max
  for(j=k+1;j<dim;j++)if(abs(A[j][k])>max)   // Si l'élément A[j][k] k étant fixé  est plus grand que l'élément précédent
    {
      max=abs(A[j][k]);                     //max vaut alors l'élément actuel
      *i=j;                                 //et l'indice est écrit dans *i
    }
  return max;
}

void changeligne(int dim, int j, int i,double A[][valmax],double b[]) //change deux lignes i et j d'une matrice et deux éléments d'un vecteur
{
  int k;
  double v;   // variable de transition
  for(k=0; k<dim; k++)
    {
      v=A[j][k];
      A[j][k]=A[i][k];//on intervertit les éléments colonnes par colonnes
      A[i][k]=v;
    }
  v=b[i];
  b[i]=b[j]; // on intervertit les éléments du vecteur
  b[j]=v;
}

void metzero(int dim,double A[][valmax], double b[],int k) // Applique le pivot situé en A[k][k]
{
  int i, j;
  double r, pivot=A[k][k];  // r est le rapport A[i][k]/pivot qui sera fixé pour chaque ligne
  for(i=k+1 ; i<dim ; i++)  //on effectue le calcul pour les lignes et colonnes inférieures  à l'élément pivot
    {
      r=A[i][k]/pivot; //calcul du rapport pour la ligne actuelle
      for(j=k; j<dim; j++)
 {
   A[i][j]=A[i][j]-A[k][j]*r; //on effectue Ligne[i]-r*Ligne[k]
 }
      b[i]=b[i]-b[k]*r;  //De même pour les solutions afin de garder l'équivalence entre le nouveau système et l'ancien
    }
}
 

int triangularise(int dim,double A[][valmax],double b[])//Permet de triangularisé le système Ax=b (A devient triangulaire supérieure)
{                               // la fonction retourne 0 le cas où il est impossible de calculer les solutions (déterminant nul)
  int j=0,i;
  double max=-1;
 
  while(j<dim && max!=0) //l'opération s'effectue sur toutes les colonnes et la fonction doit s'arrêter dans le cas où le déterminant est nul
    {                      // car dans ce cas il est impossible de trouver de solution
      if(A[j][j]==0)       //dans le cas où le pivot de la colonne est nul on lance la recherche du maximum qui renvoie le maximum
 {
   max=maximum(dim,&i,A,j);  // si tous les élément de la colonne sont nul max =0 et on sort du programme
   if(max!=0)changeligne(dim,j,i,A,b);  // Dans le cas ou il y a un maximum on intervertit les lignes
 }
      if(A[j][j]!=0)metzero(dim,A,b,j); // si le déterminant n'est pas nul,et le pivot n'est plus nul,on applique la fonction metzero
      j++; //on change de colonne
    }
  if(max==0) return 0;
  else return 1;
}
 
 

void separation(int dim,double E[][valmax], double F[][valmax], double A[][valmax])//séparation de A en E et F
{
  int i, j;//    si i>=j  E[i][j]=A[i][j] et F[i][j]=0;        si i<j  E[i][j]=0 et F[i][j]= A[i][j]
  for(i=0;i<dim; i++)
    {
      for(j=0;j<i+1;j++)
 {
   F[i][j]=0;
   E[i][j]=A[i][j];
 }
      for(j=i+1; j<dim; j++)
 {
   F[i][j]=-A[i][j];
   E[i][j]=0;
 }
    }
}
 

void systeme()//demande des données à l'utilisateur
{
  int val;
  double A[valmax][valmax];//matrica rentrée par l'utilisateur
  double b[valmax], x[valmax];//b->vecteur rentré par l'utilisateur      x-> vecteur solution du système
  printf("Nombre d'inconnues?\n");//le nombre d'inconnues sera égale à la dimension de la matrice et du vecteur
  scanf("%d",&val);
  printf("Saisie Matrice en ligne\n");
  lectureMat(val,A);
  printf("Saisie du vecteur(solution) \n");
  lectureVect(val,b);
  printf("\n");
  printf("Votre matrice:\n");//affiche la matrice rentrée
  afficheMat(val,A);
  printf("Votre vecteur:\n");//affiche le vecteur rentré
  afficheVect(val,b);
  printf("\n");
  if(triangularise(val, A, b))
    {
      resolsys(val, A,b,x);
      printf("Solution:\n");
      afficheVect(val,x);
    }
  else printf("Le déterminant est nul, le système n'est pas de Cramer\n");
}
 

void determinant()//donnée d'une matrice pour le calcul du déterminant
{
  int val,i;
  double det=1;
  double Mat[valmax][valmax];
  double t[valmax];
  printf("Taille de la matrice?\n");
  scanf("%d",&val);
  for(i=0; i<val ; i++)t[i]=0;
  printf("Saisie Matrice en ligne\n");
  lectureMat(val,Mat);
  printf("Votre Matrice:\n");
  afficheMat(val,Mat);
  printf("\n");
  if(triangularise(val,Mat,t))for(i=0;i<val;i++)det=det*Mat[i][i];//on applique la fonction triangularise pour que le déterminant soit
  else det=0;//                                                     égal au produit des éléments de la diagonale
  printf("déterminant=%lf\n",det);

}
 
 

void gauss_seidel()
{
  int val,it=30;
  double A[valmax][valmax];//matrice rentrée par l'utilisateur
  double b[valmax];//vecteur rentré par l'utilisateur
  double E[valmax][valmax];//décomposition inférieur de la matrice A
  double F[valmax][valmax];//décomposition suérieure de A
  double s[valmax];//vecteur somme de Fx+b
  double x[valmax];//initialement x0 puis prend la valeur de la solution trouvée
  double y[valmax];//vecteur solution du système
  double z[valmax];//vecteur différence entre les solutions successives
  double e;//e=précision
  int i=0,l=0;
  char choix=0;
  printf("Algorithme à %d itérations!\n",it);//nombre d'itération initialisé à 30
  printf("Précision de l'approximation?\n");//précision = différence entre deux termes consécutifs que l'on veut obtenir
  scanf("%lf",&e);
  printf("Nombre d'inconnues?\n");
  scanf("%d",&val);
  printf("Saisie Matrice en ligne\n");
  lectureMat(val,A);
  printf("Saisie du vecteur(solution) \n");
  lectureVect(val,b);
  printf("\n");
  printf("Votre matrice:\n");
  afficheMat(val,A);
  printf("Votre vecteur:\n");
  afficheVect(val,b);
  printf("\n");

  do
    {
      choix=0;
      for(i=0;i<val;i++)x[i]=1;//on initialise le vecteur x0 avec seulement des 1 puisque il est quelconque
      separation(val, E, F, A);//on sépare A en E et F
  do
    {
      produitMatVect(val,val,x, F,s);//on fait le produit de la matrice F et du vecteur x que l'on met dans le vecteur s
      sommeVect(val,s,b,s);//on somme le vecteur s et le vecteur b(rentré par l utilisateur)que l'on remet dans le vecteur s
      resolsys2(val,E,s,y);// on résout Ey=s (E étant triangulaire inférieure)
      diffVect(val,x,y,z);//on fait la différence entre y et x0
      l++;
      for(i=0;i<val;i++)x[i]=y[i];//on met le vecteur y dans le vecteur x0
    }while( normeVect(val,z)>e && l<it);//on effectue le programme tant que la différence entre 2 solutions successives est inférieure
                                       // à la précision demandée et tant que le nombre d'itérations n'est pas de 30
  if(l>=it)// on sort de la boucle parce que l'on a atteint 30 itérations
    {
      printf ("La suite n'a pas l'air de converger\n");
      printf ("Voulez vous recommencer en augmentant le nombre d'itération?(O/N)\n");//possibilité d'augmenter le nombre d'itérations
      scanf("\n %c",&choix);
      if(choix=='O' || choix=='o')
 {
   printf("Nombre d'itération?\n");
   scanf("\n %d",&it);
 }
    }
 
  else
    {
      printf("Solution:\n");//affiche le dernier terme calculé( ce sera une approximation du vecteur solution)
      afficheVect(val,y);
    }
    }while(choix=='o' || choix=='O');
}
 

void inverse()//calcul de l'inverse d' une matrice
{
  double A[valmax][valmax];//matrice rentrée
  double APrime[valmax][valmax];//recopie de A pour garder la matrice A intacte
  double AI[valmax][valmax];//matrice inverse de A
  double I[valmax][valmax];//matrice Identité
  double Id[valmax];//vectur avec un 1 à la ligne i et puis que des zéros (décomposition de la matrice Identité en vecteur colonne)
  double x[valmax];
  int k=0,i,val,pascramer=1,j;//i->ligne    j->colonne     k->numéro de la colonne de la matrice Identité concernée
  char choix=0;
  printf("Taille de la matrice?\n");
  scanf("%d",&val);
  printf("Saisie de la matrice:\n");
  lectureMat(val,A);
  printf("Votre matrice:\n");
  afficheMat(val,A);
  printf("\n");
  do
    {
      for(i=0; i<val;i++)Id[i]=0;//on initialise le vecteur à 0
      Id[k]=1;//on met un 1 en position i (ce qui forme la matrice Identité)
      for(i=0;i<val;i++)
 {
   for(j=0;j<val;j++)APrime[i][j]=A[i][j];//on recopie la matrice A dans A'
 }
      if(triangularise(val, APrime, Id))//on triangularise A' ( on effectue les mêmes opérations sur le vecteur)
 {
 resolsys(val, APrime,Id,x);//on résout pour chaque vecteur colonne de la matrice Identité modifiée le système A'x=Id
 }
 
      else
 {
   printf("Le déterminant est nul, le système n'est pas de Cramer, il est impossible de calculer la matrice inverse!\n");
   pascramer=0;
 }
      for(j=0;j<val;j++)AI[j][k]=x[j];//on remplit la matrice inverse de A où chaque colonne est le vecteur colonne solution du système
      k++;//on passe au vecteur colonne de la matrice Identité modifiée suivant
    }while(k<val && pascramer);
  if(pascramer)
    {
      printf("La matrice inverse trouvée:\n");
      afficheMat(val,AI);
      printf("Voulez-vous vérifier la matrice inverse?(O/N)\n");
      scanf("\n %c",&choix);
      if(choix=='o' || choix=='O')
 {
   printf("votre matrice:\n");
   afficheMat(val,A);
   printf("Matrice inverse calculée:\n");
   afficheMat(val, AI);
   printf("En effectuant le produit des matrices précédentes:\n");//vérification de la matrice inverse
   produitMat(val,A,AI,I);
   afficheMat(val,I);
 }
    }
}
 
 
 
 
 

main()//menu pour que l'utilisateur choisisse ce qu'il veut faire
{
 
  char choix;
 do
   {
     printf("\n");
     printf("1->Résolution d'un système linéaire par méthode du pivot de Gauss\n");
     printf("\n");
     printf("2->Résolution d'un système linéaire par méthode de Gauss-Seidel\n");
     printf("\n");
     printf("3->Calcul du Déterminant d'une matrice\n");
     printf("\n");
     printf("4->Calcul de Matrice Inverse\n");
     printf("\n");
     printf("Q->Quitter le Programme\n");
     scanf("\n %c",&choix);
     switch(choix)
       {
       case '1': systeme();
      break;
       case '2': gauss_seidel();
  break;
       case '3':determinant();
  break;
       case '4':inverse();
  break;
    }
   }while((choix !='q') && (choix!='Q'));
}

A voir également

Ajouter un commentaire Commentaires
Messages postés
1
Date d'inscription
mardi 7 mars 2006
Statut
Membre
Dernière intervention
5 mai 2006

Coucou "bonac" thank u for yr code , sympa .
mais dis moi tu aurais pas dû metre des "fabs" au lieux de "abs" car tu calcul des valeurs absolues de réels ?
et ya un p'tit pb .
-----------------------------
g fais :
cc -Wall -o exe tonfichier.c -lm
./exe
-----------------------------
il mafiche :

1->Resolution d'un systeme lineaire par methode du pivot de Gauss

2->Resolution d'un systeme lineaire par methode de Gauss-Seidel

3->Calcul du Determinant d'une matrice

4->Calcul de Matrice Inverse

Q->Quitter le Programme
-----------------------------
je fais : 4
-----------------------------
j'introduit la matrice :
5.000000 15.000000
7.000000 21.000000
---------------------------
il me dit :
3152519739159347.000000 -2251799813685248.000000
-1050839913053115.625000 750599937895082.625000

???????????????????????????????????????

tu devrai changer ta fonction " metzero " en metant une variable qui compte le nombre de zero dans une ligne (il ne faut pas qu'il existe une ligne plaine de zeros) puis tu fais une sortie de fonction " int metzero " au lieu de " void metzero"
------------------------------------
moi j'ais fais :
/*______________________________________________________*/
/* Applique le pivot situe en A[k][k]*/
int metzero(int dim,double A[][valmax], double b[],int k)
{
int i, j , verifzero , sortiedefonction = 1;
/* r est le rapport A[i][k]/pivot qui sera fixe pour chaque ligne*/
double r, rtemp, pivot=A[k][k];
/*on effectue le calcul pour les lignes et colonnes inferieures a l'element pivot*/
for(i=k+1 ; i<dim ; i++)
{
verifzero = 0 ;
r=A[i][k]/pivot; /*calcul du rapport pour la ligne actuelle*/
for(j=k; j<dim; j++)
{
rtemp = (A[i][j] / A[k][j]) ;
A[i][j]=A[i][j]-(A[k][j]*r); /*on effectue Ligne[i]-r*Ligne[k]*/
if (rtemp == r) verifzero++ ;
}
if (verifzero == (dim-k))
{
sortiedefonction = 0 ;
break ;
}
/*De meme pour les solutions afin de garder l'equivalence entre le nouveau systeme et l'ancien*/
b[i]=b[i]-(b[k]*r);
}
return sortiedefonction ;
}
/*______________________________________________________*/

-----------------------------
-----------------------------

en tout cas un grand merci pour ton programme , il ma beaucoup aidé , etant débutant en C , j'ais apris des truc avec ton programme

un re merci
Messages postés
1
Date d'inscription
samedi 25 septembre 2004
Statut
Membre
Dernière intervention
11 mars 2005

Merci, sans pointeurs c'est une source assez comprehensible pour un neophite de mon espèce,
ce soir je vais tester si ça marche :)
Messages postés
921
Date d'inscription
vendredi 20 décembre 2002
Statut
Membre
Dernière intervention
23 septembre 2010

Dis tu voudrais pas plutôt mettre ça dans un fichier C ? Je veux pas être exigeant ni quoi que ce soit, mais là...
Il y a aussi les accents à arranger, mais à part ça, je trouve parfait...
8/10
(J'ai juste trouvé que le programme était un peu 'froid', car il n'y a aucune instruction claire sur la saisie ou les calculs opérés)

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.