Skip to content
GitLab
Commit d17d9ab7 authored by E-Girl-Mira's avatar E-Girl-Mira
Browse files

job done

parent c69445d0
Hide whitespace changes
Inline Side-by-side
Rapport/Rapport_Gauss_Seidel_Anghel_Pierre.tex 0 → 100644
View file @ d17d9ab7
\documentclass[a4paper, 12pt]{report}
%\documentclass[a4paper, 12pt,draft]{report}
% Cette option permet de compiler plus rapidement car elle n'insère que des cadres au lieu des images
\usepackage[utf8]{inputenc} %package pour le français sous ubuntu : à vous d'adapter
\usepackage[french]{babel} %pour le français
\usepackage[T1]{fontenc} %pour les polices
\usepackage{amsmath} %pour des maths
\usepackage{amsfonts} %pour des maths
\usepackage{amssymb} %pour des maths
\usepackage{graphicx} %pour les inclusions de graphique
\usepackage{palatino} %choix personnel de police
\usepackage{lscape} %si vous voulez des images en paysage
\usepackage{esint} %pour les triples intégrales
\usepackage[left=2cm, right=2cm, top=2cm, bottom=2cm]{geometry} % choix personnel de marges
\usepackage{hyperref} %pour les liens croisés à travers le fichier .pdf
\usepackage{fancyhdr} %pour les marges
\usepackage{float, caption} %pour le positionnement et légende des images
\usepackage{titlesec} %pour redéfinir le titre des sections
\usepackage{color} %pour la couleur dans le code
\usepackage{listings} %pour mettre du code, plutôt en Annexe
\usepackage{lastpage} %pour les références aux pages
\usepackage{epic,eepic} %pour le positionnement de l'image de garde
\usepackage{wrapfig} %pour les images sur le coté dans le texte
\usepackage{calc,ifthen,xspace} %pour redéfinir les espaces et distances
\usepackage{amsmath}
\pagestyle{fancy}
\hypersetup{
dvips,
backref=true, %permet d'ajouter des liens dans...
pagebackref=true,%...les bibliographies
hyperindex=true, %ajoute des liens dans les index.
colorlinks=true, %colorise les liens
breaklinks=true, %permet le retour à la ligne dans les liens trop longs
urlcolor= blue, %couleur des hyperliens
linkcolor= black, %couleur des liens internes
bookmarks=blue, %créé des signets pour Acrobat
bookmarksopen=true}
\definecolor{hellgelb}{rgb}{1,1,0.8} % couleur pour le code
\definecolor{colKeys}{rgb}{0,0,1}
\definecolor{colIdentifier}{rgb}{0,0,0}
\definecolor{colComments}{rgb}{0,0.5,0}
\definecolor{colString}{rgb}{0.62,0.12,0.94}
\definecolor{INSA_GM}{cmyk}{0.6,0,0,0} % et la page de garde
\definecolor{INSA_GRIS}{cmyk}{0.7,0.6,0.5,0.3}
\definecolor{INSA_BLEU}{cmyk}{1,0.9,0.1,0}
%\titleformat{\chapter}[display]%
% {}
% {}
% {0pt}
% {\Huge\bfseries \thechapter. \quad}
% {}
\newcommand{\insertrefproj}[1]{}
\newcommand{\refproj}[1]{\renewcommand{\insertrefproj}{\textbf{\color{INSA_GRIS}#1}}}
\fancypagestyle{courant}{%
\fancyhf{}
\setlength{\headheight}{27pt}
\fancyhead[L]{\raisebox{-2mm}{\includegraphics[width=30mm]{images/insalogo2.png}}}
\fancyhead[C]{}
\fancyhead[R]{\color{INSA_GRIS}\thepage}%\scshape\leftmark}
\fancyfoot[L]{\insertrefproj}
\fancyfoot[R]{}
%\fancyfoot[CE,CO]{\color{INSA_GRIS}\thepage}
\renewcommand{\headrulewidth}{0pt}
\renewcommand{\footrulewidth}{0.2pt}
}
\fancypagestyle{special}{%
\pagestyle{courant}
\fancyfoot{}
%\addtolength{\footskip}{-20pt}
% \setlength{\footskip}{0pt}
% \fancyfoot[C]{}
\renewcommand{\footrulewidth}{0pt}
}
% Redefine the plain page style (for chapters and toc)
\fancypagestyle{plain}{%
\fancyhf{}%
\pagestyle{courant}
}
\newcommand{\Scilab}{
\lstset{
language=Scilab,
float=hbp,
basicstyle=\ttfamily\small,
identifierstyle=\color{colIdentifier},
keywordstyle=\bf \color{colKeys},
stringstyle=\color{colString},
commentstyle=\color{colComments},
columns=flexible,
tabsize=5,
frame=single,
%frame=shadowbox,
rulesepcolor=\color[gray]{0.5},
extendedchars=true,
showspaces=false,
showstringspaces=false,
numbers=left,
stepnumber=5,
firstnumber=1,
numberstyle=\tiny,
breaklines=true,
%backgroundcolor=\color{hellgelb},
captionpos=b,
}
}
\title{Projet MMSN}
\refproj{MMSN V2 - 7} % Entrer ici la reference du projet
\author{ANGHEL Miruna}
\date{}
\begin{document}
\pagenumbering{Alph}
\thispagestyle{empty}
\vspace{4cm}
\begin{picture}(0,0)%(hauteur,largeur)
\put(10,-400){\includegraphics[scale=0.75]{images/matrices.png}}
\put(0,-20){\includegraphics[width=0.4\textwidth]{images/insalogo2.png}}
\put(220,0){\color{INSA_GM}{\begin{minipage}{12cm}\centering \Large \textbf{Projet d'Analyse Numérique}\\\textbf{GM3 - Vague 2}\end{minipage}}}
\put(-10,-200){\color{black}{\begin{minipage}{\textwidth}\centering \Huge \textbf{Résolution de Système Linéaire par la Méthode de Gauss-Seidel}\end{minipage}}}
\newsavebox{\noms}
\savebox{\noms}(250, 300)[tl]{
\put(0,98){\color{INSA_GRIS}{\begin{minipage}{6cm} \textbf{Étudiants :} \end{minipage}}}
\put(0,60){\color{INSA_GRIS}{\begin{minipage}{6cm} Miruna ANGHEL \end{minipage}}}
\put(185,60){\color{INSA_GRIS}{\begin{minipage}{9cm} Edwin PIERRE \end{minipage}}}
\put(0,){\color{INSA_GRIS}\begin{minipage}{9cm} \textbf{Enseignant-responsable du projet :}\\ Jean-Guy CAPUTO
\end{minipage}}
}
\put(100,-850){\usebox{\noms}}
\end{picture}
%\raggedleft
%%%%%%%% Fin de la page de Présentation %%%%%%
\newpage
\pagenumbering{arabic}
\setcounter{page}{1}
\thispagestyle{empty}
\null %rempli la page de rien
\newpage
\pagestyle{special}
\textbf{Date de remise du rapport :} 1/12/2019\\
\textbf{Référence du projet:} GM3 MMSN - Vague 2 - Projet 7 \\
\textbf{Intitulé du projet :} Résolution de Système Linéaire par la Méthode de Gauss-Seidel \\
\textbf{Type de projet :} Programmation en Frotran et Matlab \\
\textbf{Objectifs du projet :} \\
Notre objectif est d'expliciter la méthode de Gauss-Seidel, ses conditions de convergence, l'implémenter en Fortran et en Matlab et tester sur deux cas précis.\\
\vfill
\vfill
\begin{center}
\color{INSA_BLEU}\scshape institut national des sciences appliquées de rouen \\
Département Génie Mathématique \\
685 Avenue de l'Université BP 08- 76801 Saint-Etienne-du-Rouvray \\
\end{center}
\newpage
\pagestyle{courant}
\setcounter{tocdepth}{2}
\tableofcontents
\newpage
\chapter*{Introduction - Description de la Méthode de Gauss-Seidel et ses conditions de convergence}
\addcontentsline{toc}{chapter}{Introduction}
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$. Cette méthode génère donc une suite de vecteurs x qui converge vers une solution de cette équation. L'algorithme suppose que la diagonale de A est formée d'éléments non nuls, car ils sont utilisés comme pivots, donc des divisions par ces termes doivent être effectuées. La convergence de la méthode va donc dépendre du spectre de la matrice $(L+D)^{-1}U$ soit l'ensemble de ses valeurs propres. Les conditions de convergence de cette méthode pour tout vecteur résultat b et tout vecteur initial x sont donc très simples, mais non triviales à avoir :\\
\begin{itemize}
\item La matrice A doit être symétrique définie positive
\item La matrice A doit être à diagonale strictement dominante
\end{itemize}
\vspace{\baselineskip}
La méthode de Gauss-Seidel se base sur une décomposition de la forme $A = L + D + U$, ou L, D et U sont des matrices carrées de la même dimension que la matrice A, L est la partie triangulaire inferieure de la matrice A, D est sa diagonale et U est sa partie triangulaire supérieure. L'algorithme de Gauss-Seidel comporte une descente qui peut être vue comme une méthode de splitting, selon de cours d'Analyse Numérique. On a donc à résoudre l’équation : \\
$(L + D)x^{k+1} = b - Ux^k$\\
En somme, nous cherchons le point fixe de $x = (L + D)^{-1}(b - Ux)$ par approximations successives.\\
\chapter{Implémentation de la méthode de Gauss-Seidel en Fortran}
\begin{lstlisting}
program main
use matop
implicit none
integer i, j, n, nmax
real, allocatable, dimension(:,:) :: A,F,D1,L,LI,U,Eye,tmp
real, allocatable, dimension(:) :: b, x, G
real :: prec, res
read (*, *) n, prec, nmax
allocate(A(n,n))
allocate(F(n,n))
allocate(D1(n,n))
allocate(L(n,n))
allocate(LI(n,n))
allocate(U(n,n))
allocate(Eye(n,n))
allocate(tmp(n,n))
allocate(b(n))
allocate(x(n))
allocate(G(n))
read (*, *) A
read (*, *) b
read (*, *) x
D1 = 0
Eye = 0
do i = 1, n
D1(i, i) = 1.0 / A(i, i)
Eye(i, i) = 1
enddo
tmp = matmul(D1, A)
call tril(tmp, L)
call triu(tmp, U)
U = U - Eye
call inv_trig_inf(L, LI)
F = -matmul(LI, U)
G = matmul(matmul(LI, D1), b)
deallocate(D1)
deallocate(L)
deallocate(LI)
deallocate(U)
deallocate(Eye)
deallocate(tmp)
do i = 1, nmax
x = matmul(F, x) + G
res = norm2(matmul(A, x) - b)
write(*, *) i, res
if (res < prec) then
exit
endif
enddo
write (*,*) x, matmul(A, x)
deallocate(F)
deallocate(G)
deallocate(x)
deallocate(A)
deallocate(b)
end
\end{lstlisting}
Pour simplifier l'écriture et augmenter la modularité du code nous avons créé le module ``matop'' = Matrix Operations qui contient nos algorithmes de création de la matrice triangulaire inférieure L et supérieure U comme vu en TD, ainsi que l'algorithme d'inversion de la matrice L. Ca nous permet aussi de contourner une des limitations du Fortran et ainsi utiliser size(mat).
\pagebreak
\begin{lstlisting}
module matop
implicit none
contains
subroutine tril(M, L)
implicit none
integer :: n, i
real, dimension(:,:) :: M, L
n = size(M, 1)
L = 0
do i=1,n
L(i,1:i) = M(i,1:i)
end do
end subroutine tril
subroutine triu(M, U)
implicit none
integer :: n, p, i
real :: M(:,:), U(:,:)
n = size(M, 1)
p = size(M, 2)
U = 0
do i=1,n
U(i,i:p) = M(i,i:p)
end do
end subroutine triu
subroutine inv_trig_inf(L, inv)
implicit none
integer :: i, j, n
real, allocatable, dimension(:,:) :: L(:,:), inv
real, dimension(size(L,1),size(L,2)) :: Eye
n = size(L, 1)
inv = 0;
Eye = 0
do i=1,n
Eye(i,i) = 1
enddo
do i=1,n
inv(i,:) =(Eye(i,:)-matmul(L(i,1:i-1),inv(1:i-1,:)))/L(i,i)
enddo
return
end subroutine inv_trig_inf
end module matop
\end{lstlisting}
\chapter{Implémentation de la méthode de Gauss-Seidel en Matlab}
Matlab étant un langage prévu pour le calcul matriciel, l'implémentation a été beaucoup plus simple.\\
\begin{lstlisting}
A = [1 2 3 4 5 ; 5 1 2 3 4 ; 4 5 1 2 3 ; 3 4 5 1 2 ; 2 3 4 5 1];
Abis = [1 -1 -200 0 ; 0 1 100 -100 ; 0 1 101 -101 ; 0 0 0 100];
b = [1 ; 1 ; 3 ; 4 ; 2];
bbis = [1 ; 1 ; 2 ; 4];
x = [0 ; 0 ; 0 ; 0 ; 0];
xbis = [0 ; 0 ; 0 ; 0];
tol = 1e-10;
y = x+2*tol;
n = length(b);
cont = 0;
tic;
while (cont < 200) && ((norm(x - y) > tol))
y = x;
for i=1:n
x(i) = (b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:n)*x(i+1:n)) / A(i,i);
endfor
cont = cont + 1;
endwhile
tempsCalcul = toc();
fprintf('%f\n%f\n%f\n%f\n en %d boucles',x, cont);
\end{lstlisting}
\chapter{Résultats}
\section{Résultats théoriques}
Pour la premiere équation:\\
\begin{figure}[ht]
\begin{center}
\advance\leftskip-3cm
\advance\rightskip-3cm
\includegraphics[keepaspectratio=true,scale=0.75]{images/A.png}
\end{center}\end{figure}
La matrice A n'est ni symétrique définie positive ni à diagonale strictement dominante, donc les conditions de convergence de la méthode de Gauss-Seidel ne sont pas remplies. La méthode ne donne donc aucune garantie de convergence. Ca va dépendre du spectre de $(L + D)^{-1}U$ soit:\\
$$
\begin{pmatrix}
0 & 2 & 3 & 4 & 5 \\
0 & -10 & -13 & -17 & -21 \\
0 & 42 & 53 & 71 & 88 \\
0 & -176 & -222 & -299 & -369 \\
0 & 738 & 931 & 1254 & 1546
\end{pmatrix}
$$
Dont le rayon spectral est d'environ 1289 et dont nous ne chercherons pas à regarder le spectre. La méthode de Gauss-Seidel va certainement diverger pour cette matrice.
\pagebreak
Pour la deuxieme equation:\\
\begin{figure}[ht]
\begin{center}
\advance\leftskip-3cm
\advance\rightskip-3cm
\includegraphics[keepaspectratio=true,scale=0.75]{images/Aprim.png}
\end{center}\end{figure}
Cette fois encore la matrice n'est pas symétrique définie positive ni à diagonale strictement dominante, donc on ne peut pas prévoir si l'algorithme convergera ou non. On a donc la matrice $(L + D)^{-1}U$ :\\
$$
\begin{pmatrix}
0 & -1 & -200 & 0 \\
0 & 0 & 100 & -100 \\
0 & 0 & -0.99010 & -0.00990 \\
0 & 0 & 0 & 0
\end{pmatrix}
$$
Cette fois le rayon spectral est de 0.99010, soit plus petit que 1, donc il y a des chances pour que la méthode de Gauss-Seidel converge, quoi que lentement, selon les choix de b et x.
\section{Résultats pratiques}
Toutes nos simulations sont faites avec les paramètres suivants : \\
\begin{itemize}
\item n = 5 pour la première matrice et 4 pour la deuxième
\item tol = 1.e-10
\item nmax = 200
\item Le vecteur x initial est toujours le vecteur nul
\end{itemize}
\vspace{\baselineskip}
\subsection{Fortran}
Pour la première matrice les résultats sont :\\
1 2048.69263\\
2 2643929.75\\
3 3.40974874E+09\\
4 4.39738905E+12\\
5 5.67110112E+15\\
6 7.31374574E+18\\
7 9.43218763E+21\\
8 1.21642384E+25\\
9 1.56876318E+28\\
10 2.02315839E+31\\
11 2.60917007E+34\\
12 3.36492122E+37\\
A partir de la 13e itération le résultat est NaN, ce qui confirme nos attentes quant à la divergence de la méthode de Gauss-Seidel pour cette matrice.\\
Pour la deuxième matrice, au bout de 200 itération on obtient : \\
$$
\begin{pmatrix}
95.7803574 \\
-84.6367340 \\
0.897789419 \\
3.99999991E-02
\end{pmatrix}
$$
Et le résidu est en décroissance constante. La méthode est donc entrain de converger, mais très lentement. Nous avons donc essayé d'augmenter le nombre d'itérations pour voir si on pouvait trouver un résultat satisfaisant par rapport à notre tolérance mais au vu d'un résidu de 1.17437794E-05 au bout de la 20000e itération, nous en sommes restés la.\\
\subsection{Matlab}
Tout comme pour le code Fortran, la première matrice a divergé. Le résultat renvoyé par Octave suite à l'exécution du code est :\\
\begin{figure}[ht]
\begin{center}
\advance\leftskip-3cm
\advance\rightskip-3cm
\includegraphics[keepaspectratio=true,scale=0.75]{images/matlabA.png}
\end{center}\end{figure}
La divergence a lieu à partir de la 13e itération, tout comme pour le code Fortran.
Quant a la seconde matrice :\\
\begin{figure}[ht]
\begin{center}
\advance\leftskip-3cm
\advance\rightskip-3cm
\includegraphics[keepaspectratio=true,scale=0.75]{images/matlabAprim.png}
\end{center}\end{figure}
Les résultats sont donc identiques entre les deux langages.
\pagebreak
\addcontentsline{toc}{chapter}{Conclusion}
\chapter*{Conclusion}
En vue de nos résultats expérimentaux, nous pouvons conclure qu'ils correspondent à nos attentes théoriques, conformément au cours d'Analyse Numérique ainsi qu'aux simulations déjà vues en cours de Calcul Scientifique sous Unix dans le cas du code Matlab. Par ailleurs, le fait que les résultats soient identiques entre les deux langages augmente la fiabilité des résultats ainsi que celle des codes employés. \\
La méthode de Gauss-Seidel est donc très utile et a un bon temps de calcul dans ses deux cas de convergence (matrice symétrique définie positive ou à diagonale strictement dominante) mais devrait être évitée si au moins une de ces deux conditions n'est pas remplie car elle n'offre aucune garantie et ne peut pas être adaptée pour des cas de matrices aléatoires.\\
La méthode de Gauss-Seidel est obsolette en termes de résolution d'équations linéaires mais reste une bonne introduction aux méthodes itératives.
\end{document}
gauss_seidel.m
View file @ d17d9ab7
%A = [1 2 3 4 5 ; 5 1 2 3 4 ; 4 5 1 2 3 ; 3 4 5 1 2 ; 2 3 4 5 1];
%Abis = [1 -1 -200 0 ; 0 1 100 -100 ; 0 1 101 -101 ; 0 0 0 100];
%b = [1 ; 1 ; 3 ; 4 ; 2];
%bbis = [1 ; 1 ; 2 ; 4];
%x = [0 ; 0 ; 0 ; 0 ; 0];
%xbis = [0 ; 0 ; 0 ; 0];
A = [1 2 3 4 5 ; 5 1 2 3 4 ; 4 5 1 2 3 ; 3 4 5 1 2 ; 2 3 4 5 1];
Abis = [1 -1 -200 0 ; 0 1 100 -100 ; 0 1 101 -101 ; 0 0 0 100];
b = [1 ; 1 ; 3 ; 4 ; 2];
bbis = [1 ; 1 ; 2 ; 4];
x = [0 ; 0 ; 0 ; 0 ; 0];
xbis = [0 ; 0 ; 0 ; 0];
tol = 1e-10;
y = xbis+2*tol;
n = length(bbis);
cont = 0;
tol = 1e-10;
cont = 0;
delta = 10;
A = [3 1;2 4];
b = [3;3];
x = [1;1];
n = size(b);
while delta > tol & cont < 200
ref = x;
for i = 1:n(1)
tmp = 0;
for j = 1:n(1)
if j != i
tmp = tmp + A(i,j)*x(j);
end
x(i) = (b(i) - tmp) / A(i,i);
end
end
delta = norm(ref - x);
cont = cont + 1;
end
fprintf('%f\n%f\n%f\n%f\n en %d boucles',x, cont);
\ No newline at end of file
tic;
while (cont < 200) && ((norm(xbis - y) > tol))
y = xbis;
for i=1:n
xbis(i) = (bbis(i) - Abis(i,1:i-1)*xbis(1:i-1) - Abis(i,i+1:n)*xbis(i+1:n)) / Abis(i,i);
endfor
cont = cont + 1;
endwhile
tempsCalcul = toc();
fprintf('%f\n%f\n%f\n%f\n en %d boucles',xbis, cont);
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment