1. Méthodes Itératives Stationnaires Standards

Objectifs

1. Implémenter toutes les méthodes itératives standards
2. Comparer la vitesse de convergence des différentes méthodes

Algorithme générique

Pour résoudre le problème linéaire $A x= b$, une technique consiste à réécrire la matrice $A$ sous la forme $A = M - (M-A)$, où $M$ est une matrice bien choisie, et d’appliquer le schéma itératif suivant \begin{equation}\label{eq:split} Mx_{n+1} = (M-A)x_n + b \end{equation} à convergence, nous retrouvons bien $(M-(M-A))x = Ax = b$.

L’algorithme de résolution décrit par l’équation \eqref{eq:split} s’écrit alors, pour une tolérance $\varepsilon$ et un nombre maximal d’itérations $n_{max}$ donnés :

niter = 0           // nombre d'itérations
x = 0               // vecteur solution
r = b               // vecteur résidu
resvec[0] = |r|/|b| // Résidu relatif de l'itération 0
while( (|r|/|b| > tolerance) && (niter < n_max))
  y = M^{-1}*r            // Calcul de la direction (résolution système)
  x = x + y               // Mise à jour de la solution courante
  r = b - Ax              // Mise à jour du résidu (vecteur)
  niter = niter + 1       // Mise à jour du numéro de l'itération
  resvec[niter] = |r|/|b| // Stockage du résidu relatif (à des fins d'analyse)
end

Les méthodes itératives que nous détaillons ici ne diffèrent que par le choix de la matrice $M$. Pour chacune d’entres elles, nous construirons une classe dédiée, bien qu’elles partagent toutes des codes similaires (nous pourrions utiliser l’héritage, mais cela compliquera la “templatisation”).

Méthodes standards

La matrice $M$ est choisie selon le splitting (découpe) d’une matrice $A$ comme $A = D - E - F$, où $D$, $E$ et $F$ sont des matrices de la même taille que $A$ et telles que (attention au signe des coefficents de $E$ et $F$) :

• $D = {\rm diag}(A)$ : Matrice ne contenant que les termes diagonaux de $A$
• $E$ : Matrice ne contenant que la partie triangulaire inférieure stricte de $-A$
• $F$ : Matrice ne contenant que la partie triangulaire supérieure stricte de $-A$

Ensuite, selon les choix pour la matrice $M$, la méthode obtenue sera différente :

Une Classe en Détail : Jacobi

Pour chaque méthode, nous proposons d’implémenter une classe distincte et donc, de créer deux fichiers (un en-tête et un source). Par exemple, pour Jacobi : include/jacobi.hpp et src/jacobi.cpp. Les méthodes se ressemblant, ces différentes classes se ressembleront naturellement. Nous détaillons ici la classe Jacobi, cependant, pour Gauss-Seidel et la Relaxation, la classe sera similaire à quelques modifications près (par exemple la Relaxation nécessite un paramètre $\omega$ en plus).

Données Membres (ou paramètres)

Celles-ci seront séparées en deux, les données “entrantes”, fournies par l’utilisateur, et les données “sortantes”, calculées lors de la résolution du problème linéaire. En pratique, rien ne différencie ces deux types de données qui sont de type private, seule leur utilisation permet de les distinguer :

1. Les données entrantes sont fournies par l’utilisateur et ne sont pas modifiées par la classe
2. Les données sortantes sont modifiées par la classe (vecteur solution, nombre d’itération, …) et seront récupérées par l’utilisateur à des fins d’analyses

Nous proposons les paramètres suivants, mais vous pouvez bien entendu en ajouter !

La figure suivante illustre ces paramètres et la classe Jacobi :

Constructeurs

À vous de décider ce dont vous avez besoin : un constructeur vide ? Un qui prend toutes les données d’entrée en argument, par exemple :

Jacobi(const Matrice &A, const Vecteur &b, double tol, int maxit);

Jacobi(const Matrice & A, const Vecteur & b,...)
  :A_(A), b_(b)... // Liste d'initialisation
{
    blabla
}

Méthode Solve()

Outre les accesseurs (getter) et les mutateurs (setter) habituels pour respectivement accéder et modifier les paramètres, nous avons au moins besoin d’une fonction membre qui résout le système linéaire en appliquant l’algorithme de résolution. Celle-ci aura (probablement) le prototype suivant :

void Jacobi::Solve();

D’autre part, voici quelques propositions pour définir la méthode :

• La solution sera stockée dans x_
• Le nombre d’itérations sera stocké dans niter_
• À chaque itération, la norme du résidu relatif $\frac{\|r\|}{\|b\|}$ doit être calculé pour vérifier si la convergence est atteinte ou non. Cette quantité sera stockée dans le tableau resvec_ au fur et à mesure des itérations. Ceci nous sera utile pour la partie analyse.

Implémentation et Test

Pour tester et valider votre code, utiliser la même matrice que pour la factorisation LU : $$ \begin{pmatrix} 2 & -1 & 0 & 0 &0\\\
-1 & 2 & -1 & 0 &0\\\
0 & -1 & 2 & -1 &0\\\
0 & 0& -1 & 2 & -1 \\\
0 & 0& 0 &-1 & 2 \\\
\end{pmatrix} X= \begin{pmatrix} 1 \\\
1 \\\
1 \\\
1 \\\
1 \\\
\end{pmatrix}. $$ La solution exacte de ce problème est $X = [2.5, 4,4.5, 4,2.5]^T$.

Classe Gauss-Seidel

Classe Relaxation