L'équation de Poisson

De Laplace à Poisson

Dans une page précédente, nous avons étudié l'équation de Laplace et sa résolution numérique par des méthodes aux différences finies. Cette équation, dont la forme générale est \( \Delta V = 0 \) permet, entre autres, de calculer le potentiel créé par une répartition de charges électriques externes dans un domaine fermé vide de charge. Les domaines d'application de cette EDP elliptique homogène sont multiples : mécanique des fluides, thermique et même analyse financière.

Dans la présente page, nous allons examiner une équation très proche de l'équation de Laplace : l'équation de Poisson. C'est aussi une équation aux dérivées partielles elliptique, de forme laplacienne, dont l'expression générale est \( \Delta V = f(x_0,..,x_i) \). Plus précisément, je vais aborder la résolution numérique de cette équation, dans une de ses formes particulières, qui est \( \Delta V = K \), avec K une constante non nulle bien sur !

Un peu de physique

L'équation de Poisson

Imaginons une région de l'espace où il existe une distribution de charges \( \rho(x,y) \). Cette distribution de charges produit un champ électrique dans le domaine fermé lequel nous nous positionnons pour notre étude. L'équation de Maxwell-Gauss devient donc \( div\vec{E} = \dfrac{\rho(x,y)}{\epsilon_0} \). Dans cette équation, remplaçons \( \vec{E} \) par son expression en fonction du potentiel V, nous obtenons \( -div(\vec{grad}V) = \dfrac{\rho(x,y)}{\epsilon_0} \) ou, ce qui revient au même \( div \:\vec{grad}V = -\dfrac{\rho}{\epsilon_0} \). C'est l'équation de Poisson, au encore appelée par les physiciens l'équation de Maxwell-Gauss, sous sa forme locale.

Dans la pratique, on utilise une autre notation, en employant l'opérateur laplacien et qui s'exprime par \( \Delta \: V = div(\vec{grad}V)\). Notre équation de Poisson s'écrit donc \( \Delta \: V = -\dfrac{\rho(x,y)}{\epsilon_0} \).

Son expression en coordonnées cartésiennes

Dans la suite de cette page, pour simplifier, nous nous placerons dans un plan. Dans ce plan, le laplacien d'un potentiel scalaire V, comme le potentiel électrique, s'exprime par \( \Delta V = \dfrac{\partial^2V}{\partial x^2} + \dfrac{\partial^2V}{\partial y^2} \).

L'équation de Poisson devient \( \dfrac{\partial^2V}{\partial x^2} + \dfrac{\partial^2V}{\partial y^2} = -\dfrac{\rho(x,y)}{\epsilon_0} \). C'est cette équation que nous allons résoudre numériquement.

Vous constaterez qu'il s'agit d'une équation elliptique, avec des conditions de Dirichlet, qui se résoud analytiquement assez simplement par la méthode de la séparation des variables. Ici, nous allons la résoudre numériquement avec la méthode de Gauss-Seidel déjà vue par ailleurs.

Résolution numérique de l'équation de Poisson

La physique du problème

Soit deux charges, +Q et -Q, disposées sur une surface fermée vide dont les bords sont maintenus à un potentiel constant nul. Le problème consiste à calculer le potentiel créé sur cette surface par notre distribution de charges.

La discrétisation de l'équation de Poisson 2D

La discrétisation de l'espace

Comme pour l'équation de Laplace, nous allons utiliser les méthodes aux différences finies, que j'ai abordé dans cette page.

Dans notre cas, cela revient à mailler le plan sur lequel nous voulons résoudre l'équation de Poisson, par une grille dont les mailles sont très petites, de forme rectangulaires ou carrée, de dimension \( \Delta x\) et \( \Delta y\).

La discrétisation de l'équation

Nous allons discrétiser notre équation en réalisant un développement de Taylor d'ordre de nos deux dérivées partielles. Nous obtenons:

\( \dfrac{\partial^2 V(x,y)}{\partial x^2} = \dfrac{V(x+\Delta x,y) - 2V(x,y) + V(x-\Delta x,y)}{(\Delta x)^2}\)

\( \dfrac{\partial^2 V(x,y)}{\partial y^2} = \dfrac{V(x,y + \Delta y) - 2V(x,y) + V(x,y - \Delta y)}{(\Delta y)^2}\)

ce qui nous donne la forme discrétisée de l'équation de Poisson, que nous utiliserons maintenant :

\( \dfrac{V(x+\Delta x,y) - 2V(x,y) + V(x-\Delta x,y)}{(\Delta x)^2} + \dfrac{V(x,y + \Delta y) - 2V(x,y) + V(x,y - \Delta y)}{(\Delta y)^2} = -C(x,y)\)

La matrice du terme de droite représente la distribution discréte des charges sur le domaine de résolution du problème.

