Projet MatLab

Université Sidi Mohamed Ben Abdallah Faculté des Sciences Dhar El Mahraz Fès Département d ’Informatique

Résolution des Systèmes Linéaires par la méthode du Pivot total de Gauss Sous Matlab

Réalisé par : WARME Moumouni Etudiant en Master I

Filière : Master Qualité Logiciel Année Universitaire 2010-2011

La méthode du Pivot total de Gauss

W.Moumouni

Table des matières Introduction ............................................................................................................................................. 3 A.

Les Systèmes Linéaires .................................................................................................................... 3 1.

Les notions générales .................................................................................................................. 3

2.

Les méthodes directes................................................................................................................. 4

B.

Résolution des systèmes triangulaires .................................................................................... 4

b.

Méthode d'élimination de Gauss et décomposition LU . ........................................................ 5 1.

Description de l'algorithme ................................................................................................. 5

2.

Stratégie de choix du pivot .................................................................................................. 5

Présentation de l’algorithme........................................................................................................... 7 1.

La fonction de choix du pivot et de permutation ........................................................................ 7

2.

La fonction de décomposition ..................................................................................................... 8

3.

La fonction de la méthode de montée ........................................................................................ 9

C.

D.

a.

Présentation du programme sous MatLab...................................................................................... 9 1.

Le code sous Matlab de la fonction de décomposition............................................................... 9

2.

Le code sous Matlab de la fonction de choix du pivot et de permutation ............................... 10

3.

Le code sous Matlab de la fonction de montée plus le programme......................................... 11 Test du le programme en Matlab .................................................................................................. 12

Conclusion ............................................................................................................................................. 14

2


W.Moumouni

Introduction La résolution de systèmes linéaires peut être effectuée par plusieurs méthodes. Ainsi le choix d’une des méthodes rend cette résolution plus ou moins complexe par rapport au choix d’autres méthodes. Sans oublier que la solution choisie peut donnée une solution qui s’approche ou s’écarte de la solution exacte. Les méthodes de cette résolution étant nombreuses, nous ont poussé à étudier une seule méthode en particulier celle du pivot total de Gauss, afin de fournir un algorithme de résolution et implémenter cet algorithme sous Matlab pour en faire une application. Avant de commencer nous allons faire un rappel sur les systèmes linéaires.

A. Les Systèmes Linéaires 1. Les notions générales On appelle système linéaire d'ordre n (n entier positif), une expression de la forme Ax = b; Où A = (a ij ), 1 ≤ i; j ≤ n , désigne une matrice de taille n x n de nombres réels ou complexes, b=(b i) ,1 ≤ i ≤ n , un vecteur colonne réel ou complexe et x = (x i ), 1 ≤ i ≤ n , est le vecteur des inconnues du système. La relation précédente équivaut aux équations

