<'a occupé quelques heures entre fin août 2014 et mi septembre 2014.
> Vous pouvez télécharger tous ces documents : [première archive, sans graphique](Modélisation_PDE_09_2014.zip), [seconde, avec graphiques](Modélisation_PDE_09_2014_avec_graphiques.zip) et [la dernière, avec un petit film en bonus](Modélisation_PDE_09_2014_avec_films.zip).
Et [aussi cette page en PDF](index.pdf) (compilé avec [strapdown2pdf](http://lbesson.bitbucket.io/md/strapdown2pdf.html)).
## Présentation du problème
### D'où vient le problème ?
Du monde réel ! Il s'agit d'un modèle mathématique ayant pour but de rendre compte d'une observation réellement observée lors d'une expérience par un ami, l'expérimentateur [Nicolas Vernier](http://www.ief.u-psud.fr/?page_id=2868), docteur en physique et maître de conférence à l'Université Paris-Sud, qui était en mission à la [Mahindra École Centrale](http://www.mahindraecolecentrale.edu.in/contact.html) en juillet, août et septembre avec [moi](http://besson.qc.to/) ([Lilian Besson](https://bitbucket.org/lbesson/)).
Le problème consiste à étudier le déplacement d'un front de matière (un semi-conducteur dilué magnétique de type `GaMnAs`), soumis à un courant électrique mono-directionnel.
On s'attendrait à voir le front de matière se déplacer d'un bloc dans la direction du courant électrique, mais les observations réalisées en laboratoire montrent un tout autre comportement, que nous allons justifier *mathématiquement* ici.
> Si possible, nous détaillerons ici un peu plus le contexte physiquesous-jacent, ainsi que les calculs et hypothèses qui permettent de proposer le modèle mathématique.
### Précisons les notations
Dans ce document, $y$ désigne une fonction de classe $\mathcal{C}^1$ à deux variables : une variable d'espace $x$, et une variable de temps $t$.
$k > 0$ est une constante, dont la valeur réelle n'importe peu.
On demande à $x$ de prendre ces valeurs dans l'intervalle $[0, 1]$, et le temps $t$ vérifie seulement $t \ge 0$.
### Et le profil initial
À l'instant initial, $t=0$, pour coller aux conditions de l'expérience, on demande à $y$ de suivre un *profil parabolique*, de la forme $y(x, 0) = \beta \times (1 - \left(x-\frac{1}{2}\right)^2)$, où $\beta>0$ est une certaine *constante* (encore une fois, sa valeur physique n'a pas beaucoup d'importance).
### L'équation vérifiée par $y$
Enfin, et c'est la partie intéressante, $y$ satisfait l'équation aux dérivées partielles non-linéaire suivante :
$$ \left(\frac{\partial y} {\partial x}\right)^2 = \frac{k}{\left(\frac{\partial y} {\partial t}\right)^2} - 1. \;\;\; (1)$$
Une autre façon d'écrire cette équation est $ \left(\left(\frac{\partial y} {\partial x}\right)^2 + 1\right) \left(\frac{\partial y} {\partial t}\right)^2 = k.$
Ou encore $ \left((y\_x)^2 + 1\right) \left(y\_t\right)^2 = k.$
### Le problème à résoudre
- L'objectif est de trouver la (ou *une*) solution à cette EDP, afin de confronter les relevés expérimentaux avec le modèle mathématique.
- Un objectif plus raisonnable est d'afficher l'évolution de ce $y$ au cours du temps, même si ce n'est que numériquement, cela nous suffira pour vérifier la validité du modèle.
---
## Première tentative de résolution
### Essais à la main
Assez naturellement, la première chose qu'on tente est une approche directe, à la main.
- Par séparation des variables : échec !
- Par intégration directe : échec !
### Sur *StackExchange.com* ?
Ensuite, ne sachant pas vraiment comment partir, j'ai demandé de l'aide, pour tester le site [math.stackexchange.com](http://math.stackexchange.com).
J'ai posé [une question](http://math.stackexchange.com/q/897299/), que voici (en anglais) :
> For this problem, $y$ is a function of two variables: one space variable $x$ and one time variable $t$.
$k > 0$ is some constant.
And $x$ takes is value in the interval $[0, 1]$ and $t \ge 0$.
At the initial time, $y$ follows a **parabolic profile**, like $y(x, 0) = 1 - (x-\frac{1 }{2})^2$.
> Finally, $y$ satisfies this PDE:
$$ \left(\frac{\partial y} {\partial x}\right)^2 = \frac{k}{\left(\frac{\partial y} {\partial t}\right)^2} - 1.$$
> Does anyone have an idea how to solve this problem (and find the expression of $y(x,t)$) ?
### Une réussite *partielle* pour un cas particulier
#### Via Charpit-Lagrange
[En utilisant la](http://math.stackexchange.com/a/897605/169728) [méthode de Charpit-Lagrange (résumée ici)](http://www.sci.brooklyn.cuny.edu/~mate/misc/partial_diffeqs/charpits_method_compl_int.pdf), on arrive à montrer que le résultat suivant :
> Si la condition initiale suit un *profil affine* (*ie* $\forall x\in[0,1], y(x,0) = ax + b$, pour deux $a,b\in\mathbb{R}$), alors :
> - il y a exactement deux solutions, qui se propage au cours du temps à vitesse constante mais en sens inverse,
> - qui s'exprime par $\forall x\in[0,1],t>0, y\_{1,2}(x,t) = ax + b \pm \frac{\sqrt{k}}{\sqrt{1+a^2}} t$.
Mais ce résultat ne s'applique pas avec la condition initiale qui nous intéresse !
Néanmoins, on garde cet exemple en tête, car il permettra d'avoir une certaine confiance dans le schéma numérique développé ensuite.
#### Méthode des caractéristiques
En développant la méthode des caractéristiques, comme expliquée sur [cet article Wikipédia](https://en.wikipedia.org/wiki/Method_of_characteristics#Example), on arrive à montrer que si $y$ est solution, alors elle vérifie :
$$ (y\_x+2x-1)(y\_x^2+1)^\frac{3}{2}=\pm2\sqrt kty\_x \;\;\; (2)$$
Cette forme a l'avantage d'avoir fait disparaître $y\_t$ et ainsi $(2)$ est implicite en la seule inconnue $y\_x$, mais étant polynomiale de degré $6$, sa résolution n'est pas facile. En pratique, [Maple][] va échouer.
---
## Une seconde *tentative* de résolution
### De façon formelle, avec [Maple][]
Dans [cette feuille de calcul Maple v12](pde_09_2014.mw), on essaie d'utiliser les différents moyens que [Maple][] propose pour résoudre *formellement* une EDP.
En bref, c'est encore un échec. Pour plus de détails, veuillez ouvrir la veuille `.mw`, qui contient quelques commentaires.
[Maple]: http://www.maplesoft.com/products/maple/ "Célèbre logiciel de calcul formel"
### De façon numérique, avec [Maple][]
Après l'échec prévisible d'une résolution *magique* formelle, j'ai aussi demander à [Maple][] de le traiter *numériquement*.
Second échec.
Même en utilisant la forme "simplifiée" $(2)$, je n'ai pas réussi à tirer grand chose de [Maple][], dommage.
---
## Quelques résultats théoriques
Nous avons réussi à prouver les différents résultats, qui n'ont pas grand intérêts.
- Conservation d'une *symétrie* initiale en $x$. Si on cherche à résoudre pour $x\in[-1,1]$, et si $x \mapsto y(x,t)$ est paire, elle le restera pour tout $t>0$,
- Le **Hessien** est toujours nul ! $$\mathbb{H}(y) (x,t) = \mathrm{det} \left[\begin{matrix} y\_{xx}(x,t) & y\_{xt}(x,t) \\\\
y\_{tx}(x,t) & y\_{tt}(x,t) \end{matrix}\right] = y\_{xx}(x,t) \times y\_{tt}(x,t) - y\_{xt}(x,t) \times y\_{tx}(x,t) = 0, (\forall x \in [0,1], \forall t \geq 0).$$
Bien que j'ai un souvenir d'une conséquence *intéressante* de la nullité du hessien, je n'ai pas réussi à le retrouver. Toute idée est la bienvenue !
- Il n'y a pas unicité de la solution. Du fait du carré sur le terme $\frac{\partial y}{\partial t}$ dans l'EDP $(1)$, il y aura toujours deux solutions distinctes, l'une se propagant "vers le haut" au cours du temps ($\frac{\partial y}{\partial t} > 0$), l'autre "vers le bas" (comme ce qu'on observe dans le cas particulier du profil initial affine).
---
## Une dernière approche, numérique « à la main »
### Notation pour le maillage
Après l'échec relatif des précédentes approches, je me suis rappelé de [mon stage de Licence](http://besson.qc.to/publis/rapportL3Maths12.pdf), et j'ai pu mettre au point un petit schéma numérique d'ordre 1, tout simple, que voici.
- On discrétise l'espace $[0,1]$ en $N>1$ points, $x\_i = i \times h$, pour un pas d'espace ("assez petit") : $h = \frac{1}{N+1}$ et pour chaque $0 \leq i \leq N$.
- On fait de même pour le temps, en écrivant $t_j = j \times \delta$, pour un certain pas de temps ("assez petit") : $\delta>0$, et pour tout $j \geq 0$.
### Intérêt du maillage
Ensuite, on dit que pour $i,j$, on pose $\overline{y}\_{i,j}$ comme étant la moyenne de la (d'une) solution $y$ sur le domaine $\mathcal{A}\_{i,j} = [ x\_i - \frac{h}{2}, x\_i + \frac{h}{2} ] \times [ t\_j, t\_j + \delta ]$.
C'est-à-dire qu'on pose $$\overline{y}\_{i,j} = \frac{1}{h\delta} \int\int\_{\mathcal{A}\_{i,j}} y(x,t) \mathrm{d}t \mathrm{d}x = \frac{1}{h\delta} \int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \int\_{t\_j}^{t\_j + \delta} y(x,t) \mathrm{d}t \mathrm{d}x.$$
### On intègre l'EDP $(1)$
Sur un volume de contrôle $\mathcal{A}\_{i,j}$, on intègre l'équation $(1)$, ce qui nous donne :
$$\int\int\_{\mathcal{A}\_{i,j}} \left(\left(\frac{\partial y} {\partial x}\right)^2 (x,t) + 1\right) \left(\frac{\partial y} {\partial t}\right)^2 (x,t) \mathrm{d}t \mathrm{d}x = k \times \mathrm{Vol}(\mathcal{A}\_{i,j}).$$
> Attention, les calculs suivants sont **complètement faux !**
Soit $$\left(\int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \int\_{t\_j}^{t\_j + \delta} \left(\left(\frac{\partial y} {\partial x}\right)^2 (x,t) + 1\right)\mathrm{d}t \mathrm{d}x\right) \times \left(\int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \int\_{t\_j}^{t\_j + \delta} \left(\frac{\partial y} {\partial t}\right)^2 (x,t) \mathrm{d}t \mathrm{d}x\right) = k \times h \times \delta.$$
Mais dès lors on peut simplifier un peu, avec deux appels au théorème de Fubini sur un compact pour une fonction continue, car le caractère $\mathcal{C}^1$ qu'on impose à la solution $y$ donne la continuité des dérivées partielles :
$$ \left( \int\_{t\_j}^{t\_j + \delta} \left( \left( \int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \frac{\partial y} {\partial x}(x,t) \mathrm{d}x \right)^2 \right) + \delta\; \mathrm{d}t \right) \times \left( \int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \left( \int\_{t\_j}^{t\_j + \delta} \frac{\partial y} {\partial t}(x,t) \mathrm{d}t \right)^2 \mathrm{d}x \right)= k \times h \times \delta.$$
Soit $$ \left( \int\_{t\_j}^{t\_j + \delta} \left( y(x\_i + \frac{h}{2}, t) - y(x\_i - \frac{h}{2}, t) \right)^2 + \delta\; \mathrm{d}t \right) \times \left( \int\_{x\_i - \frac{h}{2}}^{x\_i + \frac{h}{2}} \left( y(x, t\_{j+1}) - y(x, t\_j) \right)^2 \mathrm{d}x \right) = k \times h \times \delta.$$
> Bref, ce calcul est **complètement foireux** : on ne pas intégrer un produit en un produit d'intégrale !
### Expression des dérivées
L'approche consiste ensuite à approcher $y(x,t)$ pour n'importe quel $x \in[0,1]$ et $t \geq 0$ par la valeur numérique $\overline{y}(x\_i, t\_j)$ pour le $i,j$ qui convient (tel que $(x,t)\in\mathcal{A}\_{i,j}$). Ces notations seront conservées dans la suite du paragraphe.
Le choix réside dans la formule utilisée pour approcher les deux dérives partielles $\frac{\partial y}{\partial x}(x,t)$ et $\frac{\partial y}{\partial t}(x,t)$.
Le schéma que j'ai conçu (plus réutilisé que conçu d'ailleurs) est un schéma simple d'ordre 1, ainsi le choix pour la dérivée en temps est assez immédiat :
$$ \frac{\partial y}{\partial t}(x,t) \simeq \frac{\overline{y}(x\_{i}, t\_{j+1}) - \overline{y}(x\_i, t\_j)}{\delta} $$
### Différents schémas pour la dérivée en espace
#### Gauche
Le plus intuitif est un schéma dit schéma *gauche*.
$$ \frac{\partial y}{\partial x}(x,t) \simeq \frac{\overline{y}(x\_{i+1}, t\_j) - \overline{y}(x\_i, t\_j)}{h}. $$
Il présente un risque lorsque $i=N$, car $x\_{i+1}$ n'est plus bien défini.
Il faut alors utiliser une connaissance des *conditions de bords* à droite.
Dans notre situation, j'ai simplement fait en sorte qu'on calcule toujours une valeur de plus, $(\overline{y}(x\_{N+1}, t\_j))\_{j \geq 0}$, même si elle ne sera pas affichée.
#### Droit
Le second schéma est dit schéma *droit*.
$$ \frac{\partial y}{\partial x}(x,t) \simeq \frac{\overline{y}(x\_i, t\_j) - \overline{y}(x\_{i-1}, t\_j)}{h}. $$
Il présente un risque lorsque $i=0$, car $x\_{i-1}$ n'est plus bien défini.
Il faut alors utiliser une connaissance des *conditions de bords* à gauche.
Dans notre situation, j'ai fait l'hypothèse que $\overline{y}(x\_{-1}, t\_j)$ était toujours égal à $\overline{y}(x\_{0}, t\_j)$, ce qui a pour conséquence d'observer au cours du temps une translation pour la valeur $\overline{y}(x\_{0}, t\_j)$
#### Centré
Le dernier schéma est dit *schéma* centré.
$$ \frac{\partial y}{\partial x}(x,t) \simeq \frac{\overline{y}(x\_{i+1}, t\_j) - \overline{y}(x\_{i-1}, t\_j)}{2h}. $$
Ce dernier schéma réunit les deux risques des schéma gauche et droite. Il suffit d'utiliser les deux techniques précédemment expliquées !
### Quelles hypothèses sur $h$ et $\delta$ ?
On comprend dès lors que les hypothèses "$h$ assez petit" et "$\delta$ assez petit" prennent tout leur sens ici.
J'aimerais trouver une condition **CFL** pour garantir la stabilité du schéma, mais sans succès pour l'instant.
### Du schéma implicite... à l'explicite
Les calculs suivants utilisent le schéma gauche.
Si on injecte les expressions des dérivées partielles approchées dans $(1)$, on arrive à :
$$ \left(\frac{\overline{y}(x\_{i}, t\_{j+1}) - \overline{y}(x\_i, t\_j)}{\delta}\right)^2 . \left( 1 + \left( \frac{\overline{y}(x\_{i+1}, t\_j) - \overline{y}(x\_i, t\_j)}{h} \right)^2 \right) = k. \;\;\; (3)$$
Notre but est de calculer itérativement le vecteur $\overline{y}$ au cours du temps, soit donc de pouvoir exprimer les valeurs de $\left(\overline{y}(x\_{i}, t\_{j+1})\right)\_{0\leq i \leq N}$ en fonction des valeurs de $\left(\overline{y}(x\_{i}, t\_{j})\right)\_{0\leq i \leq N}$.
En transformant $(3)$, on arrive très vite à : $$\overline{y}(x\_{i}, t\_{j+1}) = \overline{y}(x\_{i}, t\_{j}) \pm \frac{\delta \sqrt{k}}{\sqrt{1+\left( \frac{\overline{y}(x\_{i+1}, t\_j) - \overline{y}(x\_i, t\_j)}{h} \right)^2}}.$$
On retrouve les deux solutions duales, une qui sera croissante (le front de matière va vers le haut), une décroissante (vers le bas). La seule qui soit physiquement acceptable est celle avec un $+$ en lieu et place du $\pm$.
On parvient donc à un schéma implicite d'ordre $1$, qui nécessite de connaître le vecteur entier $\left(\overline{y}(x\_{i}, t\_{j})\right)\_{0\leq i \leq N+1}$ au temps $t\_j$ pour le calculer au temps suivant :
$$\overline{y}(x\_{i}, t\_{j+1}) = \overline{y}(x\_{i}, t\_{j}) + \frac{\delta \sqrt{k}}{\sqrt{1+\left( \frac{\overline{y}(x\_{i+1}, t\_j) - \overline{y}(x\_i, t\_j)}{h} \right)^2}}. \;\;\; (4)$$
> On vérifie par analyse dimensionnelle, ouf ça a l'air bon !
Dès lors, calculer $(4)$ itérativement va être plutôt facile.
### Transformer cela en un calcul itératif
On choisit un pas d'espace *"assez petit"*, une valeur numérique à $k$.
Il faut commencer par initialiser le vecteur $\overline{y}$ à l'instant initial, grâce au profil de $y(x,0)$.
Ensuite, on règle numériquement le pas de temps, expérimentalement en le prenant assez petit pour les calculs restent stables, suffisamment longtemps pour observer un comportement similaire à celui réellement observé.
---
### Implémentation du la simulation
#### En [OCaml](http://ocaml.ocamlpro.com/) et avec [GNUPlot](http://www.gnuplot.info/)
Voir [la documentation](./doc/Pde_09_2014.html) générée via [OCamlDoc](http://caml.inria.fr/pub/docs/manual-ocaml/ocamldoc.html), ou [le code source de pde.ml](./doc/pde_09_2014.ml.html) coloré par [caml2html](http://mjambon.com/caml2html.html).
Voici un exemple du script [plotdata.gnuplot](./plotdata.gnuplot) tel que généré étape par étape par le [petit programme OCaml](./doc/pde.ml.html), et le [fichier de données](./data.txt) (valeurs de $(x\_i)\_{0 \leq i \leq N}$ et $(\overline{y}\_{i,j})\_{0 \leq i \leq N}$ à un instant $t\_j$).
#### En [Octave](http://www.gnu.org/software/octave/) (ou Matlab)
Le [script suivant](pde_09_2014.m) est plus simple à lire et à comprendre, et fonctionne plus simplement.
J'ai progressivement rafiné la simulation, en autorisant plus de paramètres. Le [code est assez clair et plein de commentaires](pde_09_2014.m).
> Mais j'avais envie de tester OCaml + GNUPlot d'abord !
### Résultat des simulations ?
La simulation correspond en gros à ce qu'on voulait observer. Mais c'est instable, complètement.
[Ce dossier montre des images successives d'une simulation en Octave](gauche/), utilisant les paramètres suivants :
$k = 1, N = 10000, dt = 0.00005, jMax = 10000, \beta = 1$.
De plus, pour l'affichage, on ne le fait qu'une fois toute les 1000 étapes, et on pré-processe la sortie en réalisant une moyenne locale sur une fenêtre de 50 cases (ça affine l'apparence).
#### Quelques étapes
> 
Au début, affichage symétrique du profil initial (parabolique).
> 
> 
> 
Au bout de 1000, 4000 et 6000 étapes.
Dès la troisième image, on voit apparaître ce qu'on observe physiquement : un pic commence à pointer au milieu.
> 
Dernière étape (ici au bout de 9000 actualisation de $\overline{y}_{j}$).
On voit vraiment le pic avancer au milieu. Le schéma est trop instable pour que la simulation foncionne bien très longtemps.
#### Aussi un film ?
[Voir ce film .avi qui résume l'autre simulation en OCaml](Simulation_Numérique__PDE_09_2014__Schéma_gauche.avi).
---
En espérant que ça puisse aider, au moins [Nicolas Vernier](http://www.ief.u-psud.fr/?page_id=2868).
> Je m'arrête là, plus vraiment le temps d'y revenir. Terminé le mardi 23 septembre 2014.
]]>