Équation de Poisson

Résolution numérique par différences finies

Y. Morel

Equation de Poisson

L'équation de Poisson (Siméon Denis - , 1781-1840) est la célèbre équation aux dérivées partielles du second ordre:

ΔΦ = f

où Δ est l'opérateur laplacien voir ci-dessous, et f une distribution donnée.
Poisson aboutit à cette équation en corrigeant celle de Laplace:

ΔΦ = 0

On retrouve l'équation de Poisson dans de nombreux domaines physiques. En électrostatique, le potentiel électrique V est solution de l'équation de Poisson

ΔΦ = -ρ / ε₀

dans laquelle ρ décrit la distribution de charges électriques dans l'espace, et ε₀ est une constante (permittivité diélectrique du vide).

En gravitation universelle, le potentiel gravitationnel Φ est relié à la masse volumique ρ et à la constante universelle de gravitation G par l'équation de Poisso

ΔΦ = 4πGρ

Laplacien

Le laplacien d'une fonction est la somme de ses dérivées partielles secondes. En 2D, le laplacien d'une fonction f de deux variables, f(x,y) est

Δf = ∂²f / ∂x² + ∂²f / ∂y²

Par exemple, les fonctions du 1er degré en chaque variable

f (x,y) = ax + bxy + cy + d

sont solutions de l'équation de Laplace Δf = 0.
Bien sûr, on n'a pas ici résolu l'équation de Laplace, car résoudre une équation, qu'elle soit algébriqe du 1er ou 2nd degré, ou différentielle (EDO) ou aux dérivées partielles (EDP), signifie toujours trouver toutes les solutions…
De plus, le domaine dans lequel on cherche à résoudre cette équation, et les conditions imposées sur le bord de ce domaine importent aussi.

Les paragraphes qui suivent ne visent pas à résoudre formellement et théoriquement cette équation (ni même d'ailleurs, ce qui constitue la première question à adresser, à déjà s'assurer que de telles solutions existent).
En général, on ne connaît pas d'expression explicite et analytique de la solution de ce type d'équation. On arrive néanmoins à en donner des valeurs approchées: c'est l'objet de ce qui suit.

Discrétisation des dérivées, première et seconde

Formules de Taylor-Young et approximation des dérivées

Les formules de Taylor sont les formules fondamentales pour approximer des dérivées et discrétiser des équations différentielles.

Pour une fonction d'une seule variable, de classe C¹, la formule de Taylor-Young à l'ordre 1 s'écrit:

f (x + h) = f (x) + hf'(x) + o(h)

soit, en isolant la dérivée,

f' (x) = f (x + h) − f (x) / h + o(1)

ce qui donne l'approximation, pour h petit,

f' (x) ≃ f (x + h) − f (x) / h

qui n'est rien d'autre que le taux d'accroissement de la fonction entre x et x + h.

À l'ordre 2, pour une fonction de classe C², la formule de Taylor-Young est

f (x + h) = f (x) + hf'(x) + h² / 2 f''(x) + o(h²)

On peut aussi écrire cette formule avec un accroissement de −h:

f (x − h) = f (x) + (− h)f'(x) + (−h)² / 2 f''(x) + o(h²) = f (x) − hf'(x) + h² / 2 f''(x) + o(h²)

et donc, en ajoutant ces deux dernières relations, on obtient

f (x + h) + f (x − h) = 2f (x) + h²f''(x) + o(h²)

ou encore, en isolant la dérivée seconde:

f'' (x) = f (x + h) + f (x − h) − 2f(x) / h² + o(1)

ce qui nous donne l'approximation, pour h petit,

f'' (x) ≃ f (x + h) + f (x − h) − 2f(x) / h²

Discrétisation et grille de calcul en 1D

Pour approximer une équation différentielle avec une fonction f d'une seule variable sur un intervalle [a, b], on discrétise l'intervalle en une suite d'abscisses x₀, x₁, …

