Table des matières PDF

Équation de Poisson : résolution numérique

1. Introduction

L'équation de Poisson à deux dimensions est :

où u(x,y) est la fonction inconnue et s(x,y) la fonction source, éventuellement nulle (équation de Laplace).

Un exemple d'équation de Poisson est celle vérifiée par le potentiel électrostatique :

ρ est la densité volumique de charge électrique et ε la permittivité électrique du milieu.

Un autre exemple est l'équation de la chaleur en régime stationnaire vérifiée par la température :

σ(x,y) est la source thermique et λ la conductivité thermique.

On recherche une solution de l'équation de Poisson sur le domaine carré [0,1]x[0,1]. Sur les bords du carré, la condition limite appliquée consiste à fixer la valeur de u (condition limite de Dirichlet). La figure suivante montre un exemple de problème qui sera traité, avec une source nulle partout, et avec deux lignes à l'intérieur du domaine où les valeurs de u sont imposées (modèle de condensateur plan). Un exemple de maillage utilisé pour la discrétisation est aussi représenté.

../../../../figures/numerique/elliptique/poisson/domaine.svgFigure pleine page

Nous allons voir comment effectuer la discrétisation de l'équation de Poisson, puis comment résoudre le système d'équations obtenu avec une méthode itérative.

2. Discrétisation

2.a. Maillage

Le domaine carré est divisé avec un maillage dont les mailles sont carrées. On note h le pas d'espace du maillage. La figure suivante représente quelques mailles.

../../../../figures/numerique/elliptique/poisson/mailles.svgFigure pleine page

On recherche une solution approchée de la fonction u(x,y) aux nœuds du maillage. La valeur recherchée au point d'indices (i,j) est notée Ui,j. La valeur de la fonction source au même point est notée Si,j.

On a donc deux matrices U et S (on verra plus loin pour leur dimension).

2.b. Laplacien

Nous allons discrétiser l'équation de Poisson en utilisant la méthode des différences finies. Cette méthode consiste à remplacer les dérivées par des différences finies, afin d'obtenir des équations vérifiées par les valeurs Ui,j aux nœuds du maillage.

Pour discrétiser une dérivée seconde, on utilise les développements de Taylor suivants, pour une fonction f à une variable :

En sommant ces deux équations, on obtient :

Cette relation permet de discrétiser les dérivées secondes, par exemple :

En procédant de même pour l'autre dérivée, on obtient l'équation suivante pour chaque nœud du maillage :

Cette équation constitue un schéma de discrétisation de l'équation de Poisson. Nous allons l'écrire sous une forme plus générale :

Les matrices D (droite), G (gauche), H (haut), B (bas), C (centre) et S (source) ont les mêmes dimensions que le maillage. Par exemple, la matrice D comporte, pour chaque nœud du maillage, le coefficient par lequel il faut multiplier la valeur du nœud voisin situé à sa droite.

Pour un point non situé sur une frontière, on a donc :

2.c. Conditions limites de Dirichlet

Sur les bords du domaine, une valeur de u est imposée. Il peut aussi y avoir des lignes à l'intérieur du domaine où la valeur de u est imposée.

Pour un nœud où la valeur imposée est U0, l'équation s'applique à condition d'utiliser :

2.d. Dimensions du maillage

L'équation a l'avantage de s'appliquer à tous les nœuds du maillage. Il existe d'ailleurs d'autres conditions limites qui peuvent se mettre sous cette forme. Elle peut même s'appliquer à d'autres équations linéaires à dérivées partielles.

On voit cependant un problème sur les bords du maillage. Par exemple sur le bord inférieur, il n'y a pas de point situé en dessous (i,j-1). Pour éviter d'avoir à faire des tests pour savoir si le point est sur un bord, on ajoute des rangées fictives sur les bords. Ces rangées ne sont pas utilisées dans le calcul mais permettent d'adresser les nœuds en dehors du domaine.

Par ailleurs, le nombre de nœuds dans chaque dimension sera choisi égal à une puissance de 2 plus 1 :

p est un entier. Cela fait un nombre de mailles puissance de 2. L'avantage de ce choix est de permettre une division du nombre de mailles par 2, 4, etc. En comptant les rangées fictives, le nombre d'éléments des matrices dans chaque dimension est :

La figure suivante montre un exemple avec p=3, soit 8x8 mailles. Dans ce cas, les matrices U,G,D,H,B,F ont pour dimension (11,11). Les indices d'accès aux éléments de la matrice sont indiqués sur la figure. On voit ainsi que les indices des nœuds du maillage vont de 1 à N.