La matrice A est dite régulière (ou non singulière) si detA≠ 0 ; on a existence et unicité de la solution x (pour n'importe quel vecteur b donné) si et seulement si la matrice associée au système linéaire est régulière. Théoriquement, si A est non singulière, la solution est donnée par la formule de Cramer :

où Ai est la matrice obtenue en replaçant la i -ème colonne de A par le vecteur b. Cependant l'application de cette formule est inacceptable pour la résolution pratique des systèmes, car son coût est de l'ordre de (n + 1)! floating-point opérations (flops : opération avec des virgules). En fait, le calcul de chaque déterminant par la formule suivante:

(où la somme est étendue à toutes les permutations sur n objets) requiert n ! flops. 9 47 Par exemple, sur un ordinateur effectuant 10 flops par seconde il faudrait 9,6 10 années pour résoudre un système linéaire de seulement 50 équations. Il faut donc développer des algorithmes alternatifs avec un coût raisonnable.

3


2.

W.Moumouni

Les méthodes directes

On appelle méthode de résolution directe d'un système linéaire un algorithme qui, si l'ordinateur faisait des calculs exacts, donnerait la solution en un nombre fini d'opérations. Il existe aussi des méthodes itératives qui consistent à construire une suite de vecteurs xn convergeant vers la solution x. Mais ces méthodes n’étant l’objet de notre présent projet , nous éviterons de les aborder. Pour résoudre Ax = b, on cherche à écrire A = LU où - L est une matrice triangulaire inférieure avec des 1 sur la diagonale, - U est une matrice triangulaire supérieure. La résolution de Ax = b est alors ramenée aux résolutions successives des systèmes échelonnés Ly = b et U x = y.

a. Résolution des systèmes triangulaires Une matrice A = (a ij ) est triangulaire supérieure quand a ij = 0 pour i>j et triangulaire inférieure si a ij = 0 pour i
Cet algorithme est appelé méthode de descente . Si A est triangulaire supérieure on a :

Cet algorithme est appelé méthode de remontée .

4


W.Moumouni

b. Méthode d'élimination de Gauss et décomposition LU .

1. Description de l'algorithme La méthode consiste à rendre la matrice A soit une matrice triangulaire supérieure. Soit A = (a ij ) une matrice non singulière n xn . L'algorithme de Gauss est le suivant : (1) (1) On pose A = A, c'est-à-dire a ij = a ij pour i, j = 1,…, n ; pour k = 1,…,n on calcule

A la fin de ce procédé, les éléments de la matrice U (triangulaire supérieure) sont définis par

tandis que les éléments de la matrice L (triangulaire inférieure) sont les termes l ij engendrés par l'algorithme de Gauss. (k ) (k +1) (k ) En fait, pour passer de A à A on modifie touts les éléments de A sauf les premières k lignes et les premières k colonnes, qui ne changent pas.

Remarque : (k ) Pour que l'algorithme de Gauss puisse terminer, il faut que tous les termes a kk , qui correspondent (k ) aux termes diagonaux u kk de la matrice U et qu'on appelle pivots, soient non nuls. Parlant de pivot nous devons bien choisir celui-ci pour que les solutions soient exactes. 2. Stratégie de choix du pivot Que se passe-t-il si on adopte la stratégie "simple" consistant à choisir le premier coefficient rencontré dans la colonne considérée, par exemple a1,1 si a1,1 ≠ 0 ? Cela peut donner lieu à de grandes erreurs d'arrondi, comme le montre l'exemple suivant :

On suppose que la machine fonctionne avec 3 chiffres significatifs en virgules flottantes Par exemple 0.0001, pi = 3:14159…., 1/3 = 0.33333… sont convertis respectivement en :

valeurs approchées sauf la première qui se trouve être exacte.

Remarque : Une autre convention souvent utilisée est l'écriture commençante par 0. Comme dans +1 pi ~ 0.314x10 . Dans 1.00 il faut remarquer que les deux 0 sont significatifs et que 1.00 est l'écriture approchée de tout nombre strictement compris entre 0.995 et 1.005. Rappelons que les calculateurs utilisent le signe « . » plutôt que la virgule de la tradition française pour séparer les entiers des -4 décimales. Dans l'exemple (S) et selon "la stratégie simple" on prend pour pivot a1,1 = 1.00x10 . +4 On transforme donc la deuxième ligne par : L2-10 L1. Le calcul exact est

5


W.Moumouni

mais la machine prend la version arrondie :

et la résolution du système triangulaire obtenu est alors : y = 1 puis 0.0001x = 1 - y = 0 donc la solution trouvée est x=0;y=1 ce qui est très éloignée de la solution exacte tirée de (S1) : x = 1/0.9999 = 1.00010001 … ~1.00

Explication : Lorsqu'on calcule x dans la première équation 0:0001x = 1 y = 0, on divise par un petit pivot 0.0001, et l'erreur sur y dans le calcul approche de la deuxième équation se trouve -3 multipliée par 10000 : l'erreur d'arrondi qui peut aller initialement jusqu' à 0.5 ×10 devient du même ordre de grandeur que les solutions. Si au contraire, on échange les lignes et on choisit pour pivot a2,1 = 1, la résolution donne le

Système triangulaire : -3

D'ou y = 1 et cette fois x = 2 - y = 1, ce qui est correct à 10 près.

Conclusion : ième Pivot total, pivot partiel. Même si à l'issue de la (k-1) étape on trouve on vient de voir qu'il est préférable de choisir un autre pivot. (k-1) ième La stratégie du pivot partiel consiste 0 choisir a i,k avec i ≥ k de module maximum dans la k colonne de la matrice A. (k-1) La méthode du pivot total , consiste à utiliser le plus grand des coefficients a i,k, parmi toutes les paires d'indices (i, j), ce qui mène à des échanges de lignes et de colonnes. La méthode est plus précise mais plus coûteuse en temps. De plus, il faut prendre garde au fait que la solution du système transformé n'est pas u, mais un vecteur obtenu en permutant les coordonnées de u, de la même façon qu'on a permute les colonnes. Notre exemple de mauvais choix de pivot a1,1 d’écrit un cas particulièrement défavorable où le rapport a2,1/a1,1 est le plus élevé possible. Pour une matrice assez générale, le choix simple de pivot à chaque étape est souvent numériquement acceptable et donne alors lieu à une factorisation de la matrice A en produit de matrices triangulaires. Après avoir expliqué les différentes méthodes et leur faiblesse nous allons nous intéresser à la méthode du pivot total qui est une méthode plus précise que les autres. Ainsi nous allons donner un algorithme informatique permettant la résolution de systèmes linéaire par cette méthode et en suite implémenter cet algorithme sous Matlab.

6


W.Moumouni

B. Présentation de l’algorithme L’analyse des différentes méthodes permettant de mieux comprendre la méthode du pivot total, puisque celle-ci n’intervient et est plus précise que les autres méthodes, dans les cas où le système initial à des pivots très petits. Rappelons brièvement cette méthode qui consiste à :

1. Choisir le pivot le plus grand dans toute la matrice A 2. Effectuer des permutations de lignes ou de colonnes en fonction de la localisation du meilleur pivot (plus grand). 3. Rendre la matrice A triangulaire supérieure ou inférieure dans notre cas triangulaire supérieur. 4. Appliquer la méthode de la montée ou descente en fonction du type de la matrice A obtenu dans l’étape ci-dessus pour avoir la résolution. 5. Mettre en ordre les solutions en fonction des permutations de colonnes effectuée. Ces différentes étapes nous amènent à penser à scinder la résolution en fonctions qui vont interagir pour résoudre un système donné. Nous allons détailler l’algorithme de ces fonctions . Ensuite nous présenterons les implémentations sous Matlab après.

1. La fonction de choix du pivot et de permutation Voici l’algorithme de cette fonction, qui reçoit en paramètre la matrice carrée A, la matrice colonne b et la matrice colonne pivot_sol qui à pour rôle de donner l’ordre des solutions . Cette fonction renvoie comme résultat les matrices A, b, et pivot_sol après permutation.

fonction [A,b,pivot_sol] :=permut(A,b,pivot_sol) Debut n:=lireNbreLigne(A); (*Lecture du nombre de ligne de la matrice*); m:=lireNbreColonne(A); (*Lecture du nombre de colonne de la matrice*); (*Initialisation de l'ordre des solutions*) Pour i=1 à n faire pivot_sol(i):=i; fait; Pour k=1 à n-1 faire pivo:=0; (*Choix du pivot le plus grand dans la matrice A*) Pour i=k à n faire Pour j=k à n faire Si abs (A(i,j))>pivo faire (*Indexation du numéro de ligne et de colonne du pivot*) pivo:=abs(A(i,j)); io:=i;jo:=j; fsi; fait; fait; Si pivo≠0 alors (*Permutation des lignes de la matrice A si le pivot est différent de zéro et est le plus grand*)

Pour j=K à n faire temp:=A(k,j); A(k,j):=A(io,j); A(io,j):=temp; fait (*Permutation des colonnes pour la matrice b*) temp:=b(k); b(k):=b(io); b(io):=temp; (*Permutation des colonnes de la matrice A si le pivot est différent de zéro et est le plus grand*) Pour i=K à n faire 7


W.Moumouni

temp:=A(i,k); A(i,k):=A(i,jo); A(i,jo):=temp;

fait; fait; (*Actualisation de l'ordre des solutions avec le nouvelle ordre après permutation des colonnes de la matrice A*)

temps:=pivot_sol(k); pivot_sol(k):=pivot_sol(jo); pivot_sol(jo):=temps;

Si A(k,k)==0 alors Afficher ('Un pivot nul ! => méthode de Gauss pivot total non applicable'); fsi; Si pivo#0 alors (*La sortie de la grande boucle pour éviter plusieurs permutations incohérent avec le principe du pivot total*)

fin boucle;

fsi; fait; Fin

2. La fonction de décomposition Elle retourne comme résultat la matrice A en une matrice triangulaire supérieure L’algorithme de cette fonction est la suivante :

fonction [A,b] :=triang(A,b) Debut n:=lireNbreLigne(A);(*Lecture du nombre de ligne de la matrice*); m:=lireNbreColonne(A);(*Lecture du nombre de colonne de la matrice*); Pour k=1 à n-1 faire pivot:=A(k,k); Si pivot≠0 alors Pour i=k+1 à n faire

b(i)=b(i)-(A(i,k)/pivot)*b(k); Pour j=k+1 à n faire

A(i,j):=A(i,j)-(A(i,k)/pivot)*A(k,j); fait fait Sinon Afficher('pivot null'); fait fait (*Mettre à Zero les elements qui doivent etre des zéro*)

Pour i=1 à n faire Pour j=1 à n faire Si(i>j) alors A(i,j) :=0; fait fait fait Fin

8


W.Moumouni

3. La fonction de la méthode de montée Cette fonction est intégrée dans la fonction globale de résolution.

fonction x :=montee(A,b) Début (*Calcul de xn*)

x(pivot_sol(n)):=b(n)/A(n,n); Pour i=n-1 à :1 s:=0; Pour j=i+1 à n faire

s:=(s+A(i,j)*(x(pivot_sol(j)))); fait (*Calcul des xi*)

x(pivot_sol(i)):=((b(i)-s)/A(i,i)); fait Fin

Toute ces fonction interagir dans une fonction globale pour permettre la résolution de système linéaire.

C. Présentation du programme sous MatLab Dans Matlab les fonctions sont dans des fichiers nous allons donc présenter le contenu de ces fichiers.

1. Le code sous Matlab de la fonction de décomposition function [A,b]=triang(A,b) [n,m]=size(A); for k=1:n-1 pivot=A(k,k); if pivot~=0 for i=k+1:n b(i)=b(i)-(A(i,k)/pivot)*b(k); for j=k+1:n A(i,j)=A(i,j)(A(i,k)/pivot)*A(k,j); end end elseif disp('pivot null'); end end % // Mettre à Zero les elements qui doivent etre des zéro for i=1:n for j=1:n if(i>j)A(i,j)=0; end end end end

9


W.Moumouni

2. Le code sous Matlab de la fonction de choix du pivot et de permutation function [A,b,pivot_sol]=permut(A,b,pivot_sol) [n,m]=size(A); %Initialisation de l'ordre des solutions for i=1:n pivot_sol(i)=i; end for k=1:n-1 pivo=0; %Choix du pivot le plus grand dans la matrice A for i=k:n for j=k:n if(abs(A(i,j))>pivo) pivo=abs(A(i,j)); io=i; jo=j; end end end if pivo~=0 %Permutation des lignes de la matrice A si le pivot est différent %de zéro et est le plus grand temp=A(k,:); A(k,:)=A(io,:); A(io,:)=temp; %Permutation des colonnes pour la matrice b temp=b(k); b(k)=b(io); b(io)=temp; %Permutation des colonnes de la matrice A si le pivot %est différent de zéro et est le plus grand temp=A(:,k); A(:,k)=A(:,jo) ; A(:,jo)=temp; end %Actualisation de l'ordre des solution avec le nouvelle ordre après %permutation des colonnes de la matrice A temps=pivot_sol(k); pivot_sol(k)=pivot_sol(jo); pivot_sol(jo)=temps; if A(k,k)==0 disp('Un pivot nul ! => methode de Gauss pivot total non applicable'); end if pivo~=0 %La sortie de la grande boucle pour éviter plusieurs %permutations incohérent avec le principe du pivot total break; end end

10


W.Moumouni

3. Le code sous Matlab de la fonction de montée plus le programme Dans ce code il y a l’interaction des différent es fonctions définies plus haut.

function y=gauss_total(A,b) [l,c]=size(b); [n,m]=size(A); if n==m if c==n pivot_sol=0; %Choix du pivot maximun et Permutation des lignes et colonnes

[A,b,pivot_sol]=permut(A,b,pivot_sol); %triangularisation

[A,b]=triang(A,b); %Presentation de la matrice A triangulaire supérieure disp('Matrice A triangulaire supérieure'); A disp(' '); %Résolution par la methode de la montée %Car la matrice A est une matrice triangulaire supérieur

x(pivot_sol(n))=b(n)/A(n,n); for i=n-1:-1:1 s=0; for j=i+1:n s=(s+A(i,j)*(x(pivot_sol(j)))); end x(pivot_sol(i))=((b(i)-s)/A(i,i)); end y=x; for i=1:n sol_ordre(pivot_sol(i))=y(i); end disp('Solution en ordre indication de l''ordre par sol_ordre'); for i=1:n disp(['X' num2str(i) '=' num2str(sol_ordre(i),'%10.4f')]); end disp(' '); disp('Solution en désordre indication de l''ordre par pivot_sol'); for i=1:n disp(['X' num2str(pivot_sol(i)) '=' num2str(x(i),'%10.4f')]); end elseif c~=n disp('=>>>>>La matrice b n ''est pas correcte '); disp('=>>>>>Veuillez saisir une matrice b avec le nombre de colonne égal au nombre de ligne de la matrice A '); end elseif n~=m disp('=>>>>>La matrice n est pas carrée '); disp('=>>>>>Veuillez saisir une matrice carrée '); end end

11


W.Moumouni

D.Test du le programme en Matlab Après avoir écrit le programme, nous allons effectuer un test qui consiste à entrer un système d’équation Ax=b et relever le résultat donné par le programme afin de faire une comparaison avec la solution exacte. Test 1 : Par exemple pour le système d'équations Ax=b suivant: 1

2

4

8

16

31

2

4

8

16

31

61

8 16

31

61

16 31 61

120

A= 4 8

16 31 61 120

236

b=

120 236 464

La solution exacte est x = ( 1 , 1 , 1 , 1 , 1) Voici les matrice que nous allons saisir dans Matlab : A=[1 2 4 8 16;2 4 8 16 31;4 8 16 31 61;8 16 31 61 120 ;16 31 61 120 236 ] b=[31 61 120 236 464] Ensuite nous allons appelé la fonction gauss_total(A,b).

12


W.Moumouni

L’appelle de la fonction gauss_total(A,b) et affichage des solutions :

A travers le résultat on constate bien qu’il y a eu une permutation de colonnes et de ligne, ceci est montré par la matrice A triangulaire supérieure .Car la dernière ligne a été permutée avec la première ligne et la dernière colonne avec la première. Nous constatons aussi que le résultat est exacte quant à l’ordre des solutions il n’y a pas grand-chose puis que les résultats sont tous égaux à 1. Pour bien voire que le programme remet les résultats en ordre en fonction des permutations de colonnes nous allons prendre un autre exemple. Test 2:

Soit le système d’équation suivant Ax=b comme suit : 4

8

12

4

A= 3

8

13

b= 5

2

9 18

11

13


W.Moumouni

Dans cette capture qui suit on voit très bien que le programme garde l’ordre des solutions même avec les permutations de colonnes.

Conclusion Nous sommes maintenant arrivés au terme de notre projet dont nous allons donner un bref bilan par rapport aux objectifs que visaient celui-ci. Ce projet avait pour but de réaliser la méthode du pivot total de Gauss pour la résolution des systèmes linéaires. Ce qui nous a permis d’une part de bien maitriser les algorithmes de cette méthode de résolution, et d’autre part une familiarisation avec l’outil de programmation mathématique MatLab version 7.6.0 de l’année 2008. Sans oublier qu’il serait très intéressant, de donner un programme qui rassemble tout les méthodes de résolution de systèmes linéaires, ce qui nous permettrait de comparer les différentes solutions et de voir exactement les erreurs introduites par certaines méthodes.

14

Projet MatLab

Recommend Documents