calculs de divergence et rotationnel avec Maxima

Le 30 septembre 2017 dernier nous apprenions la mort , à seulement 51 ans, du mathématicien Russe Vladimir Voevodsky. Médaillé Fields en 2002, pour sa démonstration de la célèbre conjecture de Milnor, une partie de son travail s'est révélée fausse près de 10 ans plus tard! Prenant conscience de la complexité croissante des preuves mathématiques il avait réorienté ses recherches vers le domaine des assistants de preuve formelle ,comme COQ , et voulait fournir aux mathématiciens des outils pour écrire et vérifier plus rigoureusement des preuves. Dans la croyance populaire on imagine qu'il suffirait de cliquer sur un bouton pour vérifier ou même démontrer un énoncé, évidement la réalité est bien différente et une preuve assistée par un logiciel  se révèle souvent très technique. On peut s'en rendre compte en s'intéressant au calcul formel avec des logiciels populaires comme Maxima. Démontrer les formules qui expriment la divergence où le rotationnel d'une fonction vectorielle dans un système de coordonnées non-cartésien   demande beaucoup d'habileté mathématique. C'est un bon exemple pour comprendre et réfléchir à ce que représente un calcul assisté par ordinateur ...

$$ \begin{align*}
{\rm div}({\bf A})
&=\partial_x {\bf A}_x+\partial_y{\bf A}_y+\partial_z{\bf A}_z
\\&={\frac {1}{r}} {\partial_r (rA_{r})}+{\frac {1}{r}}{\partial_\theta  A_{\theta }}+{\partial_z A_{z}}\\
\begin{array}{c}{\bf rot}({\bf A})\\~ \\ ~\end{array}&\begin{array}{c}= \\~ \\ ~\end{array}
\begin{pmatrix} {\partial_y \mathrm{A}_z } - {\partial_z \mathrm{A}_y} \\
{\partial_z \mathrm{A}_x } - {\partial_x \mathrm{A}_z }\\
{\partial_x \mathrm{A}_y } - {\partial_y \mathrm{A}_x } \end{pmatrix}
\begin{array}{c}= \\~ \\ ~\end{array}
 \begin{array}{r}
\left(\frac{1}{r}{\partial_\theta \mathrm{A}_z} - {\partial_z \mathrm{A}_\theta}\right) \mathbf{e_r} \\
+ \left({\partial_z \mathrm{A}_r} - {\partial_r \mathrm{A}_z}\right)\mathbf{e_\theta} \\
+\frac{1}{r}\left({\partial_r}(r \mathrm{A}_\theta) - {\partial_\theta \mathrm{A}_r}\right) \mathbf{e_z}
\end{array}
\end{align*}$$

Coordonnées cylindriques 

Pour illustrer ce qu'on peut faire avec maxima je vais prendre le cas le plus simple de changement de coordonnées qu'on voit lors d'études universitaires :  celui des coordonnées cylindriques :

\[x=r\, \mathrm{cos}\left( \theta\right) ,y=r\, \mathrm{sin}\left( \theta\right) ,z=z\]
Du point de vu mathématique ce changement de coordonnées est un difféomorphisme $\phi(r,\theta,z)=(\phi_x,\phi_y,\phi_z)$ dont le Jacobien $J\phi(r,\theta,z)$ fait apparaitre une base orthonormée :
$$ {\bf e}_r =\begin{pmatrix}\cos(\theta)\\\sin(\theta)\\0\end{pmatrix}
{\bf e}_\theta =\begin{pmatrix}-\sin(\theta)\\\cos(\theta)\\0\end{pmatrix}
u_z {\bf e}_z=\begin{pmatrix}0\\0\\1\end{pmatrix}
$$
On peut ensuite exprimer un vecteur ${\bf u}$ de coordonnées cartésiennes $(u_x,u_y,u_z)$ dans cette nouvelle base ${\bf e}_r ,{\bf e}_\theta {\bf e}_z$:
$$(1)~~~
{\bf u}=
\begin{pmatrix}u_x\\u_y\\u_z\end{pmatrix}
=
\underbrace{(\mathrm{sin}\left( \theta\right) {{u}_{y}} +\mathrm{cos}\left( \theta\right)  {{u}_{x}})}_{=u_r} {\bf e}_r
+\underbrace{(\mathrm{cos}\left( \theta\right){{u}_{y}} -\mathrm{sin}\left( \theta\right) {{u}_{x}})}_{=u_\theta}{\bf e}_\theta
+u_z {\bf e}_z$$

