Vous avez rencontré pour la première fois ce concept en 1ere S (programme de 2011). Vous avez appris que le champ électrostatique (ou champ électrique) était un champ vectoriel, et que sa définition était \( \vec{E} = \frac{\vec{F}}{q} \). Il s'agit de l'expression du champ électrique en un point M de l'espace. \( \vec{F} \) est la force qu'exerce ce champ sur un "corps d'épreuve", une petite boule de sureau traditionnellement, chargé d'une charge électrique q. Son unité usuelle est le Volt.m-1, alors que l'unité SI standard est le N.C-1.
Le cours d'électrodynamique classique de John David Jackson, la bible de l'électromagnétisme classique de mon temps, nous donne cette définition du champ électrique on doit définir le champ [électrique] comme la limite du rapport de la force exercée sur le corps d'épreuve à la charge portée par ce corps, lorsque cette charge devient de plus en plus petite. Et de préciser en note de bas de page, que physiquement la charge ne peut pas devenir infiniment petite, à cause de la nature discrète de celle-ci. Millikan a démontré cette nature discrète et a même mesuré cette charge élémentaire (voir ici).
La loi de Coulomb nous permet d'exprimer la force en question (vue aussi en 1ere S !), qui s'exerce entre 2 charges q et q1 distantes d'une distance r, par \( \vec{F} = kqq1\frac{\vec{r}}{|\vec{r}|^3} \), la constante k étant sans intérêt pour l'instant. Ce qui nous permet d'écrire, en vertue de la définition ci-dessus, qu'en un point M où est située la charge q, \( \vec{E} = kq1\frac{\vec{r}}{|\vec{r}|^3} \). C'est l'expression du champ électrique E créé par la charge q1, en un point M.
Le champ électrique, comme tous les champs qui sont en 1/r^2, possède une caractéristique essentielle : il dérive d'un potentiel, ce qui mathématiquement s'exprime par une relation fondemantale : \( \vec{E} = -\vec{grad}V\), V étant le potentiel électrique dont vous avez plus ou moins entendu parlé. Nous reviendrons sur cette notion dans une future page. En attendant, retenez cette relation, qui est valable pour tous les champs dits "à flux conservatif", comme le champ électrique ou le champ gravitationnel.
Quelques valeurs numériques pour fixer les idées : la constante k est égale à \( \frac{1}{4 \pi\epsilon_0} = 10^{-7}\) avec \( \epsilon_0\) la permittivité du vide. Le champ électrique produit par un électron en un point situé à une distance de 1m de cet électron est égal à environ 10^-9 V.m-1.
Vous noterez que je considère ici un champ constant, ne dépendant pas du temps (électro_statique). Je m'en tiendrai à cette contrainte dans la suite de cette page.
Dans cette page, nous avons découvert le théorème de Gauss dans son application au champ gravitationnel. Je vous avais dit qu'il était généralisable à tous les champs en 1/r^2, et donc au champ électrique. Si nous l'appliquons dans le cas du champ éléctrique, nous obtenons \( \iint\limits_{(S)} \vec{E}.\vec{n} dS = \frac{q}{\epsilon_0} \), q étant la charge électrique du volume limité par la surface S, qui peut être quelconque.
Allons un peu plus loin en appliquant un théorème d'analyse vectorielle que tous les étudiants de prépa et de licence doivent connaître par coeur, ainsi que son orthographe, j'ai cité le théorème d'Ostrogradsky (ou Green-Ostrogradsky), qui nous donne que sous certaines conditions toujours remplies par hypothèse, \( \iint\limits_{(S)} \vec{E}.\vec{n} dS = \iiint\limits_{(V)} div\vec{E}dV \).
Par une suite de calculs assez simples dont je vous fais grâce (ce sont des exo standards de prépa), à partir de l'application des théorèmes de Gauss et d'Ostrogradsky, nous obtenons l'équation dite de Maxwell-Gauss \( div\vec{E} = \frac{\rho}{\epsilon_0} \), qui permet de calculer le champ électrique en présence d'une densité volumique de charge, pour une région de l'espace donnée, dans le vide.
Imaginons une région de l'espace où cette densité volumique de charge est localement nulle. Le champ électrique est produit par une distribution de charges située en dehors de la région dans laquelle nous nous positionnons pour notre étude. L'équation de Maxwell-Gauss devient donc \( div\vec{E} = 0 \). Dans cette équation, remplaçons \( \vec{E} \) par son expression en fonction du potentiel V, nous obtenons \( -div(\vec{grad}V) = 0 \) ou, ce qui revient au même \( div \:\vec{grad}V = 0 \). C'est l'équation de Laplace.
Dans la pratique, on utilise une autre notation, en employant un opérateur qui s'appelle le laplacien (évidemment !) et qui s'exprime par \( \Delta \: V = div(\vec{grad}V)\). Notre équation de Laplace s'écrit donc \( \Delta \: V = 0 \), ce qui est assez élégant vous en conviendrez...
Pour la petite histoire, si la densité de charges n'est pas nulle, nous obtenons l'équation de Poisson \( \Delta \: V = -\frac{\rho}{\epsilon_0} \), que nous étudierons dans une autre page.
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 = \frac{\partial^2V}{\partial x^2} + \frac{\partial^2V}{\partial y^2} \).
L'équation de Laplace devient \( \frac{\partial^2V}{\partial x^2} + \frac{\partial^2V}{\partial y^2} = 0\). C'est cette équation que nous allons résoudre numériquement.
Si vous vous référez à la classification des équations aux dérivées partielles que j'ai listé ici, 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. Mais ce n'est pas notre objectif ! Nous allons étudier différentes méthodes pour la résoudre numériquement.
Pour étudier les différentes méthodes de résolution, je vais choisir un problème très simple d'électrostatique planaire. Je considère un carré vide, bordé par trois cotés conducteurs reliés à la terre, donc par référence au potentiel 0V. Le quatrième coté, lui aussi conducteur est isolé des trois autres cotés par deux cales isolantes de dimension négligeable. ce coté est porté à un potentiel Vi = 100 V. pour fixer les idées, notre dispositif est représenté sur le schéma ci-dessus.
Le problème consiste à visualiser le champ électrique créé dans le carré par la distribution de charges linéique portée le coté porté au potentiel Vi. Vous aurez noté que le carré dans lequel nous ferons le calcul est vide de charge, et donc que l'équation de Laplace s'applique à ce dispositif.
Comme pour l'équation de la chaleur, une autre équation aux dérivées partielles, 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 Laplace, par une grille dont les mailles sont très petites, de forme rectangulaires ou carrée, de dimension \( \Delta x\) et \( \Delta y\).
Comme pour l'équation de la chaleur, 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:
\( \frac{\partial^2 V(x,y)}{\partial x^2} = \frac{V(x+\Delta x,y) - 2V(x,y) + V(x-\Delta x,y)}{(\Delta x)^2}\)
\( \frac{\partial^2 V(x,y)}{\partial y^2} = \frac{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 notre équation de Laplace, que nous utiliserons maintenant :
\( \frac{V(x+\Delta x,y) - 2V(x,y) + V(x-\Delta x,y)}{(\Delta x)^2} + \frac{V(x,y + \Delta y) - 2V(x,y) + V(x,y - \Delta y)}{(\Delta y)^2} = 0\)
La méthode des différences finies appliquée à la résolution des EDP conduit à résoudre un système d'équations linéaires de NxN variables, N étant le nombre de points sur un coté de la grille et en supposant que pour simplifier on choisisse une grille carrée avec un pas identique en x et en y.
Il est possible de traiter le problème, i.e. de résoudre le système, en utilisant le calcul matriciel, ce qui se fait pour les grandes valeurs de N (plusieurs milliers). Mais il existe d'autres techniques très efficaces pour les petites valeurs de N, quelques centaines typiquement. Ce sont les méthodes de relaxation, comme la méthode de Jacobi, la méthode de Gauss-Seidel, la méthode SOR (Simultaneous OverRelaxation) ou d'autres encore. Ces méthodes sont très faciles à programmer. Elles ont en outre l'avantage d'être stables numériquement.
Dans notre cas, nous allons comparer la méthode de Jacobi et la méthode de Gauss-Seidel, qui sont les deux méthodes de relaxation les plus simples.
Il est très simple : pour chaque point intérieur de la grille, hors les bords, on pose que la valeur du potentiel en un point M[i,j] est égal à la moyenne des potentiels de ses quatre voisins. Ce que j'écris de manière symbolique:
V[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V[i,j+1] + V[i,j-1])
Et comme il s'agit d'une méthode de relaxation, ou dite encore méthode itérative, 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.
Bien sur, on ne pose pas cette égalité au hasard ! On l'obtient à partir de la forme discrète de l'équation de LapLace que nous avons établi plus haut en posant \( \Delta x = \Delta y\) et en réarrangeant les termes.
La spécificité de la méthode de Jacobi tient dans le fait que le potentiel V[i,j] est calculé dans l'itération k de la relaxation à partir des valeurs de potentiel calculées lors de l'itération k-1, ce que l'on exprime par :
\( V^k[i,j] = 0.25*(V^{k-1}[i-1,j] + V^{k-1}[i+1,j] + V^{k-1}[i,j+1] + V^{k-1}[i,j-1]) \)
Il s'agit du script Python LaplaceJacobi.py.
La première partie du script fixe les constantes de calcul et les constantes physiques et construit les deux grilles V et Vnew 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,:] = Vi # 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. Le bord supérieur est porté au potentiel Vi et les autres cotés 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
Suit la boucle de relaxation dont le code est :
Vnew = V.copy()
while ecart > EPS:
iteration += 1
# boucle de calcul en python pur
for i in range(1,N):
for j in range(1,N):
Vnew[i,j] = 0.25*(V[i-1,j] + V[i+1,j] + V [i,j-1] + V[i,j+1])
# critère de convergence
ecart = np.max(np.abs(Vnew - V))
# sauvegarde dans la grille V de la grille calculée
V = Vnew.copy()
Le code manipule deux grilles : la grille V qui contient les potentiels calculés lors de l'itération k-1 de la grille et la grille Vnew dans laquelle on stockera les potentiels calculés dans l'itération courante k. La première instruction de cette partie du code initialise la grille Vnew avec le contenu de la grille V, ici en l'occurence les valeurs initiales et aux limites.
Puis on entre dans la boucle de relaxation qui 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.
Vous rencontrez maintenant le coeur de l'algorithme que l'on a étudié plus haut. Vous constaterez que je calcule bien la nouvelle valeur du potentiel, que je stocke dans Vnew[i,j], à partir des données de l'itération précédente stockées dans V.
Le critère de convergence est ensuite calculé, en utilisant une norme entre les deux matrices V et Vnew. Puis, je transfère les données de Vnew dans V, avant de passer à l'itération k+1.
Le reste du code sert à l'affichage de la grille et ne présente pas grand intérêt...
La fonction time() permet d'acquérir le temps au début et à la fin de la boucle de calcul, afin de déterminer le temps de calcul.
Avec le code ci-dessus et sa double boucle, j'obtiens les résultats suivants:
Sur le plan physique, le potentiel de la planque en fonction de la position s'établit comme suit :
Comme nous pouvions nous y attendre, les régions de haut potentiel sont situées à proximité du bord supérieur de la plaque alors que les cotés, le centre et la partie inférieure de la plaque est à un potentiel faible ou nul.
Python est un langage interprété, qui donc à ce titre n'aime pas trop les boucles... Mais Python permet un autre traitement des tableaux et matrices en utilisant un traitement vectoriel, codé directement en C, qui est disponible avec NumPy. J'ai modifié mon code dans le script LaplaceJacobiNumpy.py en remplaçant la double boucle par l'instruction suivante :
V[1:-1,1:-1]= 0.25*(V[0:-2,1:-1] +V[2:,1:-1] + V[1:-1,0:-2] + V[1:-1,2:])
En exécutant ce script, j'obtiens une grille identique à celle obtenue précédemment.
Le nombre d'itérations pour atteindre la précision demandée (10-3) est identique : 6073, ce qui n'est pas surprenant. Par contre, le temps de calcul est très inférieur : 1,2 secondes sur le Precision M6400 et 0,7 secondes sur le MacBook.
L'amélioration est impressionnante : il faut toujours utiliser les particularités de notre langage de programmation. Ici, il s'agit d'utiliser les règles de vectorisation de Python au lieu de faire usage de boucles ! Cela serait aussi vrai en MatLab et Scilab.
La méthode de Gauss-Seidel est une variante de la méthode de Jacobi. Vous avez remarqué que dans la méthode de Jacobi, on calcule le potentiel au point M[i,j] dans l'itération k à partir des valeurs de potentiel calculées dans l'itération k-1. Dans la méthode de Gauss-Seidel, on mixte les données déjà calculées dans l'itération k (celles en V[i-1,j] et V[i,j-1]), qui sont plus à jour que les données de l'itération précédente, avec les données de l'itération précédente, qui ne sont pas encore calculées dans l'itération courante (celles en V[i+1,j] et V[i,j+1]). En somme on utilise les données les plus à jour !
En théorie, la méthode de Gauss-Seidel converge plus vite, mais je dois avouer que je n'ai jamais vu de grande différence pour les équations que nous traitons..
Il est identique à celui du script de Jacobi, sauf pour ce qui concerne la boucle de calcul, qui se simplifier un peu. Elle devient :
while ecart > EPS:
iteration += 1
# sauvegarde de la grille courante pour calculer l'écart
Vprec = V.copy()
# calcul de la nouvelle valeur du potentiel sur la grille
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:])
# critère de convergence
ecart = np.max(np.abs(V-Vprec))
Avec le code ci-dessus j'obtiens les résultats suivants, identiques à ceux obtenus avec Jacobi:
Le tracé du champ de potentiel est lui identique à celui obtenu par Jacobi, ce qui est rassurant:
Pour illustrer la mise en oeuvre de l'équation de Laplace et sa résolution numérique, j'ai choisi un grand classique des exercices d'électrostatique : le condensateur plan. Imaginons une enceinte cubique dont les parois conductrices sont reliées à la terre, donc au potentiel nul. Au centre de cette enceinte, on dispose deux plaques conductrices espacées d'une distance d, et portées l'une à un potentiel positif +Vp et l'autre à un potentiel négatif -Vp. L'exercice consiste à déterminer le champ de potentiel créé dans cette enceinte par ce dispositif.
Le dispositif expérimental, vu de dessus, est schématisé ainsi:
Nous sommes en présence d'un problème de Dirichlet : le potentiel se détermine à l'intérieur de l'enceinte et hors des plaques par l'équation de Laplace. Les conditions aux limites sont déterminées par le potentiel de l'enceinte et des deux plaques formant le condensateur.
J'ai choisi d'utiliser la méthode de Gauss-Seidel vectorisée, vu les performances obtenues ci-dessus. La méthode de Jacobi pourrait être également utilisée.
Le script Python s'appelle LaplaceCondensateur2D.py. Après les initialisations habituelles, nous fixons les conditions initiales et les conditions aux limites par les instructions suivantes:
# définition des conditions aux limites
V[0,:] = V0 # bord supérieur
V[:,0] = V0 # bord gauche
V[:,-1] = V0 # bord droit
V[-1,:] = V0 # bord inférieur
# initialisation de l'intérieur de la grille
V[1:N,1:N] = 0.0
# définition de la géométrie du condensateur
V[25:75,40] = -Vp # bord plaque gauche potentiel négatif
V[25:75,60] = Vp # bord plaque droit potentiel positif
Je fixe le potentiel sur les parois de l'enceinte, les bords de la grille de calcul puis à l'intérieur de la grille. Enfin, je fixe le potentiel pour chaque plaque du condensateur.
La boucle de relaxation est reprise de la méthode de Gauss-Seidel :
while ecart > EPS:
iteration += 1
# sauvegarde de la grille courante pour calculer l'écart
Vprec = V.copy()
# calcul de la nouvelle valeur du potentiel sur la grille
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:])
# je réaffecte les valeurs constantes sur les plaques
V[25:75,40] = -Vp
V[25:75,60] = Vp
# critère de convergence
ecart = np.max(np.abs(V-Vprec))
Notez toutefois une particularité, qui concerne le traitement du potentiel sur les plaques. Celui-ci ne varie pas, bien sur. Aussi, après chaque calcul du potentiel, je réaffecte la valeur initiale de ce potentiel sur les plaques. Ne pas oublier sous peine d'obtenir des résultats étranges !
L'affichage est identique aux précédents, à l'exception du code suivant que j'ai ajouté pour tracer les courbes équipotentielles :
x = linspace(0,N+1,N+1)
y = linspace(0,N+1,N+1)
X,Y = meshgrid(x,y)
equ = plt.contour(X,Y,V,20,colors='b')
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
plt.clabel(equ, fontsize=10, inline=1,fmt='%1.0f')
Voilà la forme que prend le champ de potentiel autour de notre condensateur dans sa boite:
Plusieurs choses à remarquer et à analyser sur cette figure:
A titre d'exercice, les plus courageux pourront tenter (et réussir) de calculer en chaque point le vecteur champ électrique et tracer les lignes de champ, en se souvenant que \( \vec{E} = -\vec{grad}V \).
Le nombre d'itérations pour atteindre la précision demandée (10-3) est de 3435. Le temps de calcul est de 0,76 secondes sur le Precision M6400 et 0,34 secondes sur le MacBook.
Les scripts Python étudiés dans cette page sont disponibles dans le package LaplacePython.zip :
La résolution numérique de l'équation de Laplace nous a permis d'aborder les méthodes itératives de solution des EDP elliptiques, méthodes puissantes et simples à mettre en oeuvre. Nous les retrouverons dans l'étude de l'équation de Poisson, qui est une généralisation de l'équation de Laplace pour un espace non vide de charge.
Accessoirement, le code contient des techniques d'affichage en Python qui seront utiles pour d'autres expériences.
Contenu et design par Dominique Lefebvre - tangenteX.com février 2016 -- Pour me joindre par ou PhysiqueX ou 
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.