Résolution de systèmes d'équations
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'));
}
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.