cette base permet d'exprimer facilement les quantités physiques en présence quand la configuration est invariante suivant certaines transformations (rotation/translation suivant l'axe des  z ici). Pour exploiter ces informations on a donc  besoin d'exprimer  les opérateurs de divergence et de rotationnel en fonction uniquement des composantes dans cette base et des dérivées par rapport aux coordonnées cylindriques :

$$ {\rm div}({\bf u})=F(u_r,u_\theta,u_z,\partial_r,\partial_\theta,\partial_z,r,\theta,z)$$

Pour cela on a déjà trouvé une correspondance entre les composantes cartésiennes $u_x,u_y,u_z$ et cylindriques $u_r,u_\theta,u_z$ mais il faut en plus trouver une correspondance entre les dérivations $\partial_x,\partial_y,\partial_z$ et $\partial_r,\partial_\theta,\partial_z$. cela s'obtient en utilisant la dérivation composée
$$ \begin{align*}
(2)~~~
\partial_t u(\phi(r,\theta,z)&=\partial_t \phi_x(r,\theta,z)\times \partial_x u(\phi(r,\theta,z))
\\&\phantom{=}+\partial_t \phi_y(r,\theta,z)\times \partial_x u(\phi(r,\theta,z))
\\&\phantom{=}+\partial_t \phi_z(r,\theta,z)\times \partial_x u(\phi(r,\theta,z))
,~~~t=r,\theta,z
\end{align*}
$$
ce qui revient à écrire en terme de matrice Jacobienne que  $J(u\circ\phi)(r,\theta,z) =Ju(\phi(r,\theta,z)\times J\phi(r,\theta,z)$, c'est beaucoup moins parlant et  les physiciens écrivent carrément par abus de notation:
$$ \partial_t u(x,y,z)=\partial_t x\times \partial_x u(r,\theta,z)
+\partial_t y\times \partial_x u(r,\theta,z)
+\partial_t z\times \partial_x u(r,\theta,z)
$$
Il faut donc appliquer la formule (2)  aux composantes écrites en coordonnées cylindriques dans (1)  puis les combiner pour essayer d'obtenir  les formules en coordonnées Cartésiennes!

Utilisation de Maxima

C'est ici que je vais utiliser maxima, ce moteur de calcul formel possède un langage de programmation assez simple avec lequel on va écrire un programme appliquant la forme de dérivation composée (2) à une fonction quelconque u :

load(eigen)$
derive_part(u,phi):=block([L,newvars,eq,var,F],
    newvars:[],    
    for k:1 thru length(phi) do (
        newvars:append(newvars,[rhs(phi[k])])
        ),
    newvars:listofvars(newvars),
    L:makelist(0,i,1,length(newvars)),
    for j:1 thru length(newvars) do (
        add_diff:true,
        F:0,
        for k: 1 thru length(phi) do(
            eq:rhs(phi[k]),
            if eq = newvars[j] then add_diff:false,
            var:lhs(phi[k]),
            F:F+diff(eq,newvars[j])*diff(u,var)   
        ),
        if add_diff then F:F+diff(u,newvars[j]),
        L[j]:F    
    ),
    L:covect(L)
)$

la seule subtilité dans la fonction derive_part  est qu'elle va s'appliquer à des expressions u qui peuvent elles même contenir les nouvelles variables  (en fait  $\theta$ pour les coordonnées cylindriques). On teste facilement la fonction avec les coordonnées cylindriques :

phi:[x=r*cos(theta),y=r*sin(theta),z=Z];
du:derive_part(u(x,y,z),phi)$
'diff(u(x,y,z),r)=du[1];
'diff(u(x,y,z),theta)=du[2];
'diff(u(x,y,z),Z)=du[3];

pour ceux qui ne connaissent pas maxima on notera que :

• l'apostrophe "  ' "  devant le diff fait que cet opérateur est appliqué symboliquement (sans effectuer le calcul) 
• le "$\$ $ " empêche d'afficher le résultat d'un calcul (les autres lignes finissent par ";")
• le "=" ne correspond pas à une affectation (qui est faite avec ":") ce qui permet d'écrire la signification des termes stockés dans la variable du

on va donc appliquer cette fonction aux composantes de u en coordonnées cylindriques :

'u[r]=u[r]:cos(theta)*u[x](x,y,z)+sin(theta)*u[y](x,y,z);
'u[theta]=u[theta]:-sin(theta)*u[x](x,y,z)+cos(theta)*u[y](x,y,z);
'u[Z]=u[Z]:u[z](x,y,z);
du_r:derive_part(u[r],phi)$
du_theta:derive_part(u[theta],phi)$
du_Z:derive_part(u[Z],phi)$

avec wxmaxima on a un affichage des formules assez agréable, car utilisant LaTeX , et qui reconnaît les lettre grecques !

On peut maintenant vérifier la formule de la divergence  en calculant :

div(a):='diff(a[1],x)+'diff(a[2],y)+'diff(a[3],z)$;
'div('u)=divu_cart:div([u[x](x,y,z),u[y](x,y,z),u[z](x,y,z)]);
'diff('u[r],r)=trigsimp(du_r[1]);
'diff('u[theta],theta)=trigsimp(du_theta[2]);
'diff('u[r],r)+(1/r)*'diff('u[theta],theta)+(1/r)*'u[r]+'diff('u[z],z)
=trigsimp(du_r[1]+(1/r)*du_theta[2]+(1/r)*u[r]+du_Z[3]);

ce qui donne

en comparant les expressions de $\partial_r u_r$ et de $\partial_\theta u_\theta$ on voit bien qu'il faut diviser ce dernier par r pour faire s'annuler les termes $\partial_x u_y$ et $\partial_y u_x$ en additionnant les 2 équations. On va bien retrouver $\partial_x u_x+\partial_yu_y$ et un terme supplémentaire qu'on reconnaît bien comme étant $-u_r$ qui faut retrancher pour obtenir la formule de la divergence (en ajoutant en plus $\partial_z u_z$).

Il faut faire le même type de calculs pour chaque composante du rotationnel en cherchant comment combiner successivement $\partial_r u_\theta$ et $\partial_\theta u_r$,  $\partial_\theta u_z$ et $\partial_z u_\theta$ et $\partial_\theta u_r$ puis $\partial_r u_z$ et $\partial_z u_r$. Ça n'est pas plus compliqué mais c'est techniquement difficile :




Tous ces calculs peuvent paraître bien compliqués, tout autant que de le faire "à la main", cependant il suffit de changer la variable phi  en prenant maintenant les coordonnées sphérique pour faire les calculs de la divergence et du rotationnel dans ces nouvelles coordonnées, calcul qui pour le coup est vraiment très difficile à faire "à la main". Si vous voulez essayer voici le fichier maxima au format *.mac, celui pour wxmaxima  *.wxmx, et le résultat au format *.pdf

conclusion

Difficile de résumer le travail du mathématicien sans faire des raccourcis simplificateurs mais on peut dire que l'apport initial de Voevodsky en géométrique algébrique a permis de démontrer la célèbre conjecture de Milnor  ce qui lui valu d'être récompensé par la  médaille Fields en 2002. Le travail de Voevodsky était très complexe, tellement complexe que peu de relecteurs étaient capables de l'évaluer et après sa médaille Fields des doutes sont apparus sur la validité de son travail. Ces doutes ont définitivement été confirmés en 2013: le théorème principal primé par la médaille Fields comportait des erreurs importantes que personne n'avait détecté!!!!  Mais bien avant cela Voevodsky s'était rendu compte d'une partie du problème et après une période de profonde déprime il avait décidé de se consacrer aux questions  des assistants de preuve formelles. Depuis il travaillait à la conception de programme comme COQ capable d'effectuer des preuves logiques ou de les réécrire pour les rendre humainement lisible (pour ceux intéressés lire cette  interview en Français ,  ce billet du mathématicien Johan Baez ou cet article de la revue nautilus en anglais ). Cela montre bien au final que Voevodsky est vraiment un grand mathématicien, au point d'admettre la place de l'erreur dans le processus mathématique s'intégrant dans une branche des mathématiques impliquant Hilbert , Russel, Godel ou Turing.

Cela pose aussi  la question de la place de l'informatique dans le calcul mathématique. Par exemple lors de la mise en évidence observationnelle des ondes gravitationnelles j'ai pu lire que le calcul formel a joué un rôle important dans la recherche fondamentale sur ce sujet. En effet le calcul des tenseurs de la relativité générale est techniquement extrêmement difficile, beaucoup de résultats ont été conjecturés sur des formes simplifiées des équations (linéarisation) et ont ensuite été amélioré en utilisant du calcul formel. On trouvera d'ailleurs dans maxima des packages comme itensor et ctensor pour le calcul des tenseurs.