RÉSOLUTION DE SYSTÈMES D'ÉQUATIONS

 
#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'));
}