Imaginons une fonction réelle quelconque f(p), qui associe à tout point p de l'espace physique un scalaire. Par exemple, elle pourrait associer une pression à chaque point de l'atmosphère terrestre. Cette fonction définit ce qu'on appelle un champ scalaire. On rencontre beaucoup de ces champs en physique: les champs de pression, de température, de potentiel, de densité, etc.
Imaginons maintenant que notre fonction f(p) n'associe plus un scalaire mais un vecteur à chaque point de l'espace physique. On pourra considérer par exemple, que chaque point de l'espace est le point d'application de ce vecteur. Les champs vectoriels sont tout aussi nombreux : les champs électriques, les champs de vitesses dans un fluide, etc.
Une notion importante est attachée à celle de champ vectoriel: il s'agit de la notion de ligne de champ. Une ligne de champ est toujours tangente en chaque point de l'espace au vecteur champ.
Python nous permet de visualiser un champ vectoriel, par exemple le champ électrique créé par une charge. Je vous propose le script OPDChampElectrique.py, qui trace le champ créé par une charge unitaire positive.
Après les imports de rigueur, je définis le calcul du champ dans la routine ChampElectrique :
def ChampElectrique(q,A,x,y):
Ex = q*(x - A[0])/((x - A[0])**2 + (y - A[1])**2)**(1.5)
Ey = q*(y - A[1])/((x - A[0])**2 + (y - A[1])**2)**(1.5)
return Ex,Ey
Rien d'extraordinaire, je calcule les composantes x et y du champ en appliquant la définition que vous avez tous appris en électrostatique.
Puis je définis la grille de calcul, comme vous savez le faire maintenant:
X = linspace(-1,1,20)
Y = linspace(-1,1,20)
x,y = meshgrid(X,Y)
Nous verrons plus loin les particularités de la fonction meshgrid() qu'il faut connaitre, sous peine de quelques déboires !
Je calcule le champ en utilisant la routine définie ci-dessus, en plaçant la charge à l'origine de la grille:
Ex, Ey = ChampElectrique(1,[0,0],x,y)
Enfin, j'utilise la fonction quiver() de matplotlib.pyplot pour tracer le champ, sans fioriture. A noter toutefois la normalisation des Ex et Ey, car sinon les flèches représentatives des vecteurs à proximité de la charge sont démesurément grandes !
quiver(x,y,Ex/sqrt(Ex**2 + Ey**2),Ey/sqrt(Ex**2 + Ey**2))
Et voilà ce que cela donne :
Vous pouvez bien sur utiliser cette méthode pour tracer tout type de champ vectoriel. Pour agrémenter la visualisation, je vous invite à consulter l'aide en ligne de MatPlotLib pour la fonction quiver().
Vous avez l'habitude de désigner chaque point de l'espace physique par un triplet de réels (x,y,z) qui donnent une référence de position du point par rapport à une origine et une base unitaire. (x,y,z) forment les coordonnées cartésiennes du point.
Vous savez également calculer la distance entre deux points de cet espace, repérés par leurs coordonnées cartésiennes dans le même référentiel en appliquant une formule que vous avez apprise par coeur \( d = \sqrt((x- x_0)^2 + (y - y_0)^2 + (z - z_0)^2) \). Cette formule est une métrique, sans doute la plus utilisée. Il en existe d'autres, mais ce n'est pas le sujet. Ce qu'il faut retenir, c'est que quelque soit le système de coordonnées utilisé, polaire, sphérique ou autre, la distance entre deux points de l'espace ne change pas. Plus généralement, les lois de la physique ne dépendent pas du système de coordonnées que vous utilisez, pas plus que l'existence des points de l'espace ne dépend de ses coordonnées.
Comme le choix d'un système de coordonnées relève plutôt d'aspects pratiques, essentiellement la simplification des calculs, je donnerai la définition des principaux opérateurs différentiels dans différents systèmes, cartésiens, cylindriques et sphériques.
Rassurez-vous, je ne vais pas partir dans des grandes définitions mathématiques, mais seulement vous donner une idée ce ce qu'est un opérateur différentiel. Considérons l'ensemble des fonctions réelles à variables réelles, avec un nombre indéterminé de variables. Et imaginons une application qui fasse correspondre à une fonction de cet ensemble, une autre fonction de cet ensemble. Si j'appelle cette application G, j'aurais par exemple g(x) = G(f(x)). Et bien, on appelle simplement G un opérateur. Et si en plus f(x) est une fonction différentiable, alors G sera un opérateur différentiel. Pour mémoire, une fonction différentiable est une fonction à plusieurs variables, dérivable pour chacune de ses variables. c'est l'extension de la dérivabilité appliquée à des fonctions à plusieurs variables.
Bien sur, ma définition n'est pas très rigoureuse sur le plan mathématique, mais ça donne une idée. Et d'ailleurs voyons cela sur un exemple.
Vous avez tous déjà utilisé un opérateur différentiel ! Posons par exemple \(f(x) = \sin x\) et essayez d'imaginer l'application G qui lie la fonction \( \sin x \) et la fonction \( \cos x \). Cette dernière est évidemment la dérivée de la première, ce que l'on écrit \( \frac{d}{dx} \sin x = \cos x\). Ici notre application G est l'opérateur de dérivation \( \frac{d}{dx} \). Vous utiliserez toujours cet opérateur pour calculer la dérivée par rapport à x d'une fonction réelle quelconque, sous réserve bien sur qu'elle soit dérivable ! D'ailleurs, tous les opérateurs différentiels que nous aborderons sont des opérateurs de dérivation.
Bien sur, cet exemple peut être étendu à l'opérateur de dérivée partielle noté \( \frac{\partial}{\partial x}\), que nous allons utiliser très largement !
C'est le plus simple des opérateurs différentiels. Il se note \( \nabla \), s'appelle nabla en français (en grec plutôt...) et "Del" en anglais (un Delta à l'envers).
Sa définition est très simple: si je considère une base dans \( \mathbb{R}^3\) \( (\vec{e_x}, \vec{e_y}, \vec{e_z}) \), l'opérateur nabla se définit en coordonnées cartésiennes par \( \vec{\nabla} = \vec{e_x} \frac{\partial}{\partial x} + \vec{e_y} \frac{\partial}{\partial y} + \vec{e_z} \frac{\partial}{\partial z}\).
Bon , attention : j'écris \( \vec{\nabla} \) comme un vecteur, mais en fait ce n'est pas un vrai vecteur, il n'en a pas toutes les propriétés. Par exemple, dans notre espace vectoriel habituel, le produit scalaire de deux vecteurs est commutatif, mais ce n'est pas le cas du "produit scalaire" \( \vec{\nabla}.\vec{A} \not= \vec{A}.\vec{\nabla} \) Mais c'est pratique et c'est un usage assez courant chez les physiciens ! Evidemment, comme d'habitude, les mathématiciens font des bonds ! Il suffit de garder à l'esprit nos quelques approximations...
Dans ce système de coordonnées, la base devient \( (\vec{e_r}, \vec{e_{\theta}}, \vec{e_z}) \). Et la définition de l'opérateur nabla \( \vec{\nabla} = \vec{e_r} \frac{\partial}{\partial r} + \frac{1}{r}\vec{e_{\theta}} \frac{\partial}{\partial \theta} + \vec{e_z} \frac{\partial}{\partial z}\).
Dans ce système de coordonnées, la base devient \( (\vec{e_r}, \vec{e_{\theta}}, \vec{e_{\phi}}) \). Et la définition de l'opérateur nabla \( \vec{\nabla} = \vec{e_r} \frac{\partial}{\partial r} + \frac{1}{r}\vec{e_{\theta}} \frac{\partial}{\partial \theta} + \frac{1}{r \sin\theta} \vec{e_{\phi}} \frac{\partial}{\partial \phi}\).
L'opérateur nabla a cette caractéristique intéressante que de pouvoir exprimer tous les opérateurs différentiels que nous rencontrerons plus loin. Si j'appelle f un champ scalaire et \( \vec{A} \) un champ vectoriel, j'ai :
Ces quelques égalités devraient vous aider grandement à mémoriser les définitions des opérateurs ! je vous avais dit que c'était pratique de considérer cet opérateur comme un vecteur !
Considérons un champ scalaire f(p) dérivable dans l'espace, qui associe une grandeur scalaire à chaque point p de l'espace, par exemple un champ de température. Imaginons que nous voulions connaitre la variation de cette grandeur scalaire, la différence de température par exemple, entre deux points très proches de l'espace, p et p + ds, où ds est un déplacement très petit entre les deux points (d'où la nécessité que le champ soit dérivable...).
Par définition, la variation df est donnée par \( df = f(p+ds) - f(p) = \vec{grad} \: f.d\vec{s} \), c'est à dire par le produit scalaire du vecteur \( \vec{grad} \: f\), le gradient de notre champ et de la distance parcourue. Vous noterez que la définition du gradient est tout à fait indépendante du système de coordonnées choisi pour repérer p dans l'espace.
Dans la littérature, l'opérateur gradient est aussi noté \( \vec{\nabla} \), ce qui se justifie lorsqu'on compare la définition des deux opérateurs en coordonnées cartésiennes !
En coordonnées cartésiennes, on a \( \vec{grad} \: f = \frac{\partial f}{\partial x}(x,y,z)\vec{e}_x + \frac{\partial f}{\partial y}(x,y,z)\vec{e}_y + \frac{\partial f}{\partial z}(x,y,z)\vec{e}_z \), si f est un champ scalaire dérivable dans l'espace.
En coordonnées cylindriques, avec la base cylindrique définie plus haut, on a \( \vec{grad} \: f = \frac{\partial f}{\partial r}(r,\theta,z) \vec{e_r} + \frac{1}{r}\frac{\partial f}{\partial \theta}(r,\theta,z) \vec{e_{\theta}} + \frac{\partial f}{\partial z}(r,\theta,z) \vec{e_z} \).
En coordonnées sphériques, avec la base sphérique définie plus haut, on a \( \vec{grad} \: f = \frac{\partial f}{\partial r}(r,\theta,\phi) \vec{e_r} + \frac{1}{r}\frac{\partial f}{\partial \theta}(r,\theta,\phi) \vec{e_{\theta}} + \frac{1}{r \sin\theta}\frac{\partial f}{\partial \phi}(r,\theta,\phi) \vec{e_{\phi}} \).
Le vecteur gradient indique en chaque point d'un champ scalaire, la direction et l'intensité de la variation du champ. Il peut être utile de noter que la direction du vecteur gradient est la direction de la plus forte variation du champ scalaire f.
Sur un champ scalaire, le vecteur gradient est toujours perpendiculaire à une surface de niveau, encore nommée surface équipotentielle.
En effet, une surface équipotentielle est une surface sur laquelle la valeur du champ scalaire est constante, autrement dit la variation du champ f est nulle, et donc son gradient est nul par définition (le gradient est une dérivation, et la dérivée d'une fonction constante est nulle...). Si je reprends la définition \( \vec{grad} \: f.d\vec{s} = 0 \), cela signifie par définition du produit scalaire que \( \vec{grad} \: f \) et \( d\vec{s} \) sont perpendiculaires.
Bien sur, cette remarque est vraie pour les courbes équipotentielles !
Pour illustrer l'utiisation du gradient en calcul numérique avec python, je vais m'appuyer sur une définition connue de tous les élèves de prépa et de licence, qui lie le champ électrique et le potentiel électrostatique : \( \vec{E} = - \vec{grad} \: V \). Je vais donc définir un potentiel puis calculer son gradient. Le script OPDGradientChamp.py va faire cela pour nous.
def Potentiel(q,A,x,y):
V = q/((x - A[0])**2 + (y - A[1])**2)
return V
Ey,Ex = gradient(Potentiel(1,[0,0],x,y))
Ex = -Ex
Ey = -Ey
Le reste du code n'appelle pas de commentaire.
La visualisation du champ obtenu s'effectue comme ci-dessus, et j'obtiens :
Vous constatez, oh surprise, que nous obtenons la même figure de champ que ci-dessus.
En coordonnées cartésiennes, avec les notations habituelles, la divergence d'un champ de vecteurs \( \vec{A}\) s'exprime ainsi :
\( div(\vec{A}) =\vec{\nabla}.\vec{A} = \frac{\partial A_x}{\partial x} + \frac{\partial A_y}{\partial y} + \frac{\partial A_z}{\partial z}\)
Attention à une confusion que l'on rencontre parfois, y compris sur internet : il s'agit bien de la somme des dérivées partielles de chaque composante du vecteur \( \vec{A}\) et pas des dérivées partielles du vecteur \( \vec{A}\). Pour être plus clair:
\( div(\vec{A}) \not= \frac{\partial A(x,y,z)}{\partial x} + \frac{\partial A(x,y,z)}{\partial y} + \frac{\partial A(x,y,z)}{\partial z}\)
ça l'air bête, mais le nombre de fois que j'ai rencontré la confusion est hallucinant, surtout dans les forums d'aide...
En coordonnées cylindriques, avec la base cylindrique définie plus haut, on a :
\( div(\vec{A}) = \frac{1}{r}\frac{\partial}{\partial r}(r.A_r) + \frac{1}{r}\frac{\partial A_{\theta}}{\partial \theta} + \frac{\partial A_z}{\partial z} \).
En coordonnées sphériques, avec la base sphérique définie plus haut, on a:
\( div(\vec{A}) = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2.A_r) + \frac{1}{r\sin{\theta}}\frac{\partial}{\partial \theta}(\sin\theta A_{\theta}) + \frac{1}{r\sin{\theta}}\frac{\partial A_{\Phi}}{\partial \Phi} \).
Intuitivement, vous percevez sans doute la signifie d'une grandeur divergente : elle s'éloigne de sa valeur initiale. Considérons un champ quelconque dans un volume d'espace de dimension 3. Ce champ vectoriel est créé par une source ou il peut s'agir d'un champ de vitesse, en hydrodynamique par exemple. Dans ce volume, il peut se passer trois choses :
La valeur de la divergence d'un champ en un volume élémentaire nous donne aussi une indication sur l'allure des lignes de champ en ce volume. Si la divergence est positive, nous sommes en présence d'une source de champ et les lignes divergent de ce volume. Si la divergence est négative, nous sommes en présence d'un puit et les lignes de champ convergent vers ce volume. Si la divergence est nulle, alors les lignes de champs ont parallèles entre elles, elles entrent et sortent du volume sans modification. En général...
En général, parce que on peut observer des champs de divergence non nulle, dont les lignes de champ sont parallèles : il suffit que le champ varie en module dans le sens de \( \vec{A} \). On peut aussi observer des champs de divergence nulle avec des lignes de champ qui ne sont pas parallèles entre elles : il suffit par exemple que \( \frac{\partial A_x}{\partial x} = -\frac{\partial A_y}{\partial y} \) en dimension 2.
Enfin, une dernière chose, que nous avons abordé dans la page sur la singularité électronique à propos de source "ponctuelle". Si la source était vraiment ponctuelle, sa densité serait infinie puisque son volume tendrait vers 0. Sa divergence serait alors infinie positive, ce qui ferait désordre dans nos équations. Mais heureusement, une source ponctuelle n'existe pas !
Examinons le cas des deux champs des équations de Maxwell. L'équation de Maxwell-Gauss pose que \( div \: \vec{E} = \frac{\rho}{\epsilon_0} \), c'est à dire que la divergence du champ électrique dépend de la charge volumique du volume élémentaire. Elle peut être nulle en l'absence de charge, positive en cas de charge positive ou de défaut de charges négatives et négative en présence de charge négative. Alors que l'équation de Maxwell-Thomson pose que \( div \: \vec{B} = 0 \). Il n'existe pas de source isolée de champ magnétique, ce qu'on appelle des monopôles magnétiques, mais seulement des dipôles. Le champ magnétique est un champ à flux conservatif.
Autre cas intéressant, le champ gravitationnel. En physique classique (hors relativité générale), sa divergence est donnée par l'équation \( div \: \vec{G} = -4 \pi G \rho \), proportionnelle à l'opposé de la distribution de masse. Cette divergence est toujours négative (en présence de masse) ou nulle (en l'absence de masse). Une masse placée dans un champ gravitationnel se comportera comme un puit pour ce champ.
On l'appelle encore théorème de la divergence de Gauss, et cela vous explique pourquoi j'éprouve la nécessité d'en parler ici. Il dit que pour un volume V donné (on reste dans l'espace quotidien!) fermé par une surface S, l'intégrale de la divergence du champ \( \vec{A} \) sur V est égale à l'intégrale du flux de \( \vec{A} \) à travers la surface fermée S. Exprimé mathématiquement, cela nous dit:
\( \iiint\limits_V \vec{\nabla}.\vec{A} dV = \iint\limits_S \vec{A} dS\)
On retrouve dans ce théorème la définition de la divergence : il nous dit que si la divergence d'un champ dans un volume V donné est non nulle, alors le flux à travers la surface du volume est égale à cette divergence. Par exemple, si le volume V renferme une charge électrique \( \rho \), alors cette charge engendrera un flux à travers S, que vous pourrez calculer à l'aide du théorème de Gauss. De la même manière, si le bilan de flux à travers S est nul (ce qui rentre sort), on peut affirmer que la divergence du champ est nulle dans ce volume.
Il nous permet ainsi de proposer une définition intrinsèque de la divergence, indépendante d'un système de coordonnées:
\( div \: \vec{A} = \lim_{\Delta \tau \rightarrow 0} \frac{\Delta \Phi}{\Delta \tau} \)
où \( \Phi \) est le flux du champ \( \vec{A} \) et \( \Delta \tau \) le volume élémentaire dans lequel on calcule la divergence.
Le calcul numérique avec python de la divergence d'un champ est relativement facile en partant de la définition de la divergence en coordonnées cartésiennes. J'ai choisi un champ simple \( \vec{E} = \frac{x^2}{2} \vec{e_x} - xy \vec{e_y}\). En faisant le calcul manuellement, vous constaterez que sa divergence est nulle. Voyons le script OPDChampDivergence.py qui effectue ce calcul et trace la divergence de ce champ.
Commençons par coder les équations de définition du champ, pour sa composante Ex et Ey, rien de sorcier :
def Champ(x,y):
Ex = x*x*0.5
Ey = -x*y
return Ex,Ey
Passons maintenant à l'essentiel, la routine de calcul de la divergence :
def Divergence(Fx,Fy):
div = (diff(Fx,axis=1)/hx)[:-1,:] + (diff(Fy,axis=0)/hy)[:,:-1]
return div
Cette routine tient sur une ligne mais mérite quelques explications. Je vous rappelle que d'après la définition, il s'agit de calculer la dérivée partielle \(\frac{\partial E_x}{\partial x}\), puis de la dérivée partielle \(\frac{\partial E_y}{\partial y}\) et de les sommer.
Pour calculer chaque dérivée partielle, je vais utiliser la fonction diff() de python, qui calcule la différence F[n+1] - F[n] pour un vecteur et j'applique la définition de la dérivée en posant diff(Fx)/hx, avec hx petit. Idem pour Fy. Vous avez sans doute reconnu la méthode d'Euler au premier ordre. La dérivation est donc un peu frustre mais conviendra bien à ce que nous voulons faire ici. Vous pourrez toujours l'améliorer en utilisant une RK4.
Je reviendrai plus loin sur le paramètre axis, du à l'usage de la fonction meshgrid() qui fait que Fx et Fy sont des matrices carrées.
Mais pourquoi [:-1,:] et [:,:-1] ? Bien que Ex et Ey soient des matrices carrées, je n'utiliserai que la première colonne pour Ex et la seconde pour Ey. Mais pourquoi ne prendre que les n-1 premiers termes de chaque colonne, ce que signifie :-1 ? Et bien parce que la fonction diff() appliquée à un vecteur de n termes ne retourne que n-1 termes, ce qui vu son travail, n'est pas surprenant... Pour rappel, elle calcule pour chaque n variant entre 0 et n-1, la différence F[n+1]-F[n] !
Je traite ici un champ 2D mais le calcul est directement transposable en 3D, voire sur des dimensions supérieures.
Le reste du code consiste à construire la grille de calcul nécessaire à l'affichage et au lancement des calculs du champ puis de sa divergence.
Je commence par définir les deux vecteurs des coordonnées x et y du domaine de calcul. Ces vecteurs sont uniquement destinés à construire la grille de calcul. Vous pouvez bien sur changer les bornes du domaine et faire varier les pas de calcul hx et hy. Rien n'oblige à ce qu'ils soient égaux, mais par contre il est nécessaire qu'ils soient petits devant l'étendue du domaine: souvenez-vous qu'ils servent au calcul des dérivées partielles...
hx = 0.01
hy = 0.01
X = arange(-2.,2.,hx)
Y = arange(-2.,2.,hy)
Puis je construis la grille de calcul à l'aide de la fonction meshgrid():
x,y = meshgrid(X,Y)
Cette fonction mérite quelques commentaires, tant son usage est parfois déroutant. Son objet est la construction d'une matrice carrée, une grille, à partir des vecteurs de coordonnées X et Y que nous avons défini ci-dessus. meshgrid() retourne curieusement deux matrices carrées que j'ai appelé x et y de dimension X*X et Y*Y. et encore plus curieusement on a x[i,j] = X[j] et y[i,j]= Y[i] pour tout i et j.
Quand dans ma définition de la divergence, j'ai écrit diff(Fx,axis=1), le paramètre axis = 1 signifie que je calcule sur la colonne 1, c'est à dire à la seconde car python commence les indices à 0, de la matrice Fx. De même, dans diff(Fy,axis=0), je calcule sur la colonne 0 de la matrice Fy.
Le lancement des calculs est trivial:
Ex, Ey = Champ(x,y)
divE = Divergence(Ex,Ey)
Les instructions de tracé n'appellent pas de commentaire. Et voilà le résultat obtenu :
Sans surprise, nous constatons que la divergence est nulle.
Vous pouvez changer la définition du champ, pour calculer la divergence d'un champ non nulle, par exemple le champ \( \vec{E} = \frac{x^2}{2} \vec{e_x} + xy \vec{e_y}\). Nous obtenons la divergence suivante:
Un rapide calcul manuel vous permettra de vérifier le résultat.
Nous avons vu plus haut que l'on pouvait définir le rotationnel (curl en anglais) d'un champ de vecteur \( \vec{A}\) comme \( \vec{rot} \: \vec{A} = \vec{\nabla} \wedge \vec{A} \). Comme tous les produits vectoriels, le rotationnel est un pseudo-vecteur, ou vecteur axial. Cela signifie que son orientation dépend du sens du triède choisi pour définir la base vectorielle. Attention donc aux histoires de symétrie...
On peut aussi le définir de façon intrinsèque en faisant appel à la notion de circulation d'un vecteur sur un contour fermé, notion que vous avez abordé en électrostatique. Imaginons un élement de surface \( \Delta S \), limité par un contour C. On définit le rotationnel par \( \vec{rot} \: \vec{A} = \lim_{\Delta S \rightarrow 0} \frac{1}{\Delta S}(\vec{a_n} \oint\! \vec{A}.\vec{dl}) \), en d'autres termes, la quantité vectorielle de direction normale à \( \Delta S \), orientée telle que la circulation de \( \vec{A}\) sur le contour C soit maximale.
Les coordonnées cartésiennes du rotationnel de \( \vec{A}\) se calculent aisément à partir du déterminant suivant :
\( \vec{rot} \: \vec{A} = \vec{\nabla} \wedge \vec{A} = \left| \begin{array}{cc} \vec{e_x}&\vec{e_y}&\vec{e_z}\\ \frac{\partial}{\partial x}&\frac{\partial}{\partial y}&\frac{\partial}{\partial z}\\ A_x&A_y&A_z \end{array} \right| \)
Il en est de même pour les coordonnées cylindriques :
\( \vec{rot} \: \vec{A} = \vec{\nabla} \wedge \vec{A} = \frac{1}{r}\left| \begin{array}{cc} \vec{e_r}&\vec{ r e_{\theta}}&\vec{e_z}\\ \frac{\partial}{\partial r}&\frac{\partial}{\partial \theta}&\frac{\partial}{\partial z}\\ A_r&A_{\theta}&A_z \end{array} \right| \)
Enfin, en coordonnées sphériques, selon le même principe :
\( \vec{rot} \: \vec{A} = \vec{\nabla} \wedge \vec{A} = \frac{1}{r^2 \sin(\theta)}\left| \begin{array}{cc} \vec{e_r}&\vec{ r e_{\theta}}&\vec{e_{\Phi}}\\ \frac{\partial}{\partial r}&\frac{\partial}{\partial \theta}&\frac{\partial}{\partial \Phi}\\ A_r&r A_{\theta}&r \sin(\theta) A_{\Phi} \end{array} \right| \)
Essayons peut-être plus intuitivement : il est sans doute question de rotation lorsqu'on parle de rotationnel ! Effectivement, cet opérateur nous vient de la mécanique des fluides, où il désigne l'axe de rotation d'un toubillon du champ de vitesse d'un fluide. Donc lorsque vous verrez sur une carte de champ vectoriel un enroulement des lignes de champ autour d'un point, dites vous que le rotationnel de ce champ ne doit pas être nul... Vous verrez plus loin sur un exemple ce que je veux dire.
Un théorème important de l'électromagnétisme lie la circulation d'un champ \( \vec{A}\) autour d'un contour C fermant un élément de surface \( \vec{dS}\). Il s'agit du théorème de Stokes-Ampère qui établit une égalité entre une intégrale de circulation et une intégrale de surface (souvent utile):
\( \oint\! \vec{A}.\vec{dl} = \iint\limits_S \vec{\nabla} \wedge \vec{A}.\vec{dS} \)
Le script python OPDChampRotationnel.py est très semblable au script de calcul de la divergence, à la différence essentielle de la routine de calcul du rotationnel ! Comme pour la divergence, je travaille sur un champ 2D, mais le principe est directement transposable à la 3D.
Le rotationnel est calculé par la routine Rotationnel :
def Rotationnel(x,y,Fx,Fy):
dFx = (diff(Fx,axis=0)/hy)[:,:-1]
dFy = (diff(Fy,axis=1)/hx)[:-1,:]
return dFy - dFx
Vous constatez que j'applique la formule de calcul à l'état brut, en utilisant la fonction diff() comme pour le calcul de la divergence. Toujours la méthode d'Euler, et toujours les mêmes contraintes dues à diff() et à l'expression matricielle de Fx et Fy.
Le reste du code est connu. J'ai adapté les valeurs des bornes de arange() en fonction des besoins. Il reste maintenant à tester notre code sur des champs dont le rotationnel est facile à calculer manuellement.
Commençons par tracer le champ \( \vec{E} = -\sin(y) \vec{e_x} + \cos(x) \vec{e_y}\), dont le rotationnel est non nul. Vour pourrez vous en assurer en faisant le calcul manuel. Voici l'aspect du champ, dont vous noterez les "vortex" répartis périodiquement sur le domaine. Certains vortex "tournent" dans le sens horaire et d'autres dans le sens antihoraire. Ce qui devrait se traduire dans la valeur du rotationnel.
Voyons si l'hypothèse est vraie. Le rotationnel obtenu à l'aide du script OPDChampRotationnel.py est :
Vous pouvez constater que les vortex de sens horaire possèdent un rotationnel négatif et les autres un rotationnel positif. Notez aussi que la valeur du rotationnel augmente lorsqu'on se rapproche du point autour duquel le vortex semble tourner. Vous trouverez sans doute une source d'inspiration de ce phénomène en mécanique des fluides !
Passons maintenant à un champ de rotationnel nul, par exemple le champ \( \vec{E} = -x \vec{e_x} + y \vec{e_y}\). Il ressemble à ça :
Sans surprise, nous observons un rotationnel constant et nul :
Il existe deux types de Laplacien, un Laplacien scalaire qui s'applique à une fonction supposée continue et au moins deux fois dérivable sur son domaine, et un Laplacien vectoriel qui s'applique à un champ de vecteurs. Je me limiterai ici au Laplacien scalaire, sachant que tous les résulats présentés ici sont transposables au Laplacien vectoriel.
Le Laplacien d'une fonction f(x,y,z) s'exprime ainsi en coordonnées cartésiennes: \( \Delta f(x,y,z) = \frac{\partial ^2 f(x,y,z)}{\partial x ^2} + \frac{\partial ^2 f(x,y,z)}{\partial y ^2} + \frac{\partial ^2 f(x,y,z)}{\partial z ^2}\)
On peut aussi, et de manière plus élégante le définir par \( \Delta f(x,y,z) = div(\vec{grad} \: f)\). Pour un Laplacien vectoriel, cette définition devient \( \Delta \vec{f}(x,y,z) = \vec{grad}(div \: \vec{f}(x,y,z)) - \vec{rot}(\vec{rot} \: \vec{f}(x,y,z) ) \)
A noter : une fonction de Laplacien nul est appelée fonction harmonique.
En coordonnées cylindriques, avec la base cylindrique définie plus haut, on a :
\( \Delta \: f(r,\theta,z) = \frac{1}{r} \frac{\partial}{\partial r}( r \: \frac{\partial f(r,\theta,z)}{\partial r}) + \frac{1}{r^2}\frac{\partial ^2 f(r,\theta,z)}{\partial \theta ^2} + \frac{\partial ^2 f(r,\theta,z)}{\partial z^2} \).
En coordonnées sphériques, avec la base sphérique définie plus haut, on a:
\( \Delta \: f(r,\theta,\Phi) = \frac{1}{r^2} \frac{\partial}{\partial r}( r^2 \: \frac{\partial f(r,\theta,\Phi)}{\partial r}) + \frac{1}{r^2\sin{\theta}} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial f(r,\theta,\Phi)}{\partial \theta}) + \frac{1}{r^2 \sin^2{\theta}}\frac{\partial^2 f(r,\theta,\Phi)}{\partial \Phi^2} \).
Le Laplacien est un opérateur de moyennage. Il indique la valeur de la fonction en un point donné par rapport aux points de son voisinage. Il mesure l'anomalie locale de la fonction en un point p donné. C'est le sens de la somme des dérivées partielles dans toutes les directions depuis le point de calcul. C'est aussi le sens du calcul des différences finies.
Revenons à l'équation de Laplace que nous avons déjà étudiée \( \Delta V = 0\) où V est une fonction potentiel V(x,y). Cette équation signifie que la valeur de V en un point p(x,y) est égale à la moyenne des valeurs de V dans le voisinnage de p. Soit V est linéaire autour de p, soit les variations dans une direction compensent celles des autres directions. Dans un domaine où l'équation de Laplace est vraie, il ne peut y avoir d'extremum local.
Géométriquement, le Laplacien mesure la courbure moyenne d'une fonction en un point donné.
J'ai testé deux possibilités de calculer un Laplacien : classiquement par la méthode des différences finies, comme nous l'avons utilisé dans la page sur l'équation de Laplace, et par la méthode directe, en employant la fonction diff() qui permet de calculer des dérivées secondes.
Je me suis limité au cas du Laplacien scalaire d'une fonction 2D, sachant que les algorithmes utilisés peuvent être très facilement adaptés en 3D et pour calculer un Laplacien vectoriel.
Le script python OPDLaplacienScalaireDF.py implémente la méthode des différences finies, que l'on retrouve dans la fonction LaplacienScalaire ci dessous :
def LaplacienScalaire(V,dx,dy):
lap = (V[0:-2,1:-1] + V[2:,1:-1] + V[1:-1,0:-2] + V[1:-1,2:] - 4*V[1:-1,1:-1])/(dx*dy)
return lap
Vous retrouvez le code très classique de la méthode des DF centrées, comme dans la résolution de l'équation de Laplace.
Le script python OPDLaplacienScalaireDiff.py met en oeuvre la fonction diff() pour dériver au second ordre, comme l'indique le paramètre n=2 dans l'appel de la fonction. On retrouve les mêmes artifices techniques que précédement pour l'usage de diff() :
def LaplacienScalaire(x,y,F):
dFx = (diff(F,n=2,axis=1)/hx**2)[:-2,:]
dFy = (diff(F,n=2,axis=0)/hy**2)[:,:-2]
return dFx + dFy
Pour obtenir une précision acceptable, on veillera à choisir des valeurs assez petites pour hx et hy. Typiquement, je fixe hx = hy = 0.001. Le temps de calcul varie en fonction de votre machine et peut prendre plusieurs secondes.
Voyons ce que cela donne sur une fonction définie par :
def Fonction(x,y):
return sqrt(x**2 + y**2)
Sa surface représentative est la suivante :
Vous noterez le puit au voisinage du point (0,0), notre fonction n'étant pas définie en ce point comme vous le savez et la discrétisation du tracé évitant les pointes vertigineuses...
Le tracé du Laplacien de la fonction montre néanmoins l'anomalie locale autour de (0,0) :
Les deux méthodes donnent le même résultat, pour des temps de calcul assez semblables.
Quelques relations liant les opérateurs différentiels de la physique, qu'il est utile de connaitre pour simplifier vos calculs, particulièrement en électromagnétisme et en méca flotte :
Les scripts Python étudiés dans cette page sont disponibles dans le package OPDPython.zip :
Les opérateurs différentiels sont des outils incontournables pour le physicien. Leur maitrise est indispensable, non seulement pendant vos études, mais aussi dans la vie quotidienne du laboratoire. La pratique de leur calcul numérique est indispensable à ceux qui se destinent au labo, à l'industrie ou à la recherche. Je vous ai donné des exemples en python, mais le principe des algorithmes reste identique en scilab, matlab, maple ou C. D'ailleurs, MatLab et Maple possèdent des fonctions internes de calcul de la plupart des opérateurs. Il est aussi possible en Maple de faire du calcul symbolique sur les opérateurs. Cette possibilité existe en python en utilisant simpy, piste que j'explorerai dans une prochaine page.
Contenu et design par Dominique Lefebvre - tangenteX.com août 2016 -- Vous pouvez me joindre par ou sur PhysiqueX
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.