../../../../figures/numerique/elliptique/poisson/matrice.svgFigure pleine page

3. Méthode itérative

3.a. Système d'équations linéaires

Pour chacun des NxN nœuds du maillage, il y a une équation linéaire de la forme :

On doit donc résoudre un système linéaire de NxN équations à NxN inconnues. Si on place toutes les inconnues Ui,j dans une matrice colonne U (en les rangeant par exemple par lignes successives), le système peut se mettre sous la forme matricielle :

La matrice A est carrée, de taille NxN par NxN. Par exemple, pour N=65, c'est une matrice de 4225x4225, soit environ 18 millions d'éléments. Il est donc hors de question de chercher à résoudre ce système en utilisant un programme standard utilisant la méthode d'élimination de Gauss. Comme la matrice est très creuse, il faudrait programmer une méthode d'élimination spécifique pour ce type de matrice.

Les méthodes de résolution utilisées pour ce type de matrice ne sont pas des méthodes directes comme la méthode d'élimination de Gauss, mais des méthodes itératives.

3.b. Méthode de Gauss-Seidel

La méthode de Gauss-Seidel est une méthode itérative ([1]). Une méthode itérative consiste à partir d'une matrice U(0) comportant des valeurs initiales nulles, puis à appliquer une relation de récurrence pour calculer une suite de matrices U(k). La relation de récurrence doit être construite pour que la suite converge vers la solution approchée du système d'équations :

La méthode de Gauss-Seidel peut se définir en utilisant le formalisme matriciel, ce qui permet d'en faire l'étude théorique ([2]). D'un point de vue pratique, on utilise la relation , c'est-à-dire que l'on garde la numérotation des inconnues à deux indices, correspondant à la grille de discrétisation.

Voyons comment la matrice U(k+1) est calculée à partir de la matrice U(k). La matrice est parcourue par lignes successives, en partant de la ligne j=1 (on utilise la numérotation précisée plus haut, avec les rangées fictives). Pour chaque ligne, l'indice de colonne i va de 1 à N. Le parcours complet de la matrice est appelé itération de Gauss-Seidel. Chaque élément de la matrice est modifié en utilisant ses 4 voisins, de la manière suivante :

On voit que la modification de Ui,j fait appel à deux points déjà calculés au cours de la même itération : le point situé à sa gauche et le point situé en dessous. La relation de récurrence s'écrit donc :

Pour l'implémentation, on utilise plutôt la relation . La figure suivante montre une itération en cours. Les cercles pleins indiquent les nœuds déjà traités par l'itération. Le nœud (i,j) est en cours de traitement.

../../../../figures/numerique/elliptique/poisson/gaussSeidel.svgFigure pleine page

On admet que la méthode de Gauss-Seidel converge pour l'équation de Poisson ([2]), c'est-à-dire qu'elle tend vers une matrice limite. Cette limite est précisément la solution (approchée) recherchée.

Pour suivre la convergence de la méthode de Gauss-Seidel, on calculera la norme suivante :

qui tend vers une limite finie.

Pour un point où la source est nulle, et non situé sur un bord où la valeur est imposée, la transformation s'écrit :

L'élément de la matrice est donc remplacé par la moyenne de ses quatre voisins. Cette opération peut s'apparenter à une opération de filtrage d'image (filtre laplacien).

3.c. Accélération par sur-relaxation

La méthode de Gauss-Seidel converge très lentement : le nombre d'itérations nécessaire pour obtenir une bonne approximation est proportionnel au nombre de points du maillage Q=N2. Comme y a Q points à calculer pour chaque itération, cela donne un temps de convergence proportionnel à Q2.

La méthode de sur-relaxation permet d'accélérer considérablement la convergence. D'une manière générale, une méthode de relaxation consiste à modifier la méthode de Gauss-Seidel de la manière suivante :

ω 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 :