La méthode de Gauss-Seidel

Le principe

Dans la forme discrète de l'équation de Poisson, posons pour simplifier \( \Delta x = \Delta y = \Delta n \), ce qui nous donne :
\( \dfrac{V(x+\Delta x,y) - 2V(x,y) + V(x-\Delta x,y) + V(x,y + \Delta y) - 2V(x,y) + V(x,y - \Delta y)}{(\Delta n)^2} = -C(x,y)\)

soit en regroupant les termes:

\( \dfrac{ -4V(x,y) + V(x+\Delta x,y) + V(x-\Delta x,y) + V(x,y + \Delta y) + V(x,y - \Delta y)}{(\Delta n)^2} = -C(x,y)\)

Ce principe est identique à celui que nous avons retenu dans le cas de l'équation de Laplace, au terme de droite près...

Sympoliquement, je peux donc écrire :
V[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V[i,j+1] + V[i,j-1] + C[i,j])

Et comme il s'agit d'une méthode de relaxation, je parcours tous les points intérieurs de la grille autant de fois que nécessaire pour que la différence entre la valeur du potentiel en chaque point de la grille entre deux itérations soit inférieure à une quantité que j'aurais fixée, qui sera la précision de mon calcul.

Le script

La première partie du script Poisson2D.py fixe les constantes de calcul et les constantes physiques et construit la grille V dont on aura besoin pour les calculs. Cette partie n'attire aucune remarque particulère.

Puis je définie les conditions aux limites et les conditions initiales à l'intérieur de la grille, car je vous rappelle que nous sommes en présence d'un problème de Dirichlet. le code est le suivant:

V[0,:] = V0 # bord supérieur

V[:,0] = V0 # bord gauche

V[:,-1] = V0 # bord droit

V[-1,:] = V0 # bord inférieur

pour les conditions aux limites de la grille. Les cotés de la grille sont au potentiel nul. Notez la notation vectorielle utilisée pour éviter l'usage de boucles.

et pour les conditions initiales à l'intérieur de la grille, au potentiel nul :

V[1:N,1:N] = V0

La matrice C, initialisée à 0, contient la répartition des charges sur le domaine de calcul. Ici, en l'occurence, je place une charge Q positive dans le premier quadrant du domaine, et une charge négative -Q dans le troisième quadrant du domaine.

C = zeros([N+1,N+1])

C[N/4,N/4] = Q

C[3*N/4,3*N/4] = -Q

Suit la boucle de relaxation dont le code est :

while ecart > EPS:

iteration += 1

Vprec = V.copy()

V[1:-1,1:-1]= 0.25*(Vprec[0:-2,1:-1]+V[2:,1:-1]+Vprec[1:-1,0:-2]+V[1:-1,2:]+C[1:-1,1:-1])

ecart = np.max(np.abs(V-Vprec))

La boucle de relaxation tournera tant que la précision déterminée par EPS n'est pas atteinte. La variable ecart, le critère de convergence, sera calculée dans la boucle. Notez dans la boucle le compteur d'itérations et aussi, avant et après la boucle, l'acquisition de l'heure pour déterminer le temps de calcul (fonction time()).

Le reste du code sert à l'affichage de la grille et ne présente pas grand intérêt...

Les résultats

Avec le code ci-dessus, j'obtiens les résultats suivants:

Le nombre d'itérations pour atteindre la précision demandée (10-3) est de 3060.
Le temps de calcul est d'environ une seconde sur mon Precision M6400.

Sur le plan physique, le potentiel dans le domaine en fonction de la position des charges s'établit comme suit :

On pourrait vérifier par quelques calculs simples que la loi de Coulomb pour l'électrostatique est vérifiée.

Les scripts Python

Les scripts Python étudiés dans cette page sont disponibles dans le package PoissonPython.zip :

Poisson2D.py : résolution de l'équation de Poisson en utilisant la méthode de Gauss-Seidel

Pour conclure

Avec un peu de pratique, l'utilisation des méthodes aux différences finies pour résoudre numériquement des EDP se révèle souple et assez puissante, du moins dans nos cas très simples. Vous pouvez vous entrainer en modifiant la répartition des charges ou bien le maillage de la grille, par exemple en le resserrant à proximité des charges.

S'agissant du potentiel créé par un système de charges discrètes, on peut remarquer que la résolution numérique ne dit pas grand chose du potentiel à proximité des charges, surtout lorsqu'on tend vers la charge. D'après la loi Coulomb, on tendrait vers l'infini, ce qui constitue une singularité. Que se passe-t-il à proximité immédiate de la charge, d'un électron par exemple ? Et d'ailleurs, la question a-t-elle un sens, à savoir qu'est-ce que la proximité d'un électron ? Je me penche sur le sujet dans cette page.