Méthode de Gauss Seidel

1. Introduction

La discrétisation des équations différentielles à dérivées partielles linéaires à conditions limites conduisent à la résolution d'une équation linéaire de la forme :

$$AU=F\text{(1)}$$

où U est une matrice colonne contenant les inconnues, c'est-à-dire les valeurs de la fonction inconnue aux N points du maillage. Voir par exemple le cas de l'équation de Poisson bidimensionnelle en coordonnées cartésiennes.

La matrice carrée A, qui comporte NxN éléments, est très largement creuse (typiquement 5 diagonales) et comporte un très grand nombre de termes non nuls (typiquement de l'ordre de 10⁴ pour un problème 2D). Les méthodes de résolution directes, comme la méthode d'élimination de Gauss, sont très peu efficaces pour ce type de système, et peuvent conduire à des erreurs de calcul importantes. Les méthodes itératives sont largement préférables dès que N est grand. La méthode de Gauss-Seidel est une méthode itérative, que ce document présente succinctement, avant d'en étudier empiriquement la convergence sur un exemple simple. On verra aussi la méthode de Gauss-Seidel avec sur-relaxation, qui permet d'accélérer considérablement la vitesse de convergence.

2. Méthode de Gauss-Seidel

On s'intéresse plus particulièrement aux équations à deux dimensions. Les inconnues sont donc représentées sur un maillage (ou une grille) à deux dimensions. On note U_i,j l'inconnue au nœud (i,j) du maillage.

On considère un schéma numérique à 5 points. La figure suivante représente un nœud du maillage avec les 4 voisins qui interviennent dans le schéma numérique :

Figure pleine page

Le schéma numérique, c'est-à-dire l'équation qui résulte de la discrétisation du nœud (i,j) se met sous la forme :

$$D_{i,j}{U}_{i+1,j}+G_{i,j}{U}_{i-1,j}+H_{i,j}{U}_{i,j+1}+B_{i,j}{U}_{i,j-1}+C_{i,j}{U}_{i,j}=F_{i,j}\text{(2)}$$

Les lettres utilisées signifient Droite, Gauche, Haut, Bas et Centre.

Cette équation générale s'applique également aux points situés sur des frontières. Pour cela on entoure les bords du domaine de rangées de points fictifs qui n'interviennent pas dans le calcul. Par exemple, pour le bord gauche du domaine (indice i=0), on introduit une rangée de points fictifs d'indice i=-1 et on pose G_0,j=0. Il y a donc 6 matrices de (N_x)x(N_y) points à stocker en mémoire.

On a donc l'équations suivante pour chaque nœud du maillage :

$$U_{i,j}=\frac{1}{C_{i,j}}(F_{i,j}-D_{i,j}{U}_{i+1,j}-G_{i,j}{U}_{i-1,j}-H_{i,j}{U}_{i,j+1}-B_{i,j}{U}_{i,j-1})\text{(3)}$$

La méthode de Gauss-Seidel ([1][2]) consiste à effectuer des itérations. Au départ, les valeurs U_i,j sont initialisées à une valeur nulle, où aux valeurs correspondant à une fonction proche de la solution si cela est possible. Lors d'une itération, le maillage est parcouru par lignes successives. Le point (i,j) est modifié comme suit :

$$U_{i,j}\leftarrow \frac{1}{C_{i,j}}(F_{i,j}-D_{i,j}{U}_{i+1,j}-G_{i,j}{U}_{i-1,j}-H_{i,j}{U}_{i,j+1}-B_{i,j}{U}_{i,j-1})\text{(4)}$$

Cette transformation consitue le schéma itératif de Gauss-Seidel.

Si le maillage est parcouru par i croissant et par j croissant, le changement du point (i,j) fait intervenir les points (i-1,j) et (i,j-1) déjà calculés au cours de l'itération. L'itération fonctionne donc sur un seul tableau de stockage. Les valeurs obtenues à l'itération k en fonction de celles de l'itération k-1 sont alors :

$$U_{i,j}^{\left(k\right)}=\frac{1}{C_{i,j}}(F_{i,j}-D_{i,j}{U}_{i+1,j}^{(k-1)}-G_{i,j}{U}_{i-1,j}^{\left(k\right)}-H_{i,j}{U}_{i,j+1}^{(k-1)}-B_{i,j}{U}_{i,j-1}^{\left(k\right)})\text{(5)}$$

Le schéma de Gauss-Seidel peut s'écrire sous forme matricielle. Si l'on note U^(k) la matrice colonne des valeurs des points du maillage obtenus à la k-ième itération, on a :

$$U^{(k+1)}=P{U}^{\left(k\right)}+G\text{(6)}$$

où P est la matrice d'itération. En introduisant l'erreur définie par :

$$E^{\left(k\right)}=A{U}^{\left(k\right)}-F\text{(7)}$$

on obtient :

$$E^{\left(k\right)}=P{E}^{(k-1)}=P^2{E}^{(k-2)}=\cdots =P^k{E}^{\left(0\right)}\text{(8)}$$

La méthode de Gauss-Seidel converge si l'erreur tend vers zéro. Il y a convergence si les valeurs propres de la matrice d'itération ont toutes un module strictement inférieur à 1 ([3]). Lorsqu'il a convergence, la valeur limite vérifie bien l'équation que l'on cherche à résoudre :

$$A{U}^{\left(\infty \right)}=F\text{(9)}$$

Une variante de la méthode consiste à traiter séparément les nœuds pairs et les nœuds impairs. Pour cela, on utilise un nombre de mailles égal à une puissance de 2 dans chaque direction. Par exemple, pour 2³ mailles dans chaque direction, on obtient le maillage suivant :

Figure pleine page

Chaque itération est constituée de deux passes : une pour les nœuds pairs, l'autre pour les nœuds impairs. Cela revient à renuméroter les inconnues de la matrice U. Cette méthode est appelée méthode de Gauss-Seidel rouge-noir (par référence aux couleurs d'un échiquier).

