Résolution numérique d’une équation aux dérivées partielles

Bonjour !

Dans le cadre de mon TIPE, je m’intéresse à la modélisation des phénomènes de transport dans une dune de sable. Après plusieurs recherches, je suis tombé sur une équation régissant les phénomènes de reptation et de saltation¹. Cette équation est² $\partial_t h + A \left[\frac{B + \partial_x h}{\sqrt{1 + (\partial_x h)^2}}\right]_{x - \ell}^x = 0$ où $h \colon (x, t) \in \mathbb R \times \left[0, +\infty\right[ \longmapsto h(x, t)$ est la fonction solution et où $A$, $B$ et $\ell$ sont trois réels positifs. Pour remettre en contexte, cette équation s’applique dans un milieu granulaire (typiquement une dune de sable) où un vent arrive à un angle constant sur le sol et elle propose de modéliser les rides de sable observées. La quantité $h(x, t)$ est la hauteur de sable à la position $x$ et à l’instant $t$.

Afin de vérifier si elle décrit bien le phénomène réel, je voudrais la simuler numériquement. J’ai essayé de mettre en place un schéma d’Euler explicite, mais le résultat n’est pas concluant. De plus, je ne sais pas comment gérer la différence entre crochet au niveau des bords. En effet, pour $x = 0$, il va falloir évaluer $\partial_x h$ en $x - \ell = -\ell$ et cette quantité n’existera pas. J’ai contacté l’auteur et il n’a jamais simulé l’équation. De plus, aucune condition aux bords n’est précisée.

D’où mes questions. Avez-vous déjà fait face à une telle équation ? Et comment peut-on la simuler, i. e. comment peut-on gérer le terme entre crochet au niveau des bords ? Je précise que $\ell$ est normalement assez petite (de l’ordre de quelques diamètres d’un grain).

Merci beaucoup ! <3

Euler explicite pour une EDP ? Pour la résolution des EDP il y a 3 grandes méthodes :

• les différences finies,
• les éléments finis,
• les volumes finis

Je ne suis pas sûr qu’il soit très accessible à ton niveau mais tu peux essayer de voir pour utiliser freefem pour résoudre ton équation. Le faire toi même sera encore plus compliqué.

Oui Euler explicite pour des ÉDP, c’est pas l’idéal. Je vais essayer de me documenter sur ces méthodes pour voir si je peux en faire quelque chose. Sinon, je ne connaissais pas FreeFem, donc pourquoi pas. Merci pour ton aide !

Euler explicite pour une EDP ?

J’imagine qu’il parle de l’intégration en temps, et c’est pas forcément un problème si on fait gaffe à la stabilité.

FreeFEM est une galère à utiliser si tu ne connais pas un minimum les éléments finis (d’ailleurs bonjour la galère pour écrire la forme faible du problème), par contre implémenter un schéma aux différences finis est très simple.

@Teguad : tu peux nous écrire ton schéma numérique et nous expliciter ce que tu as pris comme conditions aux limites ?

EDIT : l’autre option est de laisser tomber cette équation et de simuler un tas de sable avec des éléments discrets.

C’est sûrement pas bien fait, mais voici comment j’ai mis en place la méthode explicite. Un peu au hasard, j’ai pris les conditions de Neumann aux bords. J’ai d’abord discrétisé l’espace puis le temps. Comme je ne savais pas comment gérer le terme entre crochet, j’ai pris $\ell$ (notée $\ell_\text r$ dans mon texte) égale au pas de l’espace et j’ai bidouillé ensuite avec les conditions aux bords. Je te mets ci-dessous ce que j’ai fait.

J’ai déjà simulé la formation d’un tas de simple grâce aux éléments distincts.

OK, ça m’a pas l’air absurde. Deux conseils :

• fais des différences finies centrées pour la dérivée spatiale $\partial_xh\approx\dfrac{h_{i+1}-h_{i-1}}{2\Delta x}$, tu gagnes un ordre de précision sans calculer plus et ça évite de biaiser ton estimation dans une direction ou l’autre ;
• ne prend pas $\ell=\Delta x$, $\ell$ est un paramètre physique en soit et tu n’as pas intérêt à l’avoir dépendant du pas que tu as choisis.

Par ailleurs, Neumann comme condition au bord me parait pas forcément judicieux, ça t’empêche d’avoir un flux de sable. Tu pourrais imposer un niveau, ou alors éventuellement un flux non nul qui dépend du vent.

Merci pour ces conseils ! Pour ton deuxième point, le truc c’est que justement je ne sais pas gérer l’évaluation de $\partial_x h$ au point $(x - \ell, t)$… L’avantage de prendre $\ell = \Delta x$ + Neumann était que je pouvais bien gérer aux bords. Si je prends un $\ell$ quelconque, comme je l’ai dit, pour $x = 0$, il va falloir que j’évalue en $x - \ell$ ce qui va peut-être sortir de l’espace. T’as une idée pour gérer ça ? Je sais pas si je suis très clair.

Je vais essayer ces différentes conditions aux bords. En modifiant vite fait, j’ai à peu près le même résultat qu’avant, mais il faut sûrement que je revoie mon implémentation en Python.

Merci beaucoup pour ces solutions. Je n’avais pas du tout pensé à circulariser l’espace. Pour ton deuxième point, c’est ce que je voulais un peu faire avec les conditions de Neumann. Je vais essayer de mettre en place les deux solutions et voir laquelle est la meilleure.

Pour ce qui est de l’évaluation de $\partial_x h$, je vais tenter de faire les DL moi-même et, si je n’y arrive pas, je te demanderai (ce qui sera fort probable). Encore merci ! Je te tiens au courant de mes avancées.

Tu auras besoin de trois points d’ailleurs, pas deux. Sinon tu vas perdre en précision.

Salut ! J’ai un peu commencé à mettre en place ce que tu m’as dit, mais je bloque au moment de l’approximation de $\partial_x h$ en $x - \ell$. Je n’ai réussi qu’à trouver deux développements de Taylor, mais pas le troisième. En l’état, je n’arrive pas à faire une combinaison linéaire pour me débarrasser de la dérivée seconde. De plus, si on s’arrête au premier degré, j’ai un truc, mais je ne pense pas que ça fonctionne. Je te met ci-dessous ce que j’ai fait. Merci de ton aide !

Oups désolé pour le temps de réponse. Merci, je vais essayer de me débrouiller avec ça. C’est ce que je voulais faire, en disant que mon $x_{j - p}$ était le plus proche de $x_j - \ell$. Mais en effet, ça peut être $x_{j - p - 1}$. Je vais redéfinir mon $\varepsilon$.

Ah oui, c’est vrai qu’on a même un $O(\Delta x^2)$, je n’avais pas vu, mais ça ne changera pas grand-chose à mon implémentation.

Edit : Je n’arrive toujours pas à faire des combinaisons linéaires pour obtenir $\partial_x h(x - \ell)$. Y’a un moyen simple de l’obtenir ? Est-ce qu’il faut négliger des termes à un moment ?

Edit 2 : Je suis sur une piste.

Y’a un moyen simple de l’obtenir ?

Oui. Je te mets une aide en secret si tu coinces.

Merci beaucoup ! Je viens de réussir à l’exprimer et je pense que c’est la bonne relation car ça marche pour $\varepsilon = 0$ et $\varepsilon = \pm\frac12\Delta x$. Je vais essayer d’implémenter ça en Python. Et si ça marche pas, je passerai à la version implicite.

Hmm, ça me parait bizarre que ton expression soit valide à la fois pour $\varepsilon=\frac12\delta x$ et $\varepsilon=-\frac12\delta x$. Ces cas là, bien que correspondant à la même configuration, sont en fait différent en pratique (le nœud $x_k$ sera d’un côté ou de l’autre), c’est pour ça qu’il est important de définir si l’intervalle de valeurs que peut prendre $\varepsilon$ est ouvert côté $\frac12\delta x$ ou côté $-\frac12\delta x$. J’ai du mal à voir comment tu peux obtenir une expression valide sans que l’on ait à se soucier de la borne que l’on autorise comme valeur possible.

Au temps pour moi, on obtient pas le bon résultat pour $\varepsilon = -\frac12\Delta x$. Voici donc la relation que j’obtiens : $\partial_x h(x_i - \ell, t) = \frac{1}{2\delta}\left(h(x_{k + 1}, t) - \frac{2\varepsilon}{\delta}[h(x_{k + 1}, t) - 2h(x_k, t) + h(x_{k - 1}, t)] - h(x_{k - 1}, t)\right) + O(\delta).$ (J’ai noté $\delta = \Delta x$ pour simplifier.) Je crois que c’est la bonne !

J’ai réussi à implémenter le schéma explicite et les résultats ne sont pas trop mal. Je vais regarder plus en détail le problème physique pour trouver de meilleures conditions aux bords (mon prof me conseille les flux constants aux bords). Je vais voir aussi pour le schéma implicite, mais le calcule de la jacobienne me paraît lourd. Du coup, problème résolu ! Merci @adri1 d’avoir pris le temps de m’aider !!!