$\psset{unit=1.2cm,arrowsize=6pt} \begin{pspicture}(-1,-1.5)(11,.5) \newcommand{\addu}[1]{#1 1 add} \newcommand{\addoc}[1]{#1 .5 add} \psline[linewidth=1.4pt,arrowsize=9pt]{->}(-.5,0)(10.5,0) \multido{\i=0+1}{3}{\psline(\i,-.2)(\i,.2)} \multido{\i=0+1}{2}{\psline(\i,-.2)(\i,.2) \psline{<->}(\i,-1)(!\addu{\i}\space-1) \rput(!\addoc{\i}\space-1.2){$h$}} \psline(8,-.2)(8,.2)\psline(9,-.2)(9,.2) \psline{<->}(8,-1)(9,-1)\rput(8.5,-1.2){$h$} \psline(4,-.2)(4,.2)\psline(5,-.2)(5,.2) \psline{<->}(4,-1)(5,-1) \rput(4.5,-1.2){$h$} \rput(4,-.4){$x_i$} \rput(5,-.4){$x_{i+1}$} \rput[r](0.2,-.4){$a=x_0$} \rput(1,-.4){$x_1$} \rput(2,-.4){$x_2$} \rput(3,-.4){$\dots$}\rput(3,-1){$\dots$} \rput(6.6,-.4){$\dots$}\rput(6.6,-1){$\dots$} \rput(8,-.4){$x_{n-1}$} \rput[l](8.8,-.4){$b=x_n$} \end{pspicture}$

On note alors f_i = f (x_i) et on a donc, avec la formule de Taylor-Young, une approximation de la dérivée

f' (x_i) ≃ f (x_i + h) − f (x_i) / h = f_i+1 − f_i / h

et de même de la dérivée seconde

f'' (x_i) ≃ f (x_i + h) + f (x_i − h) − 2f (x_i) / h² = f_i+1 + f_i−1 − 2f_i / h²

Discrétisation et grille de calcul en 2D

On procède de même en 2D, en discrétisant comme précédemment abscisses et ordonnées: avec un pas h en abscisse et un pas k en ordonnée:

$\psset{arrowsize=6pt}\begin{pspicture}(-2.3,-2.3)(10.5,10.5) \multido{\i=0+1}{2}{ \psline(\i,0)(\i,2.2) \psline[linestyle=dashed](\i,-1)(\i,8) \psline(\i,8)(\i,10) \psline(0,\i)(2,\i) \psline[linestyle=dashed](2,\i)(8,\i) \psline(8,\i)(10,\i)} \multido{\i=9+1}{2}{ \psline(\i,0)(\i,2.2) \psline[linestyle=dashed](\i,-1)(\i,8) \psline(\i,8)(\i,10) \psline(0,\i)(2,\i) \psline[linestyle=dashed](2,\i)(8,\i) \psline(8,\i)(10,\i)} % \psline(2,0)(2,1.2)\psline[linestyle=dashed](2,-1)(2,10) \psline(0,2)(1.2,2)\psline[linestyle=dashed](1.2,2)(10,2) \psline(5,0)(5,10)\psline(0,5)(10,5) \psline(4,0)(4,10)\psline(0,4)(10,4) \psline(6,0)(6,10)\psline(0,6)(10,6) \psline{<->}(0,-1.6)(1,-1.6)\rput(0.5,-2){$h$} \rput(0,-1.3){$x_0$}\rput(1,-1.3){$x_1$} \psline{<->}(1,-1.6)(2,-1.6)\rput(1.5,-2){$h$} \rput(2,-1.3){$x_2$} \psline{<->}(4,-1.6)(5,-1.6)\rput(4.5,-2){$h$} \psline{<->}(5,-1.6)(6,-1.6)\rput(5.5,-2){$h$} \rput(10,-1.3){$x_n$} \psline[linestyle=dashed](5,-1)(5,0)\rput(5,-1.3){$x_i$} \psline[linestyle=dashed](4,-1)(4,0)\rput(4,-1.3){$x_{i-1}$} \psline[linestyle=dashed](6,-1)(6,0)\rput(6,-1.3){$x_{i+1}$} % \psline[linestyle=dashed](-1,0)(10,0) \psline[linestyle=dashed](-1,1)(10,1) \psline[linestyle=dashed](-1,2)(10,2) \psline{<->}(-1.8,0)(-1.8,1)\rput(-2,.5){$k$} \rput(-1.3,0){$y_0$}\rput(-1.3,1){$y_1$} \psline{<->}(-1.8,1)(-1.8,2)\rput(-2,1.5){$k$} \rput(-1.3,2){$y_2$} \psline{<->}(-1.8,4)(-1.8,5)\rput(-2,4.5){$k$} \psline{<->}(-1.8,5)(-1.8,6)\rput(-2,5.5){$k$} \psline[linestyle=dashed](-1,10)(1,10)\rput(-1.3,10){$y_n$} \psline[linestyle=dashed](-1,5)(0,5)\rput(-1.3,5){$y_j$} \psline[linestyle=dashed](-1,4)(0,4)\rput(-1.3,4){$y_{j-1}$} \psline[linestyle=dashed](-1,6)(0,6)\rput(-1.3,6){$y_{j+1}$} \end{pspicture}$

Chaque point de la grille est maitenant repéré par un couple ( x_i , y_j ) ou encore, plus simplement par un couple d'entiers ( i , j )

On note maintenant, en utilisant cette grille

f_{i, j} = f (x_i , y_j)

et on approxime les dérivées secondes comme dans le cas à une variable, f(x,y) est

Δf (x_i , y_j) = ∂²f / ∂x² (x_i , y_j) + ∂²f / ∂y² (x_i , y_j) ≃ f_{i+1, j} + f_{i−1, j} − 2f_{i, j} / h² + f_{i, j+1} + f_{i, j−1} − 2f_{i, j} / k²

Exemple: potentiel électrostatique dans un condensateur

La charge électrostatique V dans un condensateur est solution du problème schématisé par la situation suivante:

$\begin{pspicture}(-.5,-.5)(4.5,2.5) \psline[linewidth=2.5pt](0,2)(3,2) \rput[l](3.2,2){$V=0$} \psline[linewidth=2.5pt,linecolor=red](0,0)(3,0) \rput[l](3.2,0){$V=1$} \rput(1.5,1){$\Delta V=0$} \end{pspicture}$

En reprenant la grille de discrétisation précédente, on cherche donc les valeurs

V_{i, j} = V (x_i , y_j)

telles que, de plus, pour tout i,

V_{i, n} = 0 V_{i, 0} = 1

pour les conditions à l'anode et la cathode du condensateur, et dans le vide (de charge) à l'intérieur du condensateur, l'équation de Laplace discrétisée

ΔV ≃ V_{i+1, j} + V_{i−1, j} − 2V_{i, j} / h² + V_{i, j+1} + V_{i, j−1} − 2V_{i, j} / k² = 0

Pour simplifier, on peut choisir le même pas pour les deux variables h = k ce qui permet de réécrire la relation précédente suivant:

V_i,j = 1/4 V_{i+1, j} + V_{i−1, j} + V_{i, j+1} + V_{i, j−1}

En d'autres termes, l'équation de Poission impose qu'en chaque point le potentiel soit la moyenne de ces quatre voisins ("voisins" au sens de Von Neumann, voir plus loin).

Programmation et simulation

Il s'agit donc de construire un tableau de valeurs, ou matrice, dont les éléments vérifient la relation précédente ainsi que les conditions au bord cf. .
Une méthode pour le faire est de construire une suite en itérant l'équation précédente, c'est-à-dire en définissant une suite V⁽ⁿ⁾

V_{i , j}⁽ⁿ⁺¹⁾ = 1/4 V_{i+1, j}⁽ⁿ⁾ + V_{i−1, j}⁽ⁿ⁾ + V_{i, j+1}⁽ⁿ⁾ + V_{i, j−1}⁽ⁿ⁾

En quelques itérations, la méthode converge.
Il y a donc ici 3 boucles à programmer: i et j pour balayer les éléments de la matrice V et enfin n pour le nombre d'itérations de la méthode.

Le programme peut s'écrire en python (en utilisant la bibliothèque Pylab):

from pylab import *

n=100
T=200
V=zeros([n,n])
# cathode, V=1
V[-1,:]=1

for t in range(T):
    for i in range(1,n-1):
        for j in range(1,n-1):
            V[i,j]=(V[i+1,j]+V[i-1,j]+V[i,j+1]+V[i,j-1])/4
    # cathode, V=1
    V[-1,:]=1


matshow(V)
show();

La fonction matshow utilisée finalement permet d'afficher le contenu de la matrice, c'est-à-dire exactement les valeurs du potentiel en tous les points de la grille. On obtient le résultat:

Dernier point: on remarque des "effets de bord" qui n'ont pas forcément lieu d'être. Pour modéliser une situation "infinie", on utilise souvent la condition au bord de Neumann (condition d'absorption), qui s'écrit ici, en Python,

from pylab import *

n=100
T=200
V=zeros([n,n])
# cathode, V=1
V[-1,:]=1

for t in range(T):
    for i in range(1,n-1):
        for j in range(1,n-1):
            V[i,j]=(V[i+1,j]+V[i-1,j]+V[i,j+1]+V[i,j-1])/4
    # cathode, V=1
    V[-1,:]=1
    # Neumann au bord
    V[:,0]=V[:,1]
    V[:,-1]=V[:,-2]


matshow(V)
show();

et on obtient cette fois le résultat:

Accélération de convergence

Dès que la grille de calcul devient importante, c'est-à-dire que l'on souhaite un résultat de plus en plus fin, le temps de calcul devient très important.

Deux méthodes (au moins) permettent d'accélérer les calculs.

Calculs vetoriels avec python

Tout comme en Matlab ou dans d'autres langages interprétés, les boucles imbriquées prennent un temps important de calcul: chaque instruction répétée et interprétée puis exécutée. Tout comme en Matlab (ou...) il est bien plus judicieux d'utiliser des instructions vectorielles, nécessitant alors de moins d'étapes iontermédiaires d'interprétation avant exécution.
En Python, comme en Matlab, on peut astucieusement remplacer les boucles imbriquées

for t in range(T):
    V[-1,:]=1
    tmp=V[2:n,1:n-1]+V[0:n-2,1:n-1]+V[1:n-1,2:n]+V[1:n-1,0:n-2]
    V[1:n-1,1:n-1]=1/4*tmp

Il ne resta alors plus qu'une boucle globale, et le temps de calcul, pour aboutir au même résultat est considérablement réduit:

Accélération de convergence par relaxation

On a accéléré indéniablement les calculs en utilisant des instructuions vectorielles plutôt que de multiples boucles imbriquées.
Maintenant ces calculs vectoriels restent éventuellement effectués un grand nombre de fois. On peut maintenant essayer de réduire la taille de cette dernière boucle. On se rappelle qu'on calcule ici une approximation, par itération successive: plus le nombre d'itération est grand, emilleure est l'approximation.
Ce calcul se présente comme celui de la recherche d'un point fixe, et à ce titre on peut lui appliquer la même méthode d'accélaration de convergence, par relaxation. Plus précisément, cette méthode est connue sous l'acronyme de SOR, pour Successive Over Relaxation: on utilise une méthode de relaxation qu'on applique successivement à chaque itération.
Le programme précédent, est modifié en:

w=1# paramètre de relaxation
T=100
for t in range(T):
    tmp=V[2:n,1:n-1]+V[0:n-2,1:n-1]+V[1:n-1,2:n]+V[1:n-1,0:n-2]
    V[1:n-1,1:n-1]=(1-w)*V[1:n-1,1:n-1]+w/4*tmp

Le paramètre de relaxation optimal ω_opt est difficile à déterminer, au contraire du cas plus simple décrit là.

Cette méthode est surtout utile si ce calcul doit lui-même être effectué un grand nombre de fois, c'est-à-dire lui-même imbriqué dans une boucle. Dans ce cas, on essaie simplement de déterminer une valeur efficace (même si elle n'est pas otpimale) de ce paramètre avec quelques tests.
Par exemple, dans la simulation numérique d'éclairs, ou le calcul du nouveau potentiel V doit être effectué un grand nombre de fois, cette méthode de relaxation est particulièremnet bienvenue.