Le principal avantage de cette méthode est de permettre un traitement parallèle des nœuds pairs (ou des nœuds impairs). La calcul sur les nœuds pairs peut être effectué par plusieurs processus de calculs (ou threads) sur un CPU. On peut aussi facilement programmer l'itération sur un processeur graphique (GPU), en répartissant les calculs pour ces nœuds sur autant d'unités de calcul. Il faut bien sûr penser à synchroniser les unités de calculs (ou les threads sur CPU) à la fin du traitement des nœuds pairs, avant de commencer celui des nœuds impairs.

3. Étude de la convergence

On considère l'équation de Laplace :

$$\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=0\text{(10)}$$

avec une condition de Dirichlet u=0 sur les 4 bords du domaine. La solution est u(x,y)=0.

Les dimensions de la maille sont Δx et Δy. Pour alléger les notations, on considère des mailles carrées et on pose :

$$h=\Delta x=\Delta y\text{(11)}$$

On adopte un domaine carré [0,1]x[0,1] avec le même nombre de mailles N_x=2ⁿ dans les deux directions. Le pas de discrétisation est donc :

$$h=\frac{1}{2^n}\text{(12)}$$

La condition initiale est un mode de Fourier, c'est-à-dire une sinusoïde de la forme :

$$u(x,y)=sin\left(2\pi mx\right)\text{(13)}$$

La convergence est suivie au moyen de la norme définie par :

$$\left|\right|U\left|\right|=\sqrt{\frac{1}{N}\sum _{i,j}{U}_{i,j}^2}\text{(14)}$$

from pylab import *
import numpy
import poisson.main
             

Définition du maillage, discrétisation du Laplacien et conditions limites :

n=7
p=4
levels = n-p+1
laplace=poisson.main.Poisson(p,p,levels,1.0,1.0)
laplace.laplacien()
laplace.dirichlet_borders(0.0)
             

Condition initiale : mode m=10