Il existe une valeur optimale du paramètre de relaxation, qui permet d'obtenir une convergence beaucoup plus rapide que la méthode de Gauss-Seidel, avec un nombre d'itérations proportionnel à N (la racine carrée du nombre d'inconnues). Cette valeur optimale dépend de l'équation discrétisée, des conditions limites, et de la taille du maillage. Pour l'équation de Poisson avec des conditions limites de Dirichlet sur les bords, l'expression du paramètre optimal est connue :

Lorsque N augmente, la valeur s'approche de 2. Par exemple, pour N=65, la valeur optimale est ω*=1,91. Pour cette valeur de N le gain par rapport à la méthode de Gauss-Seidel est déjà d'un facteur 65, ce qui signifie qu'un calcul complet durera une minute au lieu d'une heure. La complexité globale de la sur-relaxation optimale est en Q3/2, à comparer à la complexité en Q2 de la méthode d'élimination de Gauss.

4. Travaux pratiques

4.a. Implémentation de la méthode de Gauss-Seidel avec sur-relaxation

import numpy
import math
from  matplotlib.pyplot import *
            

Les matrices U,G,D,H,B,F sont des tableaux numpy de taille (M,M).

Pour faciliter l'implémentation, on définit une classe, dont on donne ici le constructeur :

class Poisson:
    def __init__(self,p):
        self.N = 2**p+1
        self.M = M = self.N+2
        self.h = 1.0/(self.N-1)
        self.C = numpy.zeros((M,M),dtype=numpy.float64)
        self.G = numpy.zeros((M,M),dtype=numpy.float64)
        self.D = numpy.zeros((M,M),dtype=numpy.float64)
        self.H = numpy.zeros((M,M),dtype=numpy.float64)
        self.B = numpy.zeros((M,M),dtype=numpy.float64)
        self.F = numpy.zeros((M,M),dtype=numpy.float64)
        self.U = numpy.zeros((M,M),dtype=numpy.float64)
        self.R = numpy.zeros((M,M),dtype=numpy.float64)
        self.omega = 2.0/(1.0+math.sin(math.pi*self.h))
            

Le constructeur calcul les nombres de points sur chaque dimension en fonction de l'entier p. Il crée les matrices et calcule aussi le pas d'espace h et le paramètre de relaxation optimal.

Pour les fonctions de la classe, on donne seulement l'entête :

            
    def laplacien(self):   
            

effectue la discrétisation du laplacien, c'est-à-dire remplit les matrices C,G,D,H,B,F avec une source nulle (Fi,j=0).

             
    def valeur_bords(self,v):   
             

fixe la condition limite sur les bords, ce qui consiste à attribuer les bonnes valeurs de C,G,D,H,B,F sur les bords du domaine.

              
    def valeur_ligne_horiz(self,x,y,longueur,v):        
              

fixe une condition limite de Dirichlet (valeur v imposée) sur une ligne horizontale située dans le domaine. Les coordonnées (x,y) du point de départ et la longueur sont données dans l'intervalle [0,1].

              
    def source_ligne_horiz(self,x,y,longueur,v):        
              

définit une source sur une ligne horizontale, ce qui consiste à modifier les valeurs de F sur cette ligne.

              
    def source_ponctuelle(self,x,y,v):
              

définit une source ponctuelle, qui occupe un seul nœud du maillage.

    
    def norme(self):
              

calcule et renvoie la norme de la matrice U (équation ).

              
    def gauss_seidel(self,w1,w):
              

effectue une itération de Gauss-Seidel avec sur-relaxation, c'est-à-dire un parcours complet de la matrice U pour modifier ses valeurs selon le schéma . w est la valeur du paramètre de sur-relaxation ω et w1 est 1-ω.

               
    def iterations(self,ni,w):
               

effectue ni itérations de Gauss-Seidel avec sur-relaxation. Après chaque itération, la norme de U est affichée avec print.

4.b. Exemples

Condensateur plan

La géométrie a été définie en introduction. Il s'agit de deux plaques infinies dans la direction Z.

Voici un exemple avec N=65. Une centaine d'itérations suffit.

poisson = Poisson(6)
poisson.laplacien()
poisson.valeur_bords(0.0)
poisson.valeur_ligne_horiz(0.25,0.4,0.5,1.0)
poisson.valeur_ligne_horiz(0.25,0.6,0.5,-1.0)
poisson.iterations(100,poisson.omega)
matshow(poisson.U,cmap='gray',vmin=-1.0,vmax=1.0)
        
figAfigA.pdf
figure(figsize=(7,7))
contour(poisson.U,20)
        
figBfigB.pdf

Fils chargés

Un fil chargé (parallèle à l'axe Z) est représenté par une source ponctuelle. Une source ponctuelle positive correspond à une charge négative. On pourra par exemple traiter le cas de deux fils de charges opposées :

poisson = Poisson(6)
poisson.laplacien()
poisson.valeur_bords(0.0)
poisson.source_ponctuelle(0.3,0.5,1)
poisson.source_ponctuelle(0.7,0.5,-1)
            
Références
[1]  W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery,  Numerical recipes, the art of scientific computing,  (Cambridge University Press, 2007)
[2]  P. Lascaux, R. Théodor,  Analyse numérique matricielle appliquée à l'art de l'ingénieur,  (Dunod, 1994)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.