Calcul Différentiel Élémentaire
*********************************

Version du Friday 21 May 2021 à 16H41 (le manuel n'est pas encore
complet).

Introduction pour...  une introduction !

Introduction
============

_Calcul différentiel elémentaire_ (CDE) a pour intention de
faciliter les calculs algébriques et différentiels élémentaires.
Notamment pour les champs de scalaires et vecteurs de l'espace
ordinaire de dimension 3, pris comme absolu et dans lequel des objets
se déplacent, le temps étant le paramètre de ce déplacement.

Mais si l'intention est bienveillante, sa concrétisation peut ne pas
l'être tant que cela.  CDE introduit des catégories dans les
variables de Maxima (les vecteurs et les dyades) et des fonctions
spéciales pour manipuler les éléments de ces catégories.

De la même façon le calcul différentiel est compliqué par
l'introduction des opérateurs de dérivation gradient, divergence et
rotationnel avec lesquelles sont formées des équations aux
dérivées partielles (EDP).

Comme on ne peut pas retirer de résultats quantitatifs à partir de
la contemplation de ces EDP , il faut bien en calculer sinon les
solutions du moins une approximation de celles-ci.  C'est à dire
passer par des algorithmes de calcul numérique.

À cet effet, CDE interface FreeFem++.  Ce qui engendre encore des
fonctions dont il est important de mémoriser le nom.

1 Scalaires
***********

La plupart des opérations qu'on peut souhaiter effectuer sur des
scalaires est assurée par Maxima sans grandes transformations
Maxima.  Certaines ne le sont cependant pas ou alors
nécessitent l'emploi de plusieurs commandes qui gagnent à être
regroupées.

C'est celles-ci qui sont décrites ici, par catégories.

1.1 Résolution algébrique
===========================

1.1.1 ElVars
------------

la fonction de résolution de Maxima est 'solve' dans
Functions and Variables for Equations (Maxima).  Elle présente
l'inconvénient de ne rien retourner en cas d'échec de la
résolution.  D'où l'intérêt de

*ElVars(eq,x)* qui résoud les systèmes d'équations eq (qui doit
être une liste) par rapport aux variables de x (qui doit aussi être
une liste).

Et le fait bêtement : si x=[x1,x2,...]  et eq=[eq1,eq2,...], ElVars
initialise un liste de solutions à sol=[], puis cherche à résoudre
eq1 par rapport à x1 (en utilisant 'solve').

S'il y a une solution, la liste d'équation de devient
eq=[eq2,eq3,...], la liste de solution devient sol=[x1=s1] où s1 est
le résultat obtenu qui dépend le plus souvent des variables
x2,x3,...  et la liste de variable devient x=[x2,x3,...]

S'il n'y a pas de solution, Elvars passe à eq2 puis eq3 et les autres
jusqu'au moment où une des équations peut être résolue par rapport
à x1.

Si cela n'arrive pas, la liste d'équations reste eq=[eq1,eq2,...],
celle de solutions sol:[] mais la liste des variables devient
x=[x2,x3,...]  (x1 est oublié).

Si cela arrive, disons pour eq3 alors la liste d'équations devient
eq=[eq1,eq2,eq4,...], la liste de solution sol=[x1=s3] où s3 peut
dépendre en général de x2,x3,...  et la liste de variable
x=[x2,x3,...].

Une fois cela fait ElVars continue d'effectuer les résolutions comme
décrit pour les variables de x=[x2,x3,...].  Mais en plus, à chaque
résolution effectuée avec succès la variable résolue est
remplacée par sa valeur dans sol.

Quand il ne reste plus de variables de x, le processus se termine et
ElVars renvoie [sol,eq]

Voici quelques exemples de fonctionnement de *ElVars*

Fonctionnement analogue à *solve*
     %%%; nettoie()$
     ElVars([x+y,x-y=1],[x,y])@
     solve([x+y,x-y=1],[x,y])@
     %%%;

Fonctionnement différent
     %%%; nettoie()$
     ElVars([x+y],[x,y])@
     solve([x+y],[x,y])@
     %%%;
solve introduit une paramétrisation avec une variable intermédiaire
; pas ElVars.

Encore différent
     %%%; nettoie()$
     ElVars([x+y],[y,x])@
     solve([x+y],[y,x])@
     %%%;
Pour solve, c'est là même chose, pas pour ElVars : l'ordre des
variables par rapport auxquelle se fait la résolution est important.

Fonctionnement aberrant
     %%%; nettoie()$
     ElVars([x^2+y],[x])@
     solve([x^2+y],[x])@
     %%%;
Là ElVars se trompe.  Il est à manier avec prudence dans les cas
non-linéaires.

1.1.2 canonise
--------------

*canonise(expr,x)*

Cherche des parties de 'expr' qui sont une forme linéaire ou
quadratique de v, soit
     b*v+c
ou
     a*v^2+b*v+c
et les met sous forme canonique, soit
     canonise(b*v+c,v) -> b*(v+c/b)
ou
     canonise(a*v^2+b*v+c,v) -> a*((v+b/(2*a))^2+(4*a*c-b^2)/(4*a^2))

*canonise(expr,x,y,...)*

Effectue l'opération de canonisation successivement avec les
variables x,y,...

Par exemple
     %%%; nettoie()$
     canonise(x^2+x*y+y^2+y,x,y)@
     canonise(x^2+x*y+y^2+y,y,x)@
     %%%;

1.2 Dérivation
===============

*derive(expr,ld)* où _expr_ est une expression quelconque dépendant
de x,y,...  et _ld_ une liste de la forme
     [[x,dx],[y,dy],...]
calcule la forme linéaire tangente à _expr_ en x,y,...  des
variations dx,dy,...  de x,y,...

_derive_ reconnait la linéarité des opérateurs linéaires standards
(_integrate_ par exemple) et de ceux que CDE introduit (_grad_ par
exemple).

2 Vecteurs et Dyades
********************

2.1 Déclaration et produits
============================

Les vecteurs sont les éléments de
     liste_des_vecteurs@
les dyades(1) de
     liste_des_dyades@
Les seules fonctions qui les reconnaissent comme tels sont les 'tav',
'xtav', 'rtav' qui effectuent des simplifications sur les expressions
comportant des vecteurs et/ou des dyades.

Déclaration de vecteurs et dyades
----------------------------------

Pour déclarer comme vecteur(s) un ou plusieurs symboles
     %%%; nettoie()$ addv(x,y,z)@ %%%;
c'est à dire ajouter x,y,z à liste_des_vecteurs.

Idem pour des dyades
     %%%; nettoie()$ addy(X,Y,Z)@ %%%;
ce qui ajoute X,Y,Z à liste_des_dyades

Pour enlever des vecteur de liste_des_vecteurs et une dyade de
liste_des_dyades
S'ils y sont il sont enlevés, s'ils n'y sont pas rien n'est fait.

Produits scalaire, vectoriel et mixte de vecteurs
-------------------------------------------------

Ce sont les fonctions 'ps', 'pv', 'pm'.

     %%%; nettoie()$ addv(x,y,z)$
     ps(x,y)@
     pv(x,y)@
     pm(x,y,z)@
     %%%;

Noter que les 'ps', 'pv', 'pm' ne sont en elles-même associées à
aucune opération, ce sont les 'tav's qui leur confèrent cette
signification.

   ---------- Footnotes ----------

   (1) Gibbs avait appelé dyadics les tenseurs d'ordre 2 dans
l'espace à 3 dimensions.

2.2 Calculs algébriques avec des vecteurs
==========================================

Fonction de simplification d'expressions vectorielles
-----------------------------------------------------

Il y a une seule fonction de simplification qui est 'tav'.  'xtav'
appelle deux fois 'tav' et 'rtav' fait un 'tav' puis ordonne un
expression vectorielle.

2.2.1 Simplification de base
----------------------------

Le double produit vectoriel est développé par 'tav'
     %%%; nettoie()$ addv(x,y,z)$ xyz:pv(pv(x,y),z)@ xtav(xyz)@ %%%;

Le produit scalaire avec un produit vectoriel est reconnu comme
produit mixte
     %%%; nettoie()$ addv(x,y,z)$ xyz:ps(pv(x,y),z)@ xtav(xyz)@ %%%;

Les linéarisations sont effectuées
     %%%; nettoie()$
     addv(x,y,z)@
     xtav(ps(x+mu*y,lambda*z))@
     xtav(pv(x+mu*y,lambda*z))@
     xtav(pm(x+mu*y,lambda*z,pv(x,y)))@
     %%%;

La base réciproque de la base réciproque est la base de départ à
un facteur près qui se trouve être le produit mixte des vecteurs de
cette base de départ
     %%%; nettoie()$
     addv(x,y,z)$
     [X,Y,Z]:[pv(y,z),pv(z,x),pv(x,y)]@
     xtav([pv(Y,Z),pv(Z,X),pv(X,Y)])@
     %%%;

2.2.2 Simplification avec Regles
--------------------------------

On peut demander ajouter des prescriptions pour la simplification.
Cela ce fait en les rangeant dans la variable globale 'Regles' sous la
forme 'Regle:[[expression,la valeur de l'expression],...]'.  Par
exemple
     %%%; nettoie()$ addv(x,y)@
     Regles:[[ps(x,y),cos(theta)]]@
     expr:ps(x,lambda*y)@
     xtav(expr)@
     %%%;

2.3 Calculs algébriques avec des dyades
========================================

Les dyades représent les endomorphismes linéaires agissant sur les
vecteurs pour donner d'autres vecteurs.  Et elles les représentent de
manière à s'insérer dans les expressions de calcul symbolique avec
le moins de frais syntaxiques possible.

2.3.1 Produits licites et illicites
-----------------------------------

La règle suivante est adoptée
     %%%; nettoie()$ addv(a)$ addy(U)$
     U*a@
     %%%;
Le produit 'U*a' d'une dyade par un vecteur est un vecteur.  On garde
la commutativité du produit ce qui fait que le produit 'a*U' d'un
vecteur par une dyade est ce même vecteur 'U*a'.

Si de plus il y a un scalaire dans le produit
     %%%; nettoie()$ addv(a)$ addy(U)$
     k*U*a@
     %%%;
c'est le produit de la dyade 'U' par le vecteur 'k*a' ou encore le
produit de la dyade 'k*U' par le vecteur 'a' ou encore celui du
scalaire 'k' par le vecteur 'U*a'.  Mais toutes ces interprétations
correspondent au même vecteur 'k*U*a'.

Le produit de deux dyades est un scalaire
     %%%; nettoie()$ addy(U,V)$
     U*V %%%;
ce scalaire 'U*V' est le produit contracté des deux dyades.  Et pas
la dyade qui correspondrait à la compostion des deux endomorphismes
représentés par 'U' et 'V'.  La raison de choix vient de ce qu'on
souhaite conserver la commutativité du produit, laquelle n'est pas
vraie pour la composition d'endomorphisme.

Un produit de deux dyades et d'un vecteur est un vecteur
     %%%; nettoie()$ addy(U,V)$ addv(a)$
     U*V*a %%%;
c'est le produit du scalaire 'U*V' et du vecteur 'a'.

Et sinon tout produit comportant plus de deux dyades ou plus de deux
vecteurs forme une expression qui sera incorrectement interprétée
par les fonctions 'tav'.

2.3.2 Operateurs algébriques sur les dyades
--------------------------------------------

Une dyade particulière est l'identité
     %%%; nettoie()$
     Id@
     %%%;

Son produit (contracté) avec une dyade est la trace de la dyade
     %%%; nettoie()$ addy(U)$
     tav(U*Id)@
     %%%;
qui est un opérateur qu'on peut appliquer à une dyade.  On a par
exemple
     %%%; nettoie()$ addy(U)$
     tav(Id*Id)@
     tav(Trace(Id))@
     %%%;

La transposée d'une dyade est
     %%%; nettoie()$ addy(U)$
     Trans(U)@
     %%%;
Elle est défini de manière que
     %%%; nettoie()$ addy(U)$ addv(a,b)$
     ps(Trans(U)*a,b)=ps(a,U*b)@
     %%%;
comme le reconnait 'tav'
     %%%; nettoie()$ addy(U)$ addv(a,b)$
     tav(ps(Trans(U)*a,b)-ps(a,U*b))@
     %%%;

2.3.3 Produit dyadique
----------------------

On peut construire une dyade par le produit dyadique de deux vecteurs
comme
     %%%; nettoie()$ addv(a,b)$
     C:dya(a,b)@
     %%%;
qui est défini de manière que
     %%%; nettoie()$ addv(a,b,c)$
     C:dya(a,b)@
     tav(C*c)@
     %%%;

La trace d'un produit dyadique est le produit scalaire des deux
vecteurs
     %%%; nettoie()$ addv(a,b)$
     tav(Trace(dya(a,b)))@
     tav(Id*dya(a,b))@
     %%%;

Le produit contracté de deux produits dyadiques est
     %%%; nettoie()$ addv(a,b,c,d)$
     tav(dya(a,b)*dya(c,d))@
     %%%;

La transposée d'un produit dyadique est
     %%%; nettoie()$ addv(a,b)$
     tav(Trans(dya(a,b)))@
     %%%;

2.3.4 Vecteurs axiaux et rotations élémentaires
-------------------------------------------------

La 'rotation élémentaire' d'un vecteur est la dyade
     %%%; nettoie()$ addv(a)$
     LC(a)@
     %%%;
qui est telle que
     %%%; nettoie()$ addv(a,b)$
     tav(LC(a)*b)@
     %%%;
Le produit d'une dyade rotation élémentaire 'LC(a)' et d'un vecteur
'b' est le produit vectoriel de 'a' et 'b' ; et vice versa d'ailleurs
: le produit vectoriel est le produit de la dyade rotation
élémentaire par un vecteur.