laplace.sinus_mode(10,1)
             

Itérations de Gauss-Seidel et tracé de la norme en fonction du nombre d'itérations :

result=laplace.opencl_iterations_norm(10,100)
figure(figsize=(10,5))
plot(result[0],result[1])
xlabel('niter')
ylabel('norm')
             

plotA.pdf

Même chose pour le mode m=6

laplace.sinus_mode(6,1)
result=laplace.opencl_iterations_norm(10,100)
figure(figsize=(10,5))
plot(result[0],result[1])
xlabel('niter')
ylabel('norm')
             

plotB.pdf

Pour le mode m=2

laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(10,100)
laplace.close()
figure(figsize=(10,5))
plot(result[0],result[1])
xlabel('niter')
ylabel('norm')

             

plotC.pdf

On constate que la convergence est rapide pour les modes élevés (hautes fréquences spatiales), très lente pour les modes bas. Voyons comment converge le mode m=2 avec une maille deux fois plus grande :

n=6
p=4
levels = n-p+1
laplace=poisson.main.Poisson(p,p,levels,1.0,1.0)
laplace.laplacien()
laplace.dirichlet_borders(0.0)
laplace.sinus_mode(2,1) 
result=laplace.opencl_iterations_norm(10,100)
laplace.close()
figure(figsize=(10,5))
plot(result[0],result[1]) 
xlabel('niter')
ylabel('norm')
             

plotD.pdf

La convergence est beaucoup plus rapide sur le maillage de pas double, en terme de nombre d'itérations. En terme de quantité de calculs, il faut aussi tenir compte du fait qu'il y a 4 fois moins de points, ce qui fait un gain d'au moins 10 pour ce mode.

On voit donc qu'un maillage très fin (nécessaire pour obtenir des détails) est très efficace pour la relaxation des hautes fréquences spatiales, mais rend la convergence des basses fréquences très lentes. Pour relaxer les basses fréquences, on a intérêt à utiliser un maillage plus grossier. Cette propriété est mise à profit dans les méthodes multigrilles ([1]), qui consiste à utiliser successivement des grilles de tailles variables.

De manière générale, le nombre d'itérations nécessaires augmente avec le nombre de points N du maillage. Plus précisément, il augmente proportionnellement à N. La complexité globale de la méthode est donc N², ce qui en fait une méthode très peu efficace pour les maillages de grande dimension.

4. Accélération par sur-relaxation

Les méthodes multigrilles sont très efficaces mais sont difficiles à implémenter. La méthode d'accélération par sur-relaxation ([2]), très simple à mettre en œuvre, permet d'accélérer la convergence en modifiant l'itération de Gauss-Seidel de la manière suivante.

La valeur d'un nœud est modifiée en faisant une combinaison linéaire de l'ancienne valeur et de celle donnée par le schéma de Gauss-Seidel :

$$U_{i,j}\leftarrow (1-\omega )U_{i,j}+\omega \left(\frac{1}{C_{i,j}}\right(F_{i,j}-D_{i,j}{U}_{i+1,j}-G_{i,j}{U}_{i-1,j}-H_{i,j}{U}_{i,j+1}-B_{i,j}{U}_{i,j-1}\left)\right)\text{(15)}$$

où ω est le paramètre de relaxation. Si 0<ω<1, la nouvelle valeur est une moyenne pondérée de l'ancienne et de celle donnée par le schéma de Gauss-Seidel. Dans ce cas, la convergence est ralentie (sous-relaxation). La sur-relaxation consiste à choisir une valeur de ω strictement supérieure à 1. Pour que la méthode converge, il faut que :

$$1<\omega <2\text{(16)}$$

La méthode de Gauss-Seidel avec sur-relaxation est souvent nommée méthode de sur-relaxation (SOL, succesive over-relaxation).