Le 'vecteur axial' d'une dyade est
     %%%; nettoie()$ addy(U)$
     CL(U)@
     %%%;

Le vecteur axial d'une dyade 'rotation élémentaire de vecteur a' est
le vecteur 'a'
     %%%; nettoie()$
     addv(a,b)@
     expr:CL(LC(a))@
     xtav(expr)@
     %%%;

La 'rotation élémentaire' du 'vecteur axial' d'une dyade est
     %%%; nettoie()$ addy(U)$
     xtav(CL(U))@
     xtav(LC(CL(U)))@
     %%%;
la partie antisymétrique de la dyade.

2.4 Calculs différentiels
==========================

2.4.1 Opérateurs vectoriels
----------------------------

Le gradient s'obtient par
     %%%; nettoie()$ grad(phi)@ %%%;
où 'phi' est une expression scalaire (ici un atome) ; les rotationnel
et divergence par
     %%%; nettoie()$ addv(v)$ rot(v)@ div(v)@ %%%;
où 'v' est une expression vectorielle (ici un atome) ; le gradient
d'une expression vectorielle par
     %%%; nettoie()$ addv(v)$ Grad(v)@ %%%;
c'est une dyade.  La transposée du gradient d'une expression
vectorielle par
     %%%; nettoie()$ addv(v)$ tGrad(v)@ %%%;
mais c'est destiné à disparaître, mieux vaut écrire
     %%%; nettoie()$ addv(v)$ Trans(Grad(v))@ %%%;

'tav' connait les formules d'analyse vectorielle usuelles
     %%%; nettoie()$ addv(v,w)$
     rot(grad(phi))=tav(rot(grad(phi)))@
     div(rot(v))=tav(div(rot(v)))@
     grad(phi*psi)=tav(grad(phi*psi))@
     div(phi*v)=tav(div(phi*v))@
     rot(phi*v)=tav(rot(phi*v))@
     div(pv(v,w))=tav(div(pv(v,w)))@
     %%%;
puis
     %%%; nettoie()$ addv(v,w)$
     Trace(Grad(v))=tav(Trace(Grad(v)))@
     Div(dya(u,w))=tav(Div(dya(v,w)))@
     %%%;

Il y a aussi les expressions qu'on peut souhaiter ou non de
développer.  Pour elles on dispose du switch 'PrefDyade' qui est
'false' par défaut, d'où
     %%%; nettoie()$ addv(v,w)$
     rot(pv(v,w))=tav(rot(pv(v,w)))@
     grad(ps(v,w))=tav(grad(ps(v,w)))@
     pv(rot(v),w)=tav(pv(rot(v),w))@
     %%%;
par contre
     %%%; nettoie()$ addv(v,w)$ PrefDyade:true$
      rot(pv(v,w))=tav(rot(pv(v,w)))@
      grad(ps(v,w))=tav(grad(ps(v,w)))@
      pv(rot(v),w)=tav(pv(rot(v),w))@
     %%%;

2.4.2 Spécification du vecteur position
----------------------------------------

Le vecteur par rapport auquel opèrent les opérateurs vectoriels peut
être spécifié en ajoutant à 'Regles' les simplifications ad hoc.
Cela se peut se faire manuellement en remplissant 'Regles' où plus
simplement par la fonction *FaitReglesX(x,Cste,Regs)*.

Les arguments sont 'x' le symbole représentant la variable position
(qui doit être déclaré comme vecteur) ; 'Cste' une liste de
symboles qui représentent des quantités constantes ; 'Regs' une
liste de règles additionnelles.  See Simplifications avec Regles.

Basiquement
     %%%; nettoie()$ addv(y)$
     Regles@
     FaitReglesX(y,[],[])$
     Regles@
     %%%;

Cela permet, par exemple, de calculer %%%; nettoie()$ addv(x)$
FaitReglesX(x,[],[])$ div(grad(1/nrm(x)))=xtav(div(grad(1/nrm(x))))@
%%%; '1/nrm(x)' est une fonction harmonique !  Ce n'est pas un scoop.
Mais on peut aussi faire des calculs moins faciles comme
     %%%; nettoie()$ addv(x,d)$
     FaitReglesX(x,[],[])$
     div(grad(ps(d,x)/nrm(x)^3))=xtav(div(grad(ps(d,x)/nrm(x)^3)))@
     FaitReglesX(x,[d],[])$
     div(grad(ps(d,x)/nrm(x)^3))=xtav(div(grad(ps(d,x)/nrm(x)^3)))@
     %%%;
Si 'd' est constant, 'ps(d,x)/nrm(x)^3' est une fonction harmonique.
Dans la première ligne, on n'a pas indiqué que 'd' était constant,
ce qu'on a fait pour la seconde.

On trouvera des calculs plus élaborés dans Décomposition de Helmholtz d'un champ dipolaire

2.4.3 Coordonnées
------------------

Les coordonnées sont prises en compte via le mécanisme des règles,
comme dans Spécification du vecteur position.

2.4.3.1 Coordonnées cartésiennes
..................................

*GenCart()* génère les règles ad hoc pour les coordonnées
cartésiennes.
     %%%; nettoie()$
     GenCart();
     %%%;
Celles-ci sont qu'il existe une base de trois vecteurs orthonormés
'kx,ky,kx' ; que ces trois vecteurs sont constants ; et que la
position générale est décrite par les coordonnées cartésiennes
'x,y,z' par
     %%%; x*kx+y*ky+z*kz@ %%%;

On peut utiliser ces coordonnées pour faire des calculs algébriques
comme différentiels.  Lorsque les calculs sont différentiels il est
nécessaire de déclarer la dépendance des variables utilisée par le
'depends' de Maxima Functions and Variables for Differentiation (Maxima).
Comme ceci
     %%%; nettoie()$ addvV(v)$ Indicie([x,y,z],v)$
     depends([psi,vx,vy,vz],[x,y,z])@
     %%%;
On a alors accès aux expressions des opérateurs vectoriels en
coordonnées cartésiennes
     %%%; nettoie()$ addvV(v)$ Indicie([x,y,z],v)$ GenCart()$
     depends([psi,vx,vy,vz],[x,y,z])$
     vv:vx*kx+vy*ky+vz*kz@
     grad(psi)=tav(grad(psi))@
     div(vv)=tav(div(vv))@
     rot(vv)=tav(rot(vv))@
     Grad(vv)=tav(Grad(vv))@
     %%%;
Si on préfère une notation matricielle pour 'Grad(vv)', on peut en
disposer avec
     %%%; nettoie()$ addvV(v)$ Indicie([x,y,z],v)$ GenCart()$
     depends([vx,vy,vz],[x,y,z])$
     vv:vx*kx+vy*ky+vz*kz@
     Grad(vv)=tav(Grad(vv))@
     TensMat(tav(Grad(vv)),[kx,ky,kz])@
     %%%;

2.4.3.2 Coordonnées cylindriques
.................................

*GenCart()* génère les règles ad hoc pour les coordonnées
cylindriques.
     %%%; nettoie()$
     GenCyl();
     %%%;
Celles-ci sont qu'il existe une base de trois vecteurs orthonormés
'kr,kt,kx' ; que 'kx' est constant mais pas 'kt,kt' qui sont de la
forme
     %%%; nettoie()$
     kr= cos(theta)*kx+sin(theta)*ky@
     kt=-sin(theta)*kx+cos(theta)*ky@
     %%%;
où 'kx,ky,kz' sont les vecteurs de la base cartésienne sur laquelle
sont construite les coordonnées cylindiques ; et que la position
générale est décrite par les coordonnées cylindriques 'r,theta,z'
par
     %%%; r*kr+z*kz@ %%%;

On peut utiliser ces coordonnées pour faire des calculs algébriques
comme différentiels.  Lorsque les calculs sont différentiels il est
nécessaire de déclarer la dépendance des variables utilisée par le
'depends' de Maxima Functions and Variables for Differentiation (Maxima).
Comme ceci
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,z],v)$
     depends([psi,vr,vtheta,vz],[r,theta,z])@
     %%%;
On a alors accès aux expressions des opérateurs vectoriels en
coordonnées cylindriques
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,z],v)$ GenCyl()$
     depends([psi,vr,vtheta,vz],[r,theta,z])$
     vv:vr*kr+vtheta*kt+vz*kz@
     grad(psi)=rtav(grad(psi))@
     div(vv)=tav(div(vv))@
     rot(vv)=rtav(rot(vv))@
     Grad(vv)=tav(Grad(vv))@
     %%%;
Si on préfère une notation matricielle pour 'Grad(vv)', on peut en
disposer avec
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,z],v)$ GenCyl()$
     depends([vr,vtheta,vz],[r,theta,z])$
     vv:vr*kr+vtheta*kt+vz*kz@
     Grad(vv)=tav(Grad(vv))@
     TensMat(tav(Grad(vv)),[kr,kt,kz])@
     %%%;

2.4.3.3 Coordonnées sphériques
................................

*GenSph()* génère les règles ad hoc pour les coordonnées
sphériques.
     %%%; nettoie()$
     GenSph();
     %%%;
Celles-ci sont qu'il existe une base de trois vecteurs orthonormés
'kr,kt,kp' qui sont tels que
     %%%; nettoie()$
     kr= sin(theta)*(cos(psi)*kx+sin(psi)*ky)+cos(theta)*kz@
     kt= cos(theta)*(cos(psi)*kx+sin(psi)*ky)-sin(theta)*kz@
     kp: -sin(psi)*ky+cos(psi)*ky@
     %%%;
où 'kx,ky,kz' sont les vecteurs de la base cartésienne sur laquelle
sont construite les coordonnées cylindiques ; et que la position
générale est décrite par les coordonnées sphériques 'r,theta,psi'
par
     %%%; r*kr@ %%%;

On peut utiliser ces coordonnées pour faire des calculs algébriques
comme différentiels.  Lorsque les calculs sont différentiels il est
nécessaire de déclarer la dépendance des variables utilisée par le
'depends' de Maxima Functions and Variables for Differentiation (Maxima).
Comme ceci
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,psi],v)$
     depends([phi,vr,vtheta,vpsi],[r,theta,z])@
     %%%;
On a alors accès aux expressions des opérateurs vectoriels en
coordonnées cylindriques
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,psi],v)$ GenSph()$
     depends([phi,vr,vtheta,vpsi],[r,theta,psi])$
     vv:vr*kr+vtheta*kt+vpsi*kp@
     grad(phi)=rtav(grad(phi))@
     div(vv)=tav(div(vv))@
     rot(vv)=rtav(rot(vv))@
     Grad(vv)=tav(Grad(vv))@
     %%%;
Si on préfère une notation matricielle pour 'Grad(vv)', on peut en
disposer avec
     %%%; nettoie()$ addvV(v)$ Indicie([r,theta,psi],v)$ GenSph()$
     depends([vr,vtheta,vpsi],[r,theta,z])$
     vv:vr*kr+vtheta*kt+vpsi*kp@
     Grad(vv)=tav(Grad(vv))@
     TensMat(tav(Grad(vv)),[kr,kt,kp])@
     %%%;

2.4.3.4 Coordonnées générales
................................

*GenCoor(Psi,uvw,kuvw)* génère les règles ad hoc pour des
coordonnés quelconques (orthogonales ou non), ses arguments sont
   * 'uvw' est de la forme '[u,v,w]', 'u,v,w' sont les coordonnées ;
   * 'kuvw' est de la forme '[ku,kv,kw]', 'ku,kv,kw' sont des vecteurs
     destinés à être les vecteurs de base ;
   * 'Psi' est une fonction linéaire de 'kx,ky,kz' (la base
     cartésienne sous-jacent) et pas nécessairement linéaire des
     '[u,v,w]' ; elle est a priori quelconque (sous réserve que 'Psi'
     soit un difféomorphisme) de 'u,v,w' dans E3.

Les relations de passage entre les 'ku,kv,kw' et 'kx,ky,kz' sont
rangées dans une variable globale : 'ChBs'

Toutes les régles de transformation sont calculées dans la base
'kx,ky,kz' puis mises dans la base 'ku,kv,kw' au moyen de ces
relations de passage.

Par exemple, les règles pour les coordonnées cylindriques sont
générées par
     %%%; nettoie()$
     assume(r>0)$
     Regles:GenCoor(r*(cos(theta)*kx+sin(theta)*ky)+z*kz,[r,theta,z],[kr,kt,kz])$
     PrintList(Regles,2)$
     %%%;

Les règles pour les coordonnées toriques (rayon 1), plus compliquée
même si elle restent orthogonales, sont générées par (attention
c'est un peu long)
     %%%; nettoie()$
     addv(K)$Indicie([X,Y,Z],K)$
     assume(r>0)$
     Regles:GenCoor((cos(theta)*kx+sin(theta)*ky)*(1+r*cos(phi))+r*sin(phi)*kz,[r,phi,theta],[kr,kp,kt])$
     PrintList(Regles,1)$
     %%%;

3 Fonctionnelles
****************

3.1 Qu'est ce qu'une fonctionnelle
==================================

Une fonctionnelle est une application d'un ensemble de champs de
scalaires, de vecteurs ou de dyades à valeur réelle.

Elle prend souvent la forme de l'intégrale sur un domaine du plan ou
de l'espace d'une expression qui comporte des éléments
différentiels du ou des champs dont elle est la fonctionnelle.  Par
exemple la mesure de la surface définie par
     %%%; nettoie()$ Indicie(base,S)$ addv(x)$
     DefEns(S,x+Z*kz,dans(s,sousensemble(Sbase,E2)))@
     %%%;
est la fonctionnelle de la hauteur 'Z' :
     %%%; nettoie()$ Indicie(base,S)$
     S=Int2D(sqrt(1+nrm(grad(Z))^2),Sbase)@
     %%%;

Pour considérer des expressions comportant des fonctionnelles, il
faut ajouter aux fonctions de manipulaition des expressions
algébrique d'autres fonctions qui opèrent sur les fonctionnelles.

Ce sont des telles fonctions qui sont décrites ici.  Elle laissent
tous les termes qui n'appartiennent pas à une fontionnelle
inchangés, se contentant de tranformer les fonctionnelles.

3.2 Manipulation de fonctionnelles
==================================

3.2.1 Variation de fonctionnelles
---------------------------------

Considérons cette fonctionnelle d'un champ scalaire sur un domaine D
     %%%; nettoie()$
     I:Int3D(nrm(grad(phi))^2/2,D)@
     %%%;

Si on suppose que 'phi' varie de 'dphi', pour la décoration (la
fonction 'deltise' est décrite Fonctions de décoration (NeoAppariteur))
     %%%; deltise(phi)$ dphi %%%;
alors la fonctionnelle 'I' varie de 'dI'.  Et on peut obtenir cette
variation à l'ordre 1 en 'dphi' en utilisant 'derive'
Dérivation
     %%%; nettoie()$
     I:Int3D(nrm(grad(phi))^2/2,D)@
     deltise(phi,'I)$
     dI:derive(I,[[phi,dphi]])@
     %%%;

C'est la même chose si le champ est vectoriel
     %%%; nettoie()$ addv(v)$
     I:Int3D(nrm(rot(v))^2/2,D)@
     deltise(v,'I)$
     dI:derive(I,[[v,dv]])@
     %%%;
ou même tensoriel (pour dyadique)
     %%%; nettoie()$ addy(gamma)$
     I:Int3D(gamma^2/2,D)@
     deltise(gamma,'I)$
     dI:derive(I,[[gamma,dgamma]])@
     %%%;

3.2.2 Intégration par parties
------------------------------

L'intégration par partie d'un fonctionnelle variée Variation de fonctionnelles
se fait avec la fontion *InPart3D(dI,dv)* où 'dI'
est la fonctionnelle et 'dv' est la variable qui doit être
intégrée, toutes les autres seront dérivées.  Par exemple
     %%%; nettoie()$
     I:Int3D(nrm(grad(phi))^2/2,D)@
     deltise(phi,'I)$
     dI:derive(I,[[phi,dphi]])@
     dI:IntPart3D(dI,dphi)@
     dI;
     %%%;
Cette intégration fait apparaître un terme de bord 'Bord(D)' ; comme
on suppose implicitement que le domaine 'D' est 3D, le 'Bord(D)' est
une surface qui est dotée d'une normale extérieure qu'on note
'Next'.

On a aussi
     %%%; nettoie()$ addv(v)$
     I:Int3D(nrm(rot(v))^2/2,D)@
     deltise(v,'I)$
     dI:derive(I,[[v,dv]])@
     dI:IntPart3D(dI,dv)@
     dI;
     %%%;

Mais là ça ne marche pas encore (je n'ai pas traité le cas)...
     %%%; nettoie()$ addv(v)$
     I:Int3D(Grad(v)^2/2,D)@
     deltise(v,'I)$
     dI:derive(I,[[v,dv]])@
     dI:IntPart3D(dI,dv)@
     dI;
     %%%;

3.2.3 Séparer et rassembler domaines et intégrants
----------------------------------------------------

*EclateIntegrant(expr)* sépare l'intégrant des fonctionnelles
contenues dans expr
     %%%; nettoie()$ addv(v)$
     I:Int3D(grad(phi)^2/2+q*phi,D)@
     EclateIntegrant(I)@
     %%%;

*EclateDomaine(expr)* sépare les domaines des fonctionnelles
contenues dans expr
     %%%; nettoie()$ addv(v)$ Indicie([1,2],D)$
     I:Int3D(grad(phi)^2/2+q*phi,D1+D2)@
     EclateDomaine(I)@
     %%%;

*CompacteInt(expr)* compacte les intégrales de même domaine
     %%%; nettoie()$ addv(v)$ Indicie([1,2],D)$
     I:Int3D(grad(phi)^2/2+q*phi,D1+D2)@
     I:EclateIntegrant(I)@
     CompacteInt(I)@
     %%%;

*DeleteIntDom(expr,D)* élimine les fonctionnelles de domaine 'D' dans
'expr'
     %%%; nettoie()$ Indicie([1,2],D)$
     I:Int3D(nrm(grad(phi))^2/2,D1)+Int3D(nrm(grad(phi))^2/2,D2)@
     DeleteIntDom(I,D2)@
     %%%;

*RentreCste(expr)* rentre les constantes multiplicatives des
fonctionnelle dans leur intégrant
     %%%; nettoie()$ Indicie([1,2],D)$
     lCste:[k]$
     I:3*Int3D(nrm(grad(phi))^2/2,D1)+k*Int3D(nrm(grad(phi))^2/2,D2)@
     RentreCste(I)@
     %%%;
Pour indiquer qu'un symbole est une constante, le ranger dans la liste
'lCste' (qui est vide au démarrage).

*SortCste(expr)* sort les constantes multiplicatives de l'intégrant
des fonctionnelles
     %%%; nettoie()$ Indicie([1,2],D)$
     lCste:[k]$
     I:Int3D(3*nrm(grad(phi))^2/2,D1)+Int3D(k*nrm(grad(phi))^2/2,D2)@
     SortCste(I)@
     %%%;

3.2.4 Évalue analytiquement la valeur d'une fonctionnelle
---------------------------------------------------------

*EvInt(expr,defdom1,deffom2,...)*  évalue analytiquement la valeur
des fonctionnelles dont les domaines sont spécifiées par les
'defdom's qui dont de la forme
'[D,[u,umin,umax],[v,vmin,vmax],[w,wmin,wmax],detj]' où 'D' est le
nom de domaine, '[u,v,w]' des coordonnées permettant le paramétrage
du domaine comme un pavé compris entre les bornes spécifiées et
finalement 'D' est le déteminant du jacobien pour ces coordonnées.
Par exemple
     %%%; nettoie()$
     IV:Int3D(x+y+z,D)@
     IV=EvInt(IV,[D,[[x,0,a],[y,0,b],[z,0,c],1]])@
     %%%;
ou
     %%%; nettoie()$
     IV:Int3D(1,D)@
     IV=EvInt(IV,[D,[[r,0,R],[theta,0,2*%pi],[z,0,Z],r]])@
     %%%;
ou
     %%%; nettoie()$ Indicie([i,x],R)$
     IV:Int3D(1/sqrt(r^2+z^2),D)@
     assume(Z>0,Rx>0,Ri>0,Rx-Ri>0)$
     IV=EvInt(IV,[D,[[r,Ri,Rx],[theta,0,2*%pi],[z,-Z/2,Z/2],r]])@
     %%%;

Il arrive (souvent) que l'expression analytique des fonctionnelles ne
puissent être déterminées.  Dans ce cas la sortie de 'EvInt' prend
la forme d'intégrales.  Par exemple
     %%%; nettoie()$
     IV:Int3D(1/sqrt(1+x^2+y^2),D)@
     assume(a>0,b>0)$
     IV=EvInt(IV,[D,[[x,0,a],[y,0,b],[z,0,c],1]])@
     %%%;
On peut alors obtenir les valeurs de ces intégrales (qui ne doivent
plus contenir de valeurs symbolique) avec 'Int2Quag' ou 'Int2Romberg'
Fonctions de service.  Par exemple
     %%%; nettoie()$
     IV:Int3D(1/sqrt(1+x^2+y^2),D)@
     assume(a>0,b>0)$
     res:IV=EvInt(IV,[D,[[x,0,a],[y,0,b],[z,0,c],1]])@
     res:ev(res,a=1,b=1)@
     res:Int2Quag(res)@
     res:ev(res)@
     %%%;

3.3 Domaines de l'espace et du plan
===================================

3.3.1 Quelques explications sur les domaines
--------------------------------------------

Symboliquement, un domaine est un nom comme 'D'.  On peut lui
conférer des propriétés comme la connexité et la simple
connexité.  Ou encore la convexité.  Et c'est très précieux pour
identifier les propriétés des champs de vecteurs dont un domaine
donné est le support.  Mais cela ne permet toutefois pas de
construire explicitement le domaine.

Cette construction explicite est faite à partir du mailleur Bamg de
FreeFem++ qui est un mailleur 2D, c'est à dire un mailleur de
domaines du plan.  Les maillages 2D se construisent par
'la donnée de lignes'
     FFBorderRec,FFBorderArc,FFBorderCou,FFBorderCir
'l'instruction de maillage'
     FFMesh
'la récupération des identificateurs de domaines'
     FFIds

Les lignes sont orientées, soit : pour les segments et arcs de la
1ière extrémité à la seconde ; pour les cercles, le sens
trigonométrique ; pour les courbes, le sens croissant du paramètre.

Cette orientation peut être changée en ajoutant 'sens=-1' dans la
définition des lignes.

Les lignes sont affectées d'un identificateur de ligne (un même
identificateur peut être attribué à plusieurs lignes).  Les
domaines sont définis comme des régions connexes fermées par les
lignes et chaque domaine a son identificateur.  Pour récupérer
ceux-ci on utilise la tactique de renvoyer l'identificateur des
éléments situés à gauche d'une des lignes qui le bordent, cette
ligne étant repérée par son identificateur (il faut donc faire
attention au choix des identificateurs de lignes).

Lors de la construction du maillage, on peut s'aider de tracés des
lignes (par FFPlot) et la colorisation des régions par
identificateurs (FFVisu).

Exemple simple
--------------

     %%%; nettoie()$
     FFInit()$        /* initialisation */
     FFDensite(10,5)$ /* définition de la densité de maillage et du nombre de décimales dans les sorties */
     FFDeclare('int,'isud=1,'iest=2,'inord=3,'iouest=4)$
     FFBorderRec('b1,[-1,-1],[1,-1],'isud)$
     FFBorderRec('b2,[1,-1],[1,1],'iest)$
     FFBorderRec('b3,[1,1],[-1,1],'inord)$
     FFBorderRec('b4,[-1,1],[-1,-1],'iouest)$
     FFPlot(apply("+",map(lambda([%u],apply(%u[1],[%u[2]])),lBorder)))$
     FFMesh(m)$
     FFIds(m,['idom,'isud])$
     FFDeclare('real,'nombre=1.0/3)$
     FFCout('idom,'nombre)$
     FFPlot(m)$
     FFgo()$
     %%%;

On conviendra que ce travail est pénible à faire, aussi CDE
possède-t'elle une fonction globale qui réalise toutes ces tâches.
C'est 'FaitMaillage'.

3.3.2 Explications pour FaitMaillage
------------------------------------

Voici la structure des arguments de *FaitMaillage*
     FaitMaillage('nom=nom_du_maillage,
                  'densite=densite_de_maillage,
                  'precision=precision_d_affichage,
                  'reels=[reel1=sa_valeur,reel2=...,...],
                  'entiers=[entier1=sa_valeur,entier2=...,...],
                  'points=[p1=[coord_x,coor_y],p2=...,...],
                  'bords=[[identificateur_du_bord,type_de_bord,le reste dépend du type],[...],...],
                  'epidemies=[[identificateur_de_bord,identificateur_de_domaine],[...],...],
                  'visus=[type_de_visualisation,...])
Ses arguments peuvent être placés dans n'importe quel ordre et
certains d'entre eux sont optionnels.  Voici leur description
   * 'nom_du_maillage' est le nom de la variable qui contiendra le
     maillage ;

   * 'densite_de_maillage' est un nombre tel que son inverse
     correspond à la dimension moyenne des triangles du maillage ;

   * 'precision_d_affichage' est le nombre de chiffres après la
     virgule qui seront affichés par FreeFem++ dans ses sorties
     numériques (optionnel) ;

   * pour 'reels : c'est un tableau d'affectation de variables
     réelles ; il est possible d'utiliser les réels précédemment
     définis dans les définitions, commme
          [r1=1,r2=2,r3=r1+r1]

   * Pour 'entiers : même commentaire que pour 'reels ;

   * Pour 'points : c'est un peu plus compliqué.  Les éléments de
     'points peuvent être de la forme
          pn=[a,b]
     où a et b sont soit des nombres, soit des réels définis dans
     'reels.

     Il peuvent aussi être de la forme
          pn=[a,b]+c*pq
     où pq est défini dans 'points (pas nécessairement avant
     d'ailleurs).  Et où l'expression qui définit pn peut être plus
     complexe que celle indiquée.

     Et finalement il peuvent être de la forme
          pn=Cou(identificateur,un_ou_deux)
     où identificateur est l'identificateur des bords définis
     ci-dessous et un_ou_deux prend la valeur 1 ou 2.  Dans ce cas le
     bord correspondant à identificateur doit être de type 'Cou pour
     « courbe paramétrées » et les valeurs renvoyées par
     Cou(identificateur,un_ou_deux) seront celle de la courbe
     paramétrée à son extrémité de départ ou d'arrivée suivant
     que un_ou_deux vaut 1 ou 2.

   * Pour 'bords : c'est encore un peu compliqué

     Les éléments de 'bords peuvent être de la forme
          [in,'Rec,pn,pm]
     qui correspond à un segment rectiligne compris entre les points
     pn et pm (des points de 'points ou la donnée de coordonnées
     sous la forme [a,b]).  in est l'identificateur de ce bord.

     Il peuvent encore être de la forme
          [in,'Arc,pn,R,theta2,theta2]
     qui correspond à un arc de cercle de centre pn, de rayon R et
     délimité par les angles theta1 et theta2.  in est
     l'identificateur de ce bord.  in est l'identificateur de ce bord.

     Ou de la forme
          [in,'Cer,pn,R]
     c'est comme pour 'Arc mais avec le cercle entier.

     Et finalement de la forme
          [in,'Cou,[x(s),y(s)],[s,smin,smax]]
     Dans ce cas le bord sera la partie de la courbe paramétrée
     [x(s),y(s)] telle que smin<s<smax.  La distribution des points de
     cette partie de courbe est faite de manière qu'ils soient à peu
     près équidistant de distance 1/densite_de_maillage.

     Il s'ajoute encore une petite complication.  Les listes
     définissant les bords de n'importe quel type peuvent être
     augmenté de -1 comme
          [in,'Rec,pn,pm,-1]
     et alors ce bord sera exactement comme s'il avait été écrit
          [in,'Rec,pm,pn]
     Ça sert à changer la direction dans la définition du bord.

   * Pour 'epidémie : ça reste compliqué

     Un élément de 'epidemie de la forme
          [identificateur_de_bord,identificateur_de_domaine]
     a pour action que la variable identificateur_de_domaine sera
     affectée de l'identificateur de région que FreeFem++ aura
     donné aux triangles situés sur la gauche de n'importe lequel
     des bords qui auront été affectés de l'identificateur «
     identificateur_de_bord ».  Ça sert à repérer les différentes
     régions d'un maillage.

   * Pour 'visus : c'est plus simple.

     Si on met 'bords dans la liste, FreeFem++ affichera les bords du
     maillage (avant de tenter de créer le maillage) ; si on met
     'maillage il affichera le maillage ; si on met 'domaines il
     affichera une visualisation des régions du maillage qui auront
     été repérées par 'epidemie ; et si on met 'info il affichera
     une série d'informations sur le maillage produit.

3.3.3 Exemples simples pour FaitMaillage
----------------------------------------

Pour chacun de ces maillages, on peut voir le fichier FreeFem++ par
'<CTRL>-w' : il apparait hélas dans whg, le mieux est donc de mettre
*info* dans wbg.

** Un carré
     %%%;
     FaitMaillage('nom=m,'densite=10,
     'points=[p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1]],
     'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1]],
     'visus=['info,'bords,'maillage])$
     FFgo()$
     %%%;

** Un carré avec une diagonale
     %%%;
     FaitMaillage('nom=m,'densite=10,
     'points=[p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1]],
     'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
             [diag,'Rec,p1,p3]],
     'epidemie=[[gauchediag,diag],[gauchesud,sud]],
     'visus=['info,'bords,'maillage,'domaines])$
     FFgo()$
     %%%;

** Un carré avec deux diagonales
     %%%;
     FaitMaillage('nom=m,'densite=10,
     'points=[p0=(p1+p2+p3+p4)/4,p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1]],
     'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
             [som,'Rec,p1,p0],[mne,'Rec,p0,p3],[sem,'Rec,p2,p0],[mno,'Rec,p0,p4]],
     'epidemie=[[dsud,sud],[dnord,nord],[dest,est],[douest,ouest]],
     'visus=['info,'bords,'maillage,'domaines])$
     FFgo()$
     %%%;

** Un carré avec un cercle au centre
     %%%;
     FaitMaillage('nom=m,'densite=20,
     'points=[p0=(p1+p2+p3+p4)/4,p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1]],
     'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
             [cercle,'Cer,p0,0.2]],
     'epidemie=[[dext,sud],[dint,cercle]],
     'visus=['info,'bords,'maillage,'domaines])$
     FFgo()$
     %%%;

** Un carré avec un trou au centre (le sens de parcours du cercle est
inversé)
     %%%;
     FaitMaillage('nom=m,'densite=20,
     'points=[p0=(p1+p2+p3+p4)/4,p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1]],
     'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
             [cercle,'Cer,p0,0.2,-1]],
     'epidemie=[[dext,sud]],
     'visus=['info,'bords,'maillage])$
     FFgo()$
     %%%;

** Un haricot

     %%%;
     FaitMaillage('nom=m,'densite=10,
     'reels=[R=1],
     'points=[p1=[-2*R,0],p2=[0,0],p3=[2*R,0]],
     'bords=[[ib,'Arc,p2,3*R,pi,2*pi],[ib,'Arc,p3,R,0,pi],
             [ib,'Arc,p2,R,pi,2*pi,-1],[ib,'Arc,p1,R,0,pi]],
     'visus=['info,'bords,'maillage])$
     FFgo()$
     %%%;

** Un truc (on remarque que la densité de points sur le bord reste
constante)
     %%%;
     FaitMaillage('nom=m,'densite=20,
     'reels=[R=1],
     'points=[p1=[-2*R,0],p2=[0,0],p3=[2*R,0]],
     'bords=[[ib,'XY,[cos(s)+0.4*cos(3*s),sin(s)+0.4*cos(3*s)],[s,0,2*%pi]]],
     'visus=['info,'bords,'maillage])$
     FFgo()$
     %%%;

3.4 Stationnarisation de fonctionnelles
=======================================

3.4.1 À un champ scalaire
-------------------------

On considère la fonctionnelle
     %%%; nettoie()$
     W:Int3D(nrm(grad(phi))^2/2+k*phi^2/2+q*phi,D)@
     %%%;
dont on cherche la valeur minimale (dans le cas 'k>0') sur un domaine
qui est par exemple celui du haricot Exemples simples pour FaitMaillage
et lorsque la valeur de 'phi' est nulle sur le bord de
cedomaine.

La condition nécessaire d'extréma est obtenue par variation de la
fonctionnelle Variation de fonctionnelles comme
     %%%; nettoie()$
     W:Int3D(nrm(grad(phi))^2/2+k*phi^2/2+q*phi,D)@
     deltise(phi,'W)$
     dW:derive(W,[[phi,dphi]])@
     %%%;
Il faudra que pour tout 'dphi', 'W' soit nul.

Le passage en coordonnées (cartésiennes) se fait par
Coordonnées cartésiennes
     %%%; nettoie()$
     W:Int3D(nrm(grad(phi))^2/2+k*phi^2/2+q*phi,D)@
     deltise(phi,'W)$
     dW:derive(W,[[phi,dphi]])@
     depends([phi,dphi],[x,y])@
     GenCart()$
     WFF:xtav(W)@
     dWFF:xtav(dW)@
     %%%;

Il reste à demander à FreeFem++ de réaliser la condition
nécessaire (le domaine 'D' étant celui du haricot).  Pour cela
     %%%;  nettoie()$
     W:Int3D(nrm(grad(phi))^2/2+k*phi^2/2+q*phi,D)$ deltise(phi,'W)$
     dW:derive(W,[[phi,dphi]])$ depends([phi,dphi],[x,y])$
     GenCart()$
     WFF:xtav(W)$ dWFF:xtav(dW)$

     FaitMaillage('nom=m,'densite=10,
     'reels=[R=1],
     'points=[p1=[-2*R,0],p2=[0,0],p3=[2*R,0]],
     'bords=[[ib,'Arc,p2,3*R,pi,2*pi],[ib,'Arc,p3,R,0,pi],
             [ib,'Arc,p2,R,pi,2*pi,-1],[ib,'Arc,p1,R,0,pi]],
     'visus=['info,'bords,'maillage])$

     /* à partir de là */

     cvs:[[D,[]],[detJ,1]]$
     FFDeclare(mP1,phi,dphi)$
     FFDeclare(real,k=1,q=1)$
     FFSolve(pb,[phi,dphi],FtoFF(m,dWFF,cvs),[On(ib,phi=0)])$
     FFPlot(phi)$
     FFDeclare('real,'Wmin=FtoFF(m,WFF,cvs))$
     FFCout('Wmin)$
     FFgo()$
     %%%;
Ce script remplit un buffer 'freefem' qui apparaît (ou sinon on peut
le faire apparaître) dans 'whd' avec les sorties de FreeFem++

Il fait également apparaitre la fenêtre graphique de 'ffglut'.

Les lignes situées au dessus de 'à partir de là' sont commentées
plus haut.  Celles qui viennent ansuite sont
FtoFF (NeoAppariteur) et FFSolve (NeoAppariteur)
   * 'cvs:[[D,[]],[detJ,1]]$' C'est une variable qui contient les
     éléments de conversion entre les symboles de la fonctionnelle
     et ceux du maillage.

     Ici, c'est assez simple : le domaine 'D' est associé à '[]', ce
     qui signifie que ce domaine est celui du maillage entier.  Et le
     seul élément ajouté est le déterminant du jacobien qui vaut 1
     (en coordonnées cartésiennes).

   * 'FFDeclare(mP1,phi,dphi)$' les variables 'phi' et 'dphi' de la
     fonctionnelle seront approximées par des éléments finis de
     Lagrange P1 sur le maillage 'm'.  See
     FFDeclare (NeoAppariteur)

   * 'FFDeclare(real,k=1,q=1)$' les variables 'k' et 'q' seront des
     réels de valeur 1.

   *
     'FFSolve(pb,[phi,dphi],append([FtoFF(m,dWFF,cvs)],[On(ib,phi=0)]))$'

     La fonction *FFSolve* est décrite FFSolve (NeoAppariteur)
        * 'pb' : c'est le nom du problème

        * '[phi,dphi]' le premier argument de la liste est champ à
          calculer ; le second sa variation.

        * 'FtoFF(m,dWFF,cvs)', la fonction 'FtoFF' convertit la forme
          bilinéaire (en phi et dphi) et linéaire (en phi) dWFF en
          une expression presque conforme à la syntaxe de FreeFem++
               %%%; FtoFF(m,dWFF,cvs); %%%;

        * 'On(ib,phi=0)', c'est la syntaxe de FreeFem++ pour indiquer
          une condition de Dirichlet sur 'ib' qui est l'indicateur de
          bord du domaine de calcul.

   * 'FFPlot(phi)$' on demande à FreeFem++ (plutôt 'ffglut'
     d'ailleurs) de tracer les lignes de niveau de 'phi'
     FFPlot (NeoAppariteur)

   * 'FFDeclare('real,'Wmin=FtoFF(m,WFF,cvs))$' on range dans la
     variable 'Wmin' la valeur de la fonctionnelle 'WFF'

   * 'FFCout(Wmin)' on demande à FreeFem++ d'afficher la valeur de
     cette variable.  See FFCout (NeoAppariteur)

   * 'FFgo()' c'est l'ordre de lancement de FreeFem++ (les fonctions
     ci-dessus n'ont fait que remplir 'lScript'
     Interface de FreeFem++ (NeoAppariteur))

On notera que les fonctions dont le nom commencent par 'FF' sont très
proches des fonctions du langage de FreeFem++ (de nom, celui qui suit
FF). Ces fonction ne font guère que de la convesion syntaxique.

Il est donc utile, si on ne connaît pas bien la syntaxe de FreeFem++,
de se référer à son manuel
<https://doc.freefem.org/documentation/index.html>

3.4.2 À deux champ scalaires
----------------------------

On considère la fonctionnelle
        %%%; nettoie()$
          W:Int3D(ps(grad(Omega),grad(psi))-Omega^2/2,D)@
         %%%;
où 'D' est le domaine situé entre deux cylindres emboités l'un dans
l'autre soit
     %%%; nettoie()$
     FaitMaillage('nom=m,'densite=10,
          'reels=[Ri=1,Rx=2,Dec=0.3],
          'points=[ci=[Dec,0],cx=[0,0]],
          'bords=[[ix,'Cer,cx,Rx],[ii,'Cer,ci,Ri,-1]],
          'epidemie=[[iD,ix]],
          'visus=['info,'bords,'maillage,'domaines])$
     FFgo()$
     %%%;
On appelle 'Gammai' et 'Gammax' les bords intérieur et extérieur de
'D'.

On donne des conditions sur 'psi' qui sont

     %%%; nettoie()$ Indicie([i,x],Gamma)$
     ListeTex(DANS(psi=1,Gammai), DANS(psi=0,Gammax), SUR(ps(grad(psi),Next),Bord(D)=Gammai+Gammax))@
     %%%;

Et on cherche les valeurs de 'Omega' qui stationnarisent la
fonctionnelle, soit

     %%%; nettoie()$ primise(Omega,psi)$
     W:Int3D(ps(grad(Omega),grad(psi))-Omega^2/2,D)$
     Wp:derive(W,[[psi,psip],[Omega,Omegap]])$
     DefEns(V,psi,dans(psi=0,Bord(D)))@
     PourTout(ListeTex(Omegap,dans(psip,V)),Wp=0)@
     %%%;

Si les deux cylindres étaient centrés, on trouverait une solution
analytique celle de l'écoulement de Couette en 2D. Mais ils ne le
sont pas nécessairement, ce qui fait qu'une solution analytique reste
probablement possible (les transformations conformes, par exemple font
des miracles ou encore les coordonnées bipolaires) mais est plus
difficile à obtenir.

On décide donc d'utiliser des moyens numériques.

     %%%; nettoie()$ primise(psi,Omega)$
     W:Int3D(ps(grad(Omega),grad(psi))-Omega^2/2,D)$
     Wp:derive(W,[[psi,psip],[Omega,Omegap]])$
     FaitMaillage('nom=m,'densite=10,
          'reels=[Ri=1.0,Rx=2,Dec=0.5],
          'points=[ci=[Dec,0],cx=[0,0]],
          'bords=[[ix,'Cer,cx,Rx],[ii,'Cer,ci,Ri,-1]],
          'epidemie=[[iD,ix]],
          'visus=['info,'bords,'maillage,'domaines])$
     cvs:[[D,iD],[detJ,1]]$
     FFDeclare(mP1,Omega,psi,Omegap,psip)$
     GenCart()$ depends([Omega,psi,Omegap,psip],[x,y])$
     FFWp:xtav(Wp)@
     FFSolve(prob,[[psi,Omega],[psip,Omegap]],FtoFF(m,FFWp,cvs),On(ix,psi=0),On(ii,psi=1))$
     FFPlot(Omega)$
     FFPlot(psi)$
     FFgo()$
     %%%;
...

3.4.3 À un champ de vecteurs
----------------------------

...

4 Équations aux Dérivées Partielles
*************************************

4.1 Problèmes elliptiques
==========================

4.1.1 EDP elliptiques et Fonctionnelles
---------------------------------------

Les EDP elliptiques sont emblématiques des EDP qui peuvent se
reformuler comme la recherche d'un extréma de fonctionnelles.

Si on pose par exemple le problème
     %%%; nettoie()$ depends(phi,[n,tau])$
     Etages(DANS(div(eta*grad(phi)) - k*phi + q = 0,D),SUR(phi=phi0,Gamma0),SUR(-eta*diff(phi,n)=p0+h0*phi,Gamma1),POD("\\mbox{avec }",Bord(D)=Gamma0+Gamma1))@
     %%%;
(ce n'est pas tout à fait une EDP elliptique du fait de la présence
de 'k*u' mais on la considère comme telle), la solution est
nécessairement solution de

     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$
     PourTout(phip,Int3D(phip*(div(eta*grad(phi)) - k*phi + q),D)=0)@
     %%%;
Si on pose
     %%%; nettoie()$ addv(Psi)$
     Psi=eta*grad(phi)@
     %%%;
Cette relation devient
     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$ addv(Psi)$
     PourTout(phip,Int3D(phip*(div(Psi) - k*phi + q),D)=0)@
     %%%;
soit, après intégration par parties ('Psi' intégrée),
     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$ addv(Psi)$
     IntPart3D(Int3D(phip*(div(Psi) + k*phi + q),D),Psi)@
     %%%;
et donc
     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$ addv(Psi)$
     res:IntPart3D(Int3D(phip*(div(Psi) + k*phi + q),D),Psi)$
     PourTout(phip,ev(res,Psi=eta*grad(phi))=0)@
     %%%;
Pour prendre en compte les conditions aux limites, il faut séparer le
bord de D en ses deux parties
     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$ addv(Psi)$ Indicie([0,1],Gamma)$
     res:IntPart3D(Int3D(phip*(div(Psi) + k*phi + q),D),Psi)$
     res:ev(res,Psi=eta*grad(phi))$
     EclateDomaine(ev(res,Bord(D)=Gamma0+Gamma1))@
     %%%;
On choisit de ne garder que les 'phi'' qui sont nuls sur 'Gamma0' et,
on a
     %%%; nettoie()$ depends(phi,[n,tau])$
     SUR(MultEgal(eta*ps(Next,grad(phi)),eta*diff(phi,n),-p0-h0*phi),Gamma1)@
     %%%;
d'où
     %%%; nettoie()$ depends(phi,[n,tau])$ primise(phi)$ addv(Psi)$ Indicie([0,1],Gamma)$
     res:IntPart3D(Int3D(phip*(div(Psi) + k*phi + q),D),Psi)$
     res:ev(res,Psi=eta*grad(phi))$
     res:EclateDomaine(ev(res,Bord(D)=Gamma0+Gamma1))$
     res:DeleteIntDom(res,Gamma0)$
     PourTout(dans(phip,DefEns("",psi,SUR(psi=0,Gamma0))),ev(res,ps(Next,grad(phi))=(-p0-h0*phi)/eta)=0)@
     %%%;
C'est exactement la variation à l'ordre 1 de la fonctionnelle (pour
une condition de Dirichlet sur 'Gamma0')
     %%%; nettoie()$ Indicie([0,1],Gamma)$
     Int3D(phi*q+k*phi^2/2-eta*nrm(grad(phi))^2/2,D)-Int2D(phi*p0+h0*phi^2/2,Gamma1)@
     %%%;
En effet
     %%%; nettoie()$ primise(phi)$ Indicie([0,1],Gamma)$
     W:Int3D(phi*q+k*phi^2/2-eta*nrm(grad(phi))^2/2,D)-Int2D(phi*p0+h0*phi^2/2,Gamma1)$
     derive(W,[[phi,phip]])@
     %%%;

Le calcul de la solution de l'EDP se fait donc par extrémisation de
la fontionnelle correspondante.  Et cela de la manière indiquée
Stationnarisation de fonctionnelles.

Pour cet exemple, en reprenant un haricot un peu modifié pour qu'il y
ait deux identificateurs de bord différents (correspondant aux deux
parties du bord de D)

     %%%;  nettoie()$ Indicie([0,1],Gamma)$
     W:Int3D(phi*q+k*phi^2/2-eta*nrm(grad(phi))^2/2,D)-Int2D(phi*p0+h0*phi^2/2,Gamma1)@
     dW:derive(W,[[phi,dphi]])$ depends([phi,dphi],[x,y])$
     GenCart()$
     WFF:xtav(W)$ dWFF:xtav(dW)$

     FaitMaillage('nom=m,'densite=10,
     'reels=[R=1],
     'points=[p1=[-2*R,0],p2=[0,0],p3=[2*R,0]],
     'bords=[[ib0,'Arc,p2,3*R,pi,2*pi],[ib1,'Arc,p3,R,0,pi],
             [ib1,'Arc,p2,R,pi,2*pi,-1],[ib1,'Arc,p1,R,0,pi]],
     'visus=['info,'bords,'maillage])$

     cvs:[[D,[]],[Gamma1,ib1],[Gamma0,ib0],[detJ,1]]$
     FFDeclare(mP1,phi,dphi)$
     FFDeclare(mP0,eta=4-x)$
     FFDeclare(real,k=0,q=0,h0=0,p0=1)$
     FFSolve(pb,[phi,dphi],FtoFF(m,dWFF,cvs),[On(ib0,phi=0)])$
     FFPlot(phi)$
     FFDeclare('real,'Wmin=FtoFF(m,WFF,cvs))$
     FFCout('Wmin)$
     FFgo()$
     %%%;
On a déclaré 'eta' comme un champ (approximé en Lagrange P0) pour
le faire dépendre de la position.  On aurait pu faire cela également
pour les 'k', 'q'.

4.1.2 Fonctionnelles usuelles
-----------------------------

4.1.2.1 Elasticité linéaire
.............................

Pour l'élasticité linéaire, c'est
     %%%; nettoie()$ addvV(u)$
     gu:Grad(vu)$
     gamma:(gu+Trans(gu))/2$
     sigma:lambda*Trace(gamma)*Id+2*mu*gamma$
     W:xtav(Int3D(sigma*gamma/2,D))@
     %%%;
on obtient le passage en composantes par
     %%%; nettoie()$ addvV(u)$
     gu:Grad(vu)$
     gamma:(gu+Trans(gu))/2$
     sigma:lambda*Trace(gamma)*Id+2*mu*gamma$
     W:xtav(Int3D(sigma*gamma/2,D))$
     GenCart()$ Indicie([x,y],u)$ depends([ux,uy],[x,y])$
     WFF:xtav(ev(W,vu=ux*kx+uy*ky))@
     primise(ux,uy)$ depends([uxp,uyp],[x,y])$
     dWFF:derive(WFF,[[ux,uyp],[uy,uyp]])@
     %%%;
On notera que la fonction utilisée pour déclarer 'u' en vecteur est
*addvV* elle a pour effet que c'est 'vu' qui sera le vecteur 'u'.
L'intérêt est que pour le passage en composantes avec *Indicie* les
composantes 'ux' et 'uy' n'héritent pas du caractère de vecteur de
'u' puisqu'il n'en est pas un.

On notera aussi que, contrairement aux exemples précédents, la
fonctionnelle est variée après le passage en composantes.  On peut
faire l'un ou l'autre.

Pour le passage au calcul numérique d'un exemple particulier
Déformation d'un cylindre sous son poids

4.1.2.2 Fluide Newtonien incompressible
.......................................

Un fluide Newtonien est un milieu fluide dans lequel la relation entre
les tenseurs de vitesse de déformation et de contrainte est de la
forme
     %%%; nettoie()$ addy(sigma,gamma)$ dotise(gamma)$
     sigma=-P*Id+2*eta*gamma@
     %%%;
avec
     %%%; nettoie()$ addy(gamma)$ dotise(gamma)$ addvV(v)$
     vv@
     %%%;
L'inconpressibilité correspond à la condition
     %%%; nettoie()$ addvV(v)$
     div(vv)=0@
     %%%;

Si le fluide est placé dans un domaine 'D' bordé par une paroi
'Bord(D)' solide, alors la condition aux limites sur cette paroi est
(imperméabilité et non-glissement)
     %%%; nettoie()$ addvV(v)$
     SUR(vv=vzero,Bord(D))@
     %%%;

On suppose que le fluide est sollicité par une densité de force
extérieure 'f' dans le domaine 'D' et on cherche à caractériser la
fonction (peu importe ce que sont les espaces fonctionnels)
     %%%; nettoie()$ addvV(v)$ addvV(f)$
     DefFonc("",vL2,ProduitCart(H1,vH1),vf,""(P,vv))@
     %%%;
où 'P' et 'V' sont tels que
     %%%; nettoie()$ addvV(v)$ addvV(f)$ primiseV(vv)$ addy(sigma,gamma)$ primise(gammap)$
     PourTout(dans(vvp,H1),Int3D(sigma*gammapp+ps(vf,vv),D)=0)@
     %%%;
et où 'v' satisfait la condition aux limites d'imperméabilité et de
non glissement sur 'Bord(D)'

Voilà donc bien du verbillage !  Décodons cela par le calcul.

4.2 Problèmes paraboliques
===========================

4.2.1 Le problème emblématique
--------------------------------

Prenons un problème de diffusion.  Le champ des concentrations d'un
quantité est 'c' ; son support est un domaine 'D' dans lequel la
diffusivité est 'alpha' ; son bord comporte :
   * une partie 'Gamma1' sur laquelle on impose que la concentration
     est nulle ;
   * un partie 'Gamma0' sur laquelle on impose une condition de Robin
          %%%; -alpha*ps(Next,grad(c)) = h*(c-C)@ %%%;
     où 'C' est une constante donnée ;
   * et finalement sur le reste du bord on a une condition
     d'imperméabilité
          %%%; -alpha*ps(Next,grad(c)) = 0@ %%%;
et initialement la concentration dans 'D' est 'c0'.

Voici déjà la forme faible discrétisée en temps d'une façon
implicite
     %%%; nettoie()$
     Regles:[]$
     depends(c,[t])$primise(c)$declare(alpha,constant)$
     ff:Int3D(cp*(-alpha*div(grad(c))+diff(c,t)),D)@
     ff:IntPart3D(ff,c)@
     Indicie(0,c)$deltise(t)$
     ffd:eval(ff,diff(c,t)=(c-c0)/dt)@
     Indicie([0,1,L],Gamma)$
     ffd:EclateDomaine(eval(ffd,Bord(D)=Gamma0+Gamma1+GammaL))@
     ffd:eval(DeleteIntDom(ffd,Gamma1,GammaL),eval,ps(Next,grad(c))=-h*(c-C)/alpha)@
     %%%;
le pas de temps est 'dt'.

Si le domaine de calcul a une symétrie cylindrique (avec invariance
de toutes les quantités le long de la direction d'invariance
géométrique 'kx'), le problème devient 2D et le passage en
coordonnées donne
     %%%;
     GenCart()$ depends([c,cp],[x,y])$
     FFffd:expand(xtav(ffd))@
     %%%;
Et il ne reste plus qu'à déléguer le calcul de cette forme faible
à FreeFem++ (en incorporant la condition aux limites de Dirichlet sur
'Gamma1'.

On choisit un domaine de base carrée avec un trou en son centre, les
bords 'Gamma0' et 'Gamma1' étant deux côtés opposés du carré
     %%%;
     FaitMaillage('nom=m,'densite=20,
      'reels=[R=0.3],
      'points=[p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1],p5=(p1+p3)/2],
       'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
               [trou,'Cer,p5,R,-1]],
       'visus=['info,'bords,'maillage])$
     cvs:[[D,[]],[Gamma0,nord],[detJ,1]]$
     lTableau:[]$
     FFDeclare('real,alpha=1,h=1.0e+3,C=1,dt=0.01)$
     FFDeclare('mP1,c,cp,c0=0)$
     jmx:20$
     FFBoucle(j,j>jmx,jmx,'[
      FFSolve(pb,[c,cp],FtoFF(m,FFffd,cvs),[On(sud,c=0)]),
      FFAffecte(c0=c),
      FFPlot(c,cmm="c"+j*dt,fill=true,dim=2),
      FFDeclare('real,Q=FtoFF(m,Int3D(c,D),cvs)),
      FFCout(Q)
      ])$
     FFgo()$
     %%%;
Ce script illustre l'utilisation de 'FFBoucle'
FFBoucle (NeoAppariteur) pour effectuer les pas de temps.

4.2.2 La récupération de tableaux de valeurs
----------------------------------------------

On a ajouté dans la boucle précédente Le problème emblématique
le calcul de la quantité (dont 'c' est la
concentration) présente au cours du temps dans le domaine 'D'.

Et on souhaite récupérer ces valeurs de manière qu'on puisse la
tracer en fonction du temps.

Pour cela, on peut procéder de deux façon pour obtenir le tracé
désiré

4.2.2.1 Méthode avec récupération par Emacs
..............................................

Étant entendu que 'FFffd' a été calculé précédemment Le problème emblématique
et a la valeur
     %%%;
     FFffd:Int2D(c*cp*h-C*cp*h,Gamma0)+Int3D((-(c0*cp)/dt)+(c*cp)/dt +alpha*'diff(c,y,1) *'diff(cp,y,1)+alpha*'diff(c,x,1)*'diff(cp,x,1),D)$
     %%%;

     %%%;
     ok;
     FaitMaillage('nom=m,'densite=20,
      'reels=[R=0.3],
      'points=[p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1],p5=(p1+p3)/2],
       'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
               [trou,'Cer,p5,R,-1]],
       'visus=['info,'bords,'maillage])$
     cvs:[[D,[]],[Gamma0,nord],[detJ,1]]$
     lTableau:[]$
     FFDeclare('real,alpha=1,h=1.0e+3,C=1,dt=0.01)$
     FFDeclare('mP1,c,cp,c0=0)$
     jmx:20$
     FFDeclare(real[int],Temps(jmx),Concs(jmx))$
     FFVerbatim(printf(false,"real[int] ~a,~a;",Ttemps(jmx),Tconc(jmx)))$
     FFBoucle(j,j>jmx,jmx,'[
      FFSolve(pb,[c,cp],FtoFF(m,FFffd,cvs),[On(sud,c=0)]),
      FFAffecte(c0=c),
      FFPlot(c,cmm="c"+j*dt,fill=true,dim=2),
      FFDeclare('real,Q=FtoFF(m,Int3D(c,D),cvs)),
     /* FFCout(Ct),*/
      FFAffecte(Temps(j)=j*dt,Concs(j)=Q)
      ])$
     NFFCout(Temps)$
     NFFCout(Concs)$
     FFgo()$
     %%%;
Une fois que c'est fini, on peut alors faire
     %%%;
     plot2d(append([discrete],map(second,data_tst)),[xlabel,"temps"],[ylabel,"quantite"])$
     %%%;

4.2.2.2 Méthode avec récupération par Maxima
...............................................

Étant entendu que 'FFffd' a été calculé précédemment Le problème emblématique
et a la valeur
     %%%;
     FFffd:Int2D(c*cp*h-C*cp*h,Gamma0)+Int3D((-(c0*cp)/dt)+(c*cp)/dt +alpha*'diff(c,y,1) *'diff(cp,y,1)+alpha*'diff(c,x,1)*'diff(cp,x,1),D)$
     %%%;

     %%%;
     ok;
     FaitMaillage('nom=m,'densite=20,
      'reels=[R=0.3],
      'points=[p1=[0,0],p2=[1,0],p3=p2+p4,p4=[0,1],p5=(p1+p3)/2],
       'bords=[[sud,'Rec,p1,p2],[est,'Rec,p2,p3],[nord,'Rec,p3,p4],[ouest,'Rec,p4,p1],
               [trou,'Cer,p5,R,-1]],
       'visus=['info,'bords,'maillage])$
     cvs:[[D,[]],[Gamma0,nord],[detJ,1]]$
     lTableau:[]$
     FFDeclare('real,alpha=1,h=1.0e+3,C=1,dt=0.01)$
     FFDeclare('mP1,c,cp,c0=0)$
     jmx:20$
     FFDeclare(real[int],Temps(jmx),Concs(jmx))$
     FFVerbatim(printf(false,"real[int] ~a,~a;",Ttemps(jmx),Tconc(jmx)))$
     FFBoucle(j,j>jmx,jmx,'[
      FFSolve(pb,[c,cp],FtoFF(m,FFffd,cvs),[On(sud,c=0)]),
      FFAffecte(c0=c),
      FFPlot(c,cmm="c"+j*dt,fill=true,dim=2),
      FFDeclare('real,Q=FtoFF(m,Int3D(c,D),cvs)),
     /* FFCout(Ct),*/
      FFAffecte(Temps(j)=j*dt,Concs(j)=Q)
      ])$
     NFFCout(Temps)$
     NFFCout(Concs)$
     FFScript(tst)$
     res:FreeFem(tst)$
     plot2d(append([discrete],map(second,res)),[xlabel,"temps"],[ylabel,"quantite"])$
     ffglut(tst)$
     %%%;
Noter que l'appel à 'ffglut' doit être fait à la fin.  Tant que la
fenêtre graphique est active, 'Maxima' atteint qu'elle soit
arrêtée.

5 Équations Différentielles Ordinaires
***************************************

D'une façon générale, les Équation Différentielles Ordinaires
(EDO)s sont représentées par une structure '[VD,SM]' où 'SM' est
une liste contenant les seconds membres de l'EDO et VD la liste des
variables dépendantes.  La variable indépendante étant en
générale notée 't'.

Par exemple
     %%%; nettoie()$
     EDO:[[x,y],[-y,x+e(t)]];
     %%%;
est la représentation de
     %%%;
     Systeme(EDO)@
     %%%;
La fonction de restitution visuelle 'Systeme' ne fonctionnant
correctement que pour les EDOs linéaires.

5.1 Résolution analytique
==========================

En dehors des cas qu'on connait bien, les solveurs analytiques du
calcul symbolique sont plutôt décevants.  En tout cas pour Maxima.
Par contre le calcul symbolique est intéressant pour faire des
changements de variables.

5.1.1 Packaging de desolve et ode2
----------------------------------

Maxima possède deux solveurs d'EDO : 'desolve' et 'ode2'
Functions and Variables for Differential Equations (Maxima).  Ils ont
été encapsulés.

Pour 'desolve' par *EDOLin(SM,VD,VD0,t,t0)* où
     %%%; nettoie()$ ListeTex('diff(VD,t)=SM,VD(t=t0)=VD0)@ %%%;

Par exemple la solution de l'oscillateur harmonique est
     %%%;
     EDOLin([-y,x],[x,y],[x0,y0],t,0)@
     %%%;

Pour 'ode2' par *edo2(EQ,X,t,CICL)*

Par exemple la solution de l'oscillateur harmonique est (on ne met ni
de conditions initiales CI ni de conditions aux limites CL)
     %%%;
     edo2('diff(x,t,2)+x,x,t)@
     %%%; nettoie()$
On peut mettre des conditions initiales par
     %%%; nettoie()$ Indice(0,x)$ primise(x0)$
     res:edo2('diff(x,t,2)+x,x,t,[t0,x0,x0p])@
     trigrat(res)@
     %%%;
et des conditions aux limites par
     %%%; nettoie()$
     res:edo2('diff(x,t,2)+x,x,t,[t0,x0],[t1,x1])@
     trigrat(res)@
     %%%;

5.1.2 Changement de variables
-----------------------------

*ChangeTetV(eqs,chT,chV)* transforme 'eqs' en changeant la variable
indépendante VI suivant les indication de 'chT' qui a le format
'[ancienne VI,nouvelle VI,expression de l'ancienne VI en fonction de
la nouvelle VI]' et de 'chV' qui a le format, pour chaque variable
dépendante '[...[ancienne variable dépendante, expression de
l'ancienne variable dépendante en fonction de la nouvelle]...]'

Par exemple
     %%%; nettoie()$
     depends([x,v],t)$ assume(m>0,k>0)$
     eqs:[m*diff(v,t)=-k*x,diff(x,t)=v]$
     depends([xi,nu],tau)$
     res:ChangeTetV(eqs,[t,tau,tau*sqrt(m/k)],[[x,xi/L],[v,nu*alpha]])@
     res:solve(res,[diff(nu,tau),diff(xi,tau)])@
     %%%;
on a laissé un paramètre libre alpha et on voit bien que la valeur
'alpha=1/(L/sqrt(k/m))' rend le système parfaitement adimensionnel
     %%%; ev(res,alpha=1/(L/sqrt(k/m)))@ %%%;

Ça reste très peu automatisé, mais le changement de variable est a
priori plus un art qu'une science.  D'autant que son choix dépend des
objectifs qu'on vise.  Ceux-ci pouvant être un peu plus subtils que
l'adimensionnement lorsque le système étudié est plus complexe
qu'un oscillateur harmonique.

5.2 Méthodes de Rung-Kutta emboîtées
=======================================

5.2.1 Utilisation simpliste
---------------------------

L'utilisation la plus simple des fonctions d'intégration d'EDO est
illustrée par l'exemple de l'oscillateur harmonique

     %%%; nettoie()$
     EDO:[[x,y],[-y,x],[w=x^2+y^2]];
     Prepare(EDO);
     X0:[1,0];
     tol:1.0e-6;
     TestSTA(0,X0);
     %%%;
     SEDO(X0,0,2*%pi,1000);
     %%%;

Tout d'abord, l'EDO est définie comme une liste de la forme
'[VD,SM,IV]' ; VD pour Variables dépendantes, SM pour second membre,
IV pour invariants (mais on peut mettre n'importe quoi).  L'EDO est
alors
     %%%; nettoie()$ 'diff(VD,t) = SM %%%;
VD est la liste des symboles représentant les variables dépendantes
; SM celle des seconds membres qui dépendent des éléments de VD ;
IV ajoute des variables qui dépendent des variables dépendantes,
ainsi d'ailleurs que de la variable indépendante dont le symbole est
't'.

Cette structure EDO est alors injectée à 'Prepare' qui crée les
fonctions 'SecondMembre', 'Action', 'Tolerance' décrites
Utilisation plus poussée.

La fonction Tolerance nécessite de donner une valeur à la variable
'tol' : le pas sera validé si l'écart entre les solutions calculées
avec deux ordres différents est plus petite que 'tol' en norme 'max'.

On donne ensuite une condition initiale par 'X0=[1,0]'.  Puis on teste
un premier pas par 'TesStat(temps initial,condition initiale)'.  Le
retour doit être de la forme
     Second Membre :  [1,0,1]
     Tolerance -> true :  true
     Action :  [[0,1,0,0]]
la liste après 'Second Membre' est la valeur du second membre pour la
condition initiale ; elle doit être entièrement numérique ; après
'Tolerance -> true' la valeur doit être 'true' ou 'false', si c'est
'unknown' il ne faut pas lancer le calcul parce que le mécanisme
d'adaptation de pas ne fonctionnera pas.  La liste après 'Action' est
la valeur initiale de VD à laquelle est ajoutée en début de liste
celle du temp initial.

'TestStat' a aussi pour fonction d'initialiser la variable dans
laquelle seront rangées les valeurs de VD au cours du calcul.  Cette
variable est 'lX' (détail de sa structure plus bas).

Le calcul complet est alors lancé par *SEDO(X0,temps initial,temps
final, nombre de pas)* qui remplit 'lX'.

Avec une méthode à pas variable, on peut être étonné que le
'nombre de pas' soit fixé.  c'est utile pour qu'il y ait suffisamment
de points pour que les tracés de trajectoire semblent lisses (voir
ci-dessous).  'nombre de pas' est en fait le nombre minimal de pas de
temps qu'on impose.

Lorsque le calcul est fini, on peut examiner son résultat par des
tracés de trajectoire comme
     %%%; XvsX(EDO,[x,y])$ %%%;
ou de variables en fonction du temps comme
     %%%; XvsT(EDO,x,y)$ %%%;
Les fonctions 'XvsX' et 'XvsT' peuvent comporter des arguments en plus
grand nombre et comportant des calculs sur les variables comme
     %%%; XvsX(EDO,[x,y],[x,y+x],[x,y-x])$ %%%;
ou
     %%%; XvsT(EDO,x*t,y,x^2-y^2)$ %%%;

Finalement, pour illustrer ce qui a été dit sur 'nombre de pas', on
peut refaire un calcul où il a une valeur faible
     %%%;
     TestSTA(0,X0)$
     SEDO(X0,0,2*%pi,2)$
     XvsT(EDO,x,y)$
     lX;
     %%%;
On a demandé 2 pas de temps, il y en a eu 2.  En dépit de ce que les
x et y ne ressemblent pas à des sinusoïdes leurs valeurs en pi et 2
pi sont proches des valeurs de la solution exacte.  C'est parce qu'on
utilise la paire 5-4 de Dormand-Prince qui est précise.

L'examen de lX montre qu'elle est de la forme '[...,[temps,x au
temps,y au temps,nombre]...]'.  'nombre d'itérations' est (en gros,
c'est un peu plus compliqué) le nombre de pas intermédiaires qui ont
été fait sans que leurs valeurs aient été reportées.

5.2.2 Utilisation plus poussée
-------------------------------

La fonction *SEDO* appelle les fonctions *Action* et *Tolerance* pour
effectuer les pas de temps ; celle-ci sont créées par 'Prepare'.
Elle utilise par défaut la paire de Dormand-Prince.

Il est possible également de changer cela, en utilisant
*EDORKe(SecondMembre,X0,t0,t1,N,Action,Tolerance,Paire)* plutôt que
'SEDO'.

Ces fonctions peuvent être créées manuellement.

5.2.2.1 Paire choisie
.....................

La paire de moins haut ordre est la paire de Heun Euler correspondant
à la fonction *HeunEuler(SM,Var0,h)*.  On peut la tester sur
l'oscillateur harmonique
     %%%; nettoie()$ depends([x,y],t)$
     [diff(x,t)=-y,diff(y,t)=x]@
     %%%;
On écrit plutôt ce système sous forme autonome
     %%%; nettoie()$ depends([t,x,y],s)$
     [diff(t,s)=1,diff(x,s)=-y,diff(y,s)=x]@
     %%%;
et donc le second membre est
     %%%; nettoie()$
     matrix([1],[-y],[x])@
     %%%;
La fourniture de ce second membre à 'HeunEuler', avec une valeur
initiale '[0,1,1]' fournit
     %%%; nettoie()$
     SecondMembre(X):=[1,-X[3],X[2]]$
     expand(HeunEuler(SecondMembre,[0,1,1],h))@
     %%%;
Si l'intégration avait été analytique, on aurait obtenu
     %%%; nettoie()$
     taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,3)@
     %%%;
et donc par comparaison
     %%%; nettoie()$
     SecondMembre(X):=[1,-X[3],X[2]]$
     res:expand(HeunEuler(SecondMembre,[0,1,1],h))$
     res[1]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,3)@
     res[2]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,3)@
     %%%;
on voit que la première liste renvoyée 'res[1]' est d'ordre 2 alors
que la second est d'ordre 1.  Heun Euler est une méthode d'ordre 2 et
1.

La paire Bogacki Shampine est d'ordre 3 et 2
     %%%; nettoie()$
     SecondMembre(X):=[1,-X[3],X[2]]$
     res:expand(BogackiShampine(SecondMembre,[0,1,1],h))$
     res[1]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,4)@
     res[2]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,4)@
     %%%;
et celle de Dormand Prince 5 et 4
     %%%; nettoie()$
     SecondMembre(X):=[1,-X[3],X[2]]$
     res:expand(DormandPrince54(SecondMembre,[0,1,1],h))$
     res[1]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,6)@
     res[2]-taylor([h,cos(h)-sin(h),cos(h)+sin(h)],h,0,6)@
     %%%;

On peut ajouter d'autres paires, pourvu qu'elle renvoient une liste
composée de deux estimations de la solution qui soient d'ordres
différents.

5.2.2.2 Tolérance
..................

La fonction 'Tolerance' créée par 'Prepare' est très simple, c'est
     %%%;
     Tolerance(X,tX,p):=block([dX:X-tX,ndX],ndX:apply(max,abs(dX)),is(ndX < tol))$
     %%%;
Elle calcule la valeur maximale de la valeur absolue de la différence
entre les deux estimations fournis par la paire choisie et renvoie
'true' si cette valeur est plus petite que 'tol' la tolérance qui est
un paramètre à ajuster.

Ce n'est pas une très bonne façon de faire si les variables n'ont
pas été préalablement adimensionnée.  On peut adapter la fonction
'Tolerance' à cette situation.

5.2.2.3 Action
..............

Chaque fois que le temps (la variable indépendante) correspond (ou
dépasse un peu) les subdivisions de temps prédéfinies, la fonction
'Action' est appellée.

Celle que crée Prepare est
     %%%;
     Action(X,p):=lX:append(lX,[append(X,[p])])
     %%%;
La valeur de la solution au pas de temps atteint est ajouté à la
liste 'lX'.

On peut vouloir faire des choses plus amusantes comme des animations :
ici c'est le mouvement d'une masse assujettie à rester sur un cercle
dans un plan vertical.  La flêche visualise la force qui s'exerce sur
la masse (sans mettre la force de réaction du cercle sur la masse)
     %%%;
     [X,Xp,dX,dXp]:InitVars(x,z)$ data:[m=1,g=9.81,f=x*xp+z*zp]$
     L:1/2*m*(xp^2+zp^2)-m*g*z$ e:1$ C:[x*xp+z*zp/e^2]$ D:0$F:[]$
     SM:EquationLagrange(L,d,C,F,X,Xp);
     InitSM(X,Xp,SM,data,float)$
     InitTol(X,Xp,'block([dmx:1.0e-8],is(apply(max,abs(''dX)) < dmx)))$
     scl:1.5*max(e,1)$
     set_draw_defaults(xrange=[-scl,scl],yrange=[-scl,scl],head_length=0.05,head_angle=30,
                       fill_color=white,points_joined = true)$
     InitAct(X,Xp,'block([],
       ?sleep(0.2),draw2d(line_width=1,color=green,fill_color=white,ellipse(0,0,1,e,0,360),
       color=red,fill_color=red,ellipse(x,z,0.02,0.02,0,360),color=black,
       vector([x,z],0.05*((xp^2+zp^2)*[x,z]+[0,-g]))
     )),data)$
     X0:[0,e,0.01,0];
     TestSTA(0,X0);
     IntSM(X0,0,4.8,1000);
     %%%;
Les fonctions 'SecondMembre', 'Tolerance', 'Action' sont générées
par *InitSM*, *InitTol* et *InitAcT* mais c'est aussi facile de
fabriquer ces fonctions manuellement.

5.3 Lagrangien et Hamiltonien
=============================

5.3.1 Equations de Lagrange
---------------------------

*EquationLagrange(L,D,C,F,X,Xp)* génère les équations de Lagrange
pour un problème dont
   * le lagrangien est 'L'
   * la fonction de dissipation de Rayleigh 'D'
   * les contraintes (sous forme non-holonome même si elle sont
     holonomes) 'C'
   * les forces généralisées additionnelles 'F'
   * les variables de configuration 'X'
   * leurs vitesse 'Xp'

Par exemple, pour un système masse-ressort
     %%%; nettoie()$
     SM:EquationLagrange(1/2*m*v^2-1/2*k*x^2,0,[],[],[x],[v])$
     EDO:[[x,v],ev(SM)];
     Systeme(EDO)@
     %%%;
pour un système masse-ressort avec force appliquée et dissipation
     %%%; nettoie()$
     SM:EquationLagrange(1/2*m*v^2-1/2*k*x^2+f*x,alpha*v^2/2,[],[],[x],[v])$
     EDO:[[x,v],ev(SM)];
     Systeme(EDO)@
     %%%;
pour un système masse-ressort avec force appliquée et dissipation,
la force étant mise dans F plutôt que dans le lagrangien
     %%%; nettoie()$
     SM:EquationLagrange(1/2*m*v^2-1/2*k*x^2,alpha*v^2/2,[],[f],[x],[v])$
     EDO:[[x,v],ev(SM)];
     Systeme(EDO)@
     %%%;

Maintenant, un point matériel se déplaçant sur un cercle
     %%%; nettoie()$
     LH:x^2+y^2-R^2 LNH:derive(LH,[[x,u],[y,v]])@
     SM:EquationLagrange(1/2*m*(u^2+v^2),0,[LNH],[],[x,y],[u,v]);
     %%%;
La liaison holonome LH est transformée en non-holonome LNH par
dérivation Dérivation et c'est LNH qui est fournie à
EquationLagrange.  Le second membre a la forme d'un block
Functions and Variables for Function Definition (Maxima) pour ne pas
ranger directement l'expression du multiplicateur de Lagrange '%mu1'
dans le vecteur second membre.  Mais bien entendu on peut évaluer ce
block
     %%%; factor(ev(SM,x^2=R^2-y^2))@ %%%;
pour retrouver le résultat attendu (le cercle a pour rayon R).

Comme la liaison holonome a été changée en une liaison
non-holonome, il y a perte d'une constante (ici le rayon du cercle).
C'est la condition initiale qui permet de retrouver cette constante.
Ici si on décide qu'initialement 'x=R,y=0' et cela pours des valeurs
initiales de '(u,v)' quelconques, l'intégration de l'EDO se fera de
manière que 'x^2+y^2=R^2'.

5.3.2 Circuit électrique
-------------------------

Le circuit électrique relève a priori d'un traitement par
'EquationLagrange', les positions et vitesses étant remplacées par
les charges(1) et les courants électriques.

Neanmoins les circuits électriques comportent de nombreuses liaisons
(pour la loi aux nœuds) et les équations de Lagrange fournies sont
souvent dépourvues d'inertie (les variations temporelles de flux
magnétiques).  Ce qui justifie la confection d'une fonction
particulière.

   ---------- Footnotes ----------

   (1) quantité d'électricité dans le chapitre ad hoc du Treatise
of Electricity and Magnetism

5.3.2.1 EQCirc
..............

La première fonction utile au circuit électrique est
*EQCirc(Circuit)* où Circuit est une liste de la forme
'[...,[N,P,i,u,R],...]' où
   * N et P sont des symboles représentant des nœuds ; entre ces
     nœuds il y a une liaison ;
   * i est un symbole représentant le courant électrique circulant
     danbs la liaison de N à P ;
   * u est un symbole représenant la tension correspondant à la
     différence de potentiel entre les nœuds N et P ;
   * rel est une relation entre les tensions et les courants de tout
     le circuit ; mais en général cette relation se réduit à une
     relation entre i et u

Par exemple, pour le circuit rl avec source
     %%%; nettoie()$
     circ:[[A,A,i,u,u=l*ip+r*i-e]]$
     EDO:EQCirc(circ);
     Systeme(EDO)@
     %%%;

La sortie de EqCirc est '[[i],[-(i*r-e)/l],[u = 0],ver = []]' ; c'est
une liste comportant : la liste des variables dépendantes de l'EDO,
la liste des seconds membres, la liste des variables éliminées et
leur dépendance par rapport à celles qui ne le sont pas et un
dernier élément 'ver=[]' qui cesse d'avoir cette forme lorsqu'il y a
eu une erreur dans la recherche de l'EDO.

Pour le circuit rlc avec source
     %%%; nettoie()$
     circ:[[A,B,i,u,u=l*ip+r*i-e],[B,A,j,v,j=c*vp]]$
     EDO:EQCirc(circ)@
     Systeme(EDO)@
     %%%;

Pour des bobines couplées
     %%%; nettoie()$
     circ:[[A,A,i1,u1,u1=l1*i1p+m*i2p+r1*i1-e1],[B,B,i2,u2,u2=l2*i2p+m*i1p+r2*i2-e2]]$
     EDO:EQCirc(circ)@
     Systeme(EDO)@
     %%%;

Toujours pour des bobines couplées mais avec un effet capacitif
     %%%; nettoie()$
     depends([e],t)$
     circ:[[A,C,i1,u1,u1=l1*i1p+m*i2p],[C,A,i3,u3,u3=-e],
           [B,C,i2,u2,u2=l2*i2p+m*i1p],[C,B,i4,u4,u4=R*i4],
           [A,B,i5,u5,i5=c*u5p]]$
     EDO:EQCirc(circ);
     Systeme(EDO)@
     %%%;

5.3.2.2 ChV
...........

Les exemples de EQCirc (particulièrement celui des bobines
couplées avec effet capacitif où second membre comporte une
dérivation de 'e') mettent en évidence que le choix des variables
conservées n'est pas toujours heureux.

La fonction *ChV(EDO,VARS)* permet de remédier à cela.  'EDO' est
l'EDO obtenue à partir de 'EQCirc' et 'VARS' la liste des variables
désirées.

Pour le circuit rlc avec source, on obtient
     %%%; nettoie()$
     circ:[[A,B,i,u,u=l*ip+r*i-e],[B,A,j,v,j=c*vp]]$
     EDO:EQCirc(circ)@
     Systeme(EDO)@
     EDO:ChV(EDO,[v,i])@
     Systeme(EDO)@
     %%%;
C'est cosmétique.

Par contre, pour les bobines couplées mais avec un effet capacitif
     %%%; nettoie()$
     depends([e],t)$
     circ:[[A,C,i1,u1,u1=l1*i1p+m*i2p],[C,A,i3,u3,u3=-e],
           [B,C,i2,u2,u2=l2*i2p+m*i1p],[C,B,i4,u4,u4=R*i4],
           [A,B,i5,u5,i5=c*u5p]]$
     EDO:EQCirc(circ);
     EDO:ratsimp(ChV(EDO,[i1,i2,u5]));
     Systeme(EDO)@
     %%%;
on gagne que le terme source 'e' n'apparait plus sous forme d'une
dérivation temporelle.  Ce qui n'eût pas été commode dans le cas
où 'e' eût été une fonction discontinue du temps.

5.3.2.3 Régime sinusoïdal établi
...................................

La fonction *RegSin(Circuit)* traite le cas où les variables de
'Circuit' sont les amplitudes complexes (efficaces mais c'est sans
importance) d'un régime sinusoïdal établi

Pour les bobines couplées avec effet capacitif
     %%%; nettoie()$
     circ:[[A,C,i1,u1,u1=l1*i1*%i*omega+m*i2*%i*omega],[C,A,i3,u3,u3=-e],
           [B,C,i2,u2,u2=l2*i2*%i*omega+m*i1*%i*omega],[C,B,i4,u4,u4=R*i4],
           [A,B,i5,u5,i5=c*u5*%i*omega]]$
     res:ratsimp(RegSin(circ));
     %%%;

La fonction *EDOSin(EDO)* traite le cas où 'EDO' a été obtenu par
EQCirc et qu'on veut le régime sinusoïdal établi.

6 Fonctions de service
**********************

'nettoie()'
     Remet l'état interne de Maxima à celui du démarrage.  Les
     variables qui ont été affectées sont trouvées par inspection
     de _values_ dans Functions and Variables for Command Line (Maxima).
'symbolise(x,y,...)'
     Efface les valeurs auxquelles x,y,...  ont été affectées, ce
     qui a pour effet que chacun des x,y,...  reprend comme valeur son
     symbole lui même.  Identique à 'x:'x' (pour x).
'Int2Quag(expr)'
     Cherche les intégrales dans 'expr' et les évalue numériquement
     avec 'quad_qag' Functions and Variables for QUADPACK (Maxima)
'Int2Romberg(expr)'
     Cherche les intégrales dans 'expr' et les évalue numériquement
     avec 'romberg' Functions and Variables for romberg (Maxima)
'id(x)'
     l'application identique
'eval(expr,r1,r2,...)'
     Fait une succession de 'ev' Functions and Variables for Evaluation (Maxima)
     sur 'expr' avec les règles 'r1', 'r2',...

7 Exemples plus complets
************************

Les exemples donnés servent surtout à expliquer les arguments des
fonction de CDE. Voici maintenant quelques exemples qui ont pour objet
d'illustrer ce que CDE permet de faire.

7.1 Décomposition de Helmholtz d'un champ dipolaire
====================================================

Exemple 2 : la décomposition d'Helmholtz d'une densité dipolaire.

Dans les listes v=[v1,v2], v1 est la valeur pour |x|< 1 et v2 pour
|x|>1.

D'abord le potentiel newtonien et son gradient

     %%%; nettoie()$ addv(x)$
     Regles:[[Grad(x),Id],[tGrad(x),Id],[div(x),3],[rot(x),0],[grad(ps(x,x)),2*x]]@
     PrefDyade:true@
     phi:[1/2-nrm(x)^2/6,1/(3*nrm(x))]@
     gphi:tav(map(lambda([%p],grad(%p)),phi))@
     xtav(map(lambda([%p],div(%p)),gphi))@
     saut_de_la_composante_normale=block([Regles:[[ps(x,x),1]]],tav(gphi[1]-gphi[2]))@
     saut_de_la_composante_normale=xtav(block([Regles:[[ps(x,x),1]]],tav(gphi[1]-gphi[2])))@
     %%%;

Ensuite la décomposition d'Helmholtz de [d,0] où d est un vecteur
constant.  On a [1,0]*d=b-h avec b=rot(grad(phi) x d) et
h=grad(grad(phi) .  d)

     %%%;  addv(d)$declare(d,constant)$
     h:tav(map(lambda([%p],grad(ps(%p,d))),gphi))@
     b:tav(map(lambda([%p],rot(pv(%p,d))),gphi))@
     tav(b-h)@
     %%%;

On en tire que la solution de
     %%%; nettoie()$ addv(b,h,d)$
     b-h = POD(d,cond([[nrm(x)<1,1]],0))@
     %%%;
est
     %%%; nettoie()$ addv(b,h,d)$
     ListeTex(b=rot(pv(grad(phi),d)),h=grad(ps(d,grad(phi))))@
     %%%;
avec
     %%%; nettoie()$ addv(x)$
     phi=cond([[nrm(x)<1,1/2-nrm(x)^2/6]],1/(3*nrm(x)))@
     %%%;
le potentiel newtonien.

C'est un résultat de 'Foundation of potential theory, P. D. Kellogg,
Dover, 1954' qui s'étend d'ailleurs au cas des ellipsoïdes.

7.2 Trièdre de Frenet
======================

Les éléments du trièdre de Frenet
     %%%; nettoie()$
     addv(X)$primise(X)$secoise(X)$tercise(X)$
     %appif(Xt,ListeTex1)$%appif(Xt,liste_des_vecteurs)$
     lder:[[X,Xp],[Xp,Xs],[Xs,Xt]]$
     addv(tau,n,b)$
     Norme(expr):=tav(expr/nrm(expr))$
     X2TF:[tau=Norme(derive(X,lder))]@
     X2TF:append(X2TF,[n=Norme(derive(ev(tau,X2TF),lder))])@
     X2TF:append(X2TF,[b=tav(ev(pv(tau,n),X2TF))])@
     primise(tau,n,b)$
     X2dTF:map("=",[taup,np,bp],derive(ev([tau,n,b],X2TF),lder))$
     cTF(expr):=block([%c:map(lambda([%u],%u*factor(tav(ev(ps(expr,%u),X2dTF,X2TF)))),[tau,n,b])],
                apply("+",%c))$
     dTF:map(lambda([%u],%u=cTF(%u)),[taup,np,bp])$
     assume(k>0,ps(Xp,Xp)>0)$
     dTF:expand(tav(subst(k^2*nrm(Xp)^4,tav(nrm(pv(Xp,Xs))^2),dTF)))$
     dTF:expand(tav(subst(chi*nrm(Xp)^3*k^2,pm(Xp,Xs,Xt),dTF)))@
     Cf:[k=nrm(pv(Xp,Xs))/nrm(Xp)^2,chi=pm(Xp,Xs,Xt)/(k^2*nrm(Xp)^3)]$
     Cf:[Cf[1],ev(Cf[2],Cf[1])]@
     %%%;

Pour une hélice, par exemple,
     %%%;
     GenCart()$
     block([X,Xp,Xs,Xt],
      X:R*(cos(s)*kx+sin(s)*ky)+a*s*kz,
      Xp:diff(X,s),
      Xs:diff(Xp,s),
      Xt:diff(Xs,s),
      assume(R>0),
      factor(trigsimp(xtav(ev(Cf)))))@
     %%%;

7.3 Déformation d'un cylindre sous son poids
=============================================

On veut calculer la déformée d'un cylindre axisymétrique sous
l'action de sa propre masse (le cylindre a été fabrique dans
l'espace, on le ramène sur terre et on le pose sur l'une de ses
bases).

Comme le calcul des opérateur en coordonnées a un coût, on le fait
une fois pour toute (mais cela suppose de ne plus utiliser les
'nettoie()'.

     %%%; GenCyl()$ %%%:

Le calcul complet (penser à exécuter 'Gencyl()' au préalable)
     %%%;
     ok;
     assume(r>0)$
     Indicie([r,z],u)$primise(ur,uz,psi)$
     vu:ur*kr+uz*kz$
     depends([ur,uz,urp,uzp,psi,psip],[r,z])$
     gu:xtav(Grad(vu))$
     Epsilon:xtav((gu+Trans(gu))/2)$
     TrE:xtav(Id*Epsilon)$
     Sigma:xtav(lambda*TrE*Id+2*mu*Epsilon)$
     W:expand(xtav(Int3D(Sigma*Epsilon,D)))$
     Wp:derive(W,[[ur,urp],[uz,uzp]])$
     F:expand(xtav(Int3D(ps(kz,vu),D)))$
     Fp:derive(F,[[ur,urp],[uz,uzp]])$
     'W=W@
     'F=F@
     my2cl:[lambda=nu*E/(1-2*nu)/(1+nu),mu=E/2/(1+nu)]@

     FaitMaillage('nom=m,'densite=50,'precision=14,'reels=[R=1,H=1],
      'points=[p1=[0,0],p2=[R,0],p3=[R,H],p4=[0,H]],
      'bords=[[Gammab0,Rec,p1,p2],[Gammal,Rec,p2,p3],[Gammab1,Rec,p3,p4],[axe,Rec,p4,p1]],
      'visus=['bords,'maillage,'domaines,'info])$
     cvs:[[D,[]],[r,x],[z,y],[detJ,2*pi*x]]$

     FFDeclare(real,E=1,nu=0.40,fct=1.0,lambda=ev(lambda,my2cl),mu=ev(mu,my2cl),hpsi=1.0e+10)$
     FFDeclare(mP1,ur,uz,urp,uzp)$
     FFSolve(pb,[[ur,uz],[urp,uzp]],append([FtoFF(m,Wp,cvs),FtoFF(m,Fp,cvs)],[On(Gammab0,ur=0,uz=0)]))$
     FFPlot(ur,cmm="ur")$FFPlot(uz,cmm="uz")$

     FFmovemesh(m,mo,[x+fct*ur,y+fct*uz])$
     FFPlot(m,mo,cmm="m et sa deformee")$
     FFgo()$
     %%%;
On écrit les énergies élastique 'W' et potentiel 'F' ; on les fait
varier ; on injecte cette variation à FreeFem++ et on lui demande
d'en calculer le point stationnaire (ici le minimum).

7.4 Déplacement d'une bille dans un tube rempli de fluide Newtonien
====================================================================

Voici un exemple où une bille se déplace sur l'axe d'un tube rempli
de fluide newtonien.
     %%%; GenCyl()$ %%%;

     %%%;
     Indicie([r,z],v)$depends([vr,vz,P,psi],[r,z])$
     primise(psi)$primise(v)$Indicie([r,z],vp)$depends([vpr,vpz,Pp,psip],[r,z])$

     /* Préparation des fonctionnelles */
     vV:vr*kr+vz*kz$
     gV:xtav(Grad(vV))$
     Gammap:xtav((gV+Trans(gV))/2)$
     Q:expand(xtav(Int3D(2*eta*Gammap*Gammap,D)))$
     LQ:Q+expand(xtav(Int3D(ps(grad(P),vV),D)))$
     DVG:expand(xtav(Int3D(-P*div(vV)+epsilon*P^2,D)))@
     LQ:Q+DVG$
     LQp:expand(derive(LQ,[[vr,vpr],[vz,vpz],[P,Pp]]))@
     FPsi:expand(xtav(Int3D(nrm(vV-rot(psi*kt))^2,D)+Int2D(hpsi*diff(psi,z)^2,Gtube)))@
     FPsip:derive(FPsi,[[psi,psip]])@

     /* Calcul avec FreeFem++ */

     /* Maillage */
     FaitMaillage('nom='m,'densite=30,'precision=6,
     'reels=[R=0.2,H=3,Rp=3],
     'points='[p1='[R,0],p2=[H,0],p3='[H,Rp],p4=[0,Rp],p5=[0,R]],
     'bords=[[iaxe,'Rec,p1,p2],[ibinf,'Rec,p2,p3],[itube,'Rec,p3,p4],[isym,'Rec,p4,p5],[ibille,'Arc,[0,0],R,pi/2,0]],
     'visus=[bords,maillage,info])$
     /* Conversion des identifiants de domaines et bords */
     cvs:[[D,[]],[Gbille,[ibille]],[Ginf,[ibinf]],[Gsym,[isym]],[Gtube,[itube]],[r,y],[z,x],[detJ,2*pi*y]]$
     /* Déclarations */
     FFDeclare(real,V=1,eta=100,epsilon=0.0,hpsi=1.0e+10)$
     FFDeclare(mP2,vr,vz,vpr,vpz,psi,psip)$
     FFDeclare(mP1,P,Pp)$
     /* Résolution */
     FFSolve(pb,[[vr,vz,P],[vpr,vpz,Pp]],append(FtoFF(m,LQp,cvs),
      [On(ibille,vr=0,vz='V),On(itube,vr=0,vz=0),On(isym,vr=0,P=0),On(ibinf,vr=0,vz=0)]))$
     FFPlot('vr,cmm="vr")$ FFPlot('vz,cmm="vz")$ FFPlot('P,cmm="P")$
     FFMinMax(vr,vz,P)$
     FFDeclare(real,'Q=2*FtoFF(m,Q,cvs),FSto=6*'eta*'pi*'V*'R,QSto='V*FSto,
       'FP=2*FtoFF(m,[Int2D(-P*z/R,Gbille)],cvs)[1],
       'FV=2*FtoFF(m,[xtav(Int2D(eta*ps(Gammap*(r*kr+z*kz)/R,kz),Gbille))],cvs)[1]
       )$
     /* Sortie */
     FFCout(('FP+FV)/'FSto,'Q/('FSto*'V))$
     FFCout('FP/'FSto,'FV/'FSto)$
     FFSolve(pbpsi,[[psi],[psip]],append([FtoFF(m,FPsip,cvs)],[]))$
     FFPlot('psi,cmm="psi")$
     FFMinMax('psi)$

     /* Commande synchrone d'appel de FreeFem++ */
     FFgo()$
     %%%;

7.5 Force exercée par un aimant sur une bille
==============================================

7.5.1 Calcul de base en potentiel vecteur
-----------------------------------------

     %%%; GenCyl()$ %%%;

     %%%;
     addvV(a)$ primiseV(va)$ Indicie([air,bille,aimant],D)$ Indicie([0,r],mu)$
     W:Int3D(nrm(rot(va))^2/(2*mu0),Dair+Daimant)+Int3D(nrm(rot(va))^2/(2*mu0*mur),Dbille)+Int3D(m*ps(kz,rot(va)),Daimant)@
     W:EclateDomaine(W)@
     Wp:derive(W,[[va,vap]])@
     depends([a,ap],[r,z])$ primise(a)$
     WpFF:expand(xtav(ev(Wp,va=a*kt,vap=ap*kt)))@

     FaitMaillage('nom=mail,'densite=10,
      'reels=[e=0.1,R=0.5,Ra=1.0,La=1.0,Zxd=3.0,Zxg=3.0,Rx=4.0],
      'points=[p1=[-Zxg,0],p2=[-La,0],p3=[0,0],p4=p5-pR,p6=p5+pR,p5=[R+e,0],p7=[Zxd,0],
               p8=p7+pRx,p9=p1+pRx,p10=p2+pRa,p11=p3+pRa,
               pRa=[0,Ra],pRx=[0,Rx],pR=[R,0]],
      'bords=[[iaxes,Rec,p1,p2],[iaxes,Rec,p2,p3],[iaxes,Rec,p3,p4],[iaxes,Rec,p4,p6],[iaxes,Rec,p6,p7],
              [ibx,Rec,p7,p8],[ibx,Rec,p8,p9],[ibx,Rec,p9,p1],
              [ibaN,Rec,p3,p11],[ibaL,Rec,p11,p10],[ibaS,Rec,p10,p2],
              [ibC,Arc,p5,R,0,pi]],
      'epidemie=[[Dair,ibx],[Dbille,ibC],[Daimant,ibaL]],
      'visus=['bords,maillage,info,domaines])$
     cvs:[[z,x],[r,y],[detJ,2*pi*y]]$
     FFDeclare(mailP1,a,ap)$
     FFDeclare(real,mu0=1,mur=100,m=1)$
     FFSolve(pb,[a,ap],FtoFF(mail,WpFF,cvs),On(ibx,a=0))$
     FFDeclare(mailP1,ra=y*a)$
     FFPlot(ra,cmm="ra")$
     FFgo()$
     %%%;

La fonctionnelle à stationnariser est la somme de plusieurs
intégrales dont l'intégrant est différent des unes aux autres.
Cela se traduit par plusieurs sous-domaines dans le maillage.

8 Index
*******

Concepts
--------

Fonctions et variables
----------------------