Qualitativement, la méthode de sur-relaxation consiste à modifier la valeur en allant plus loin dans le sens indiqué par le schéma de Gauss-Seidel. Il existe une valeur optimale de ω qui permet d'obtenir un nombre d'itération proportionnel à N^1/2 ([1]), ce qui fait une complexité globale N^3/2. Cette valeur optimale s'approche de 2 lorsque N augmente. La difficulté de cette méthode réside dans l'obtention de la valeur optimale de ω. Théoriquement, on peut l'obtenir si l'on connait les valeurs propres de la matrice de relaxation du schéma de Gauss-Seidel. Plus précisément, il faut connaître le rayon spectral ρ_gs, c'est-à-dire le module maximal des valeurs propres, qui doit être inférieur à 1 pour que la méthode de Gauss-Seidel converge. La valeur optimale est alors ([1]) :

$${\omega }^{*}=\frac{2}{1+\sqrt{1-{\rho }_{gs}}}\text{(17)}$$

Par exemple, pour l'équation de Poisson bidimensionnelle avec des conditions limites de Dirichlet sur les bords d'un domaine carré, le rayon spectral de la méthode de Gauss-Seidel avec parcours par lignes est ([1]):

$${\rho }_{gs}={cos}^2\left(\frac{\pi }{N_x}\right)\text{(18)}$$

La valeur optimale est dans ce cas :

$${\omega }^{*}=\frac{2}{1+sin\left(\frac{\pi }{N_x}\right)}\simeq \frac{2}{1+\frac{\pi }{N_x}}\simeq 2-\frac{2\pi }{N_x}\text{(19)}$$

Poson N_x=2ⁿ et voyons le tracé de la valeur optimale en fonction de n :

import math
n = numpy.arange(4,12)
w = numpy.zeros(n.size)
for i in range(w.size):
    w[i] = 2.0/(1.0+math.pi/(2**n[i]))
figure()
plot(n,w,'o')
xlabel('n')
ylabel('w')
grid()
                 

plotD1.pdf

On voit par exemple pour un maillage de 256x256 une valeur optimale ω^*=1.976.

Lorsque des conditions limites supplémentaires interviennent dans le domaine, ou lorsque le schéma de discrétisation est différent (par exemple en symétrie axiale), la valeur optimale est modifiée. On peut néanmoins se servir de la valeur donnée ci-dessus comme point de départ. Si l'on souhaite optimiser, il faut faire des essais avec différentes valeurs de ω.

Comme exemple, on reprend le cas précédent en testant différentes valeurs de ω :

n=6
p=4
levels = n-p+1
laplace=poisson.main.Poisson(p,p,levels,1.0,1.0)
laplace.laplacien() 
laplace.dirichlet_borders(0.0)
figure(figsize=(10,5))
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.0)
plot(result[0],result[1],label='1.0')
laplace.init_array()
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.1)
plot(result[0],result[1],label='1.1')
laplace.init_array()
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.3)
plot(result[0],result[1],label='1.3')
laplace.init_array()
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.5)
plot(result[0],result[1],label='1.5')
laplace.init_array()
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.7)
plot(result[0],result[1],label='1.7')
laplace.init_array()
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.9)
plot(result[0],result[1],label='1.9')
xlabel('niter')
ylabel('norm')
legend(loc='upper right')

               

plotE.pdf

Pour cet exemple, la valeur ω=1.9 donne la donne la convergence la plus rapide. Voyons ce qu'il arrive pour des valeurs supérieures à 1.9 :

figure(figsize=(10,5))
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.92)
plot(result[0],result[1],label='1.92')
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.95)
plot(result[0],result[1],label='1.95')
laplace.sinus_mode(2,1)
result=laplace.opencl_iterations_norm(5,100,omega=1.98)
plot(result[0],result[1],label='1.98')
xlabel('niter')
ylabel('norm')
legend(loc='upper right')
laplace.close()
               

plotF.pdf

La convergence n'est plus monotone, et plus lente que pour ω=1.9. La valeur optimale de ω est donc proche de 1.9 pour cet exemple, conformément à ce que prévoit la théorie pour ce problème.