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

Dont les dérivées sont évaluées en plusieurs points

Le problème exposé dans ce sujet a été résolu.

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 saltation1. Cette équation est2 th+A[B+xh1+(xh)2]xx=0\partial_t h + A \left[\frac{B + \partial_x h}{\sqrt{1 + (\partial_x h)^2}}\right]_{x - \ell}^x = 0h:(x,t)R×[0,+[h(x,t)h \colon (x, t) \in \mathbb R \times \left[0, +\infty\right[ \longmapsto h(x, t) est la fonction solution et où AA, BB 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)h(x, t) est la hauteur de sable à la position xx et à l’instant tt.

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=0x = 0, il va falloir évaluer xh\partial_x h en x=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


  1. Cf. https://tel.archives-ouvertes.fr/tel-00199032, page 143

  2. On a noté [f(x,t)]xx=f(x,t)f(x,t)[f(x, t)]_{x - \ell}^x = f(x, t) - f(x - \ell, t).

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é.

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.

+1 -0

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 r\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.

Schéma
Schéma

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

+0 -0

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

  • fais des différences finies centrées pour la dérivée spatiale xhhi+1hi12Δx\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 =Δx\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.

+1 -0

Merci pour ces conseils ! Pour ton deuxième point, le truc c’est que justement je ne sais pas gérer l’évaluation de xh\partial_x h au point (x,t)(x - \ell, t)… L’avantage de prendre =Δx\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=0x = 0, il va falloir que j’évalue en xx - \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.

Pour le xx-\ell qui sort, c’est pas forcément un problème, tu as deux options simples (et des choses plus subtiles éventuellement, mais je n’ai pas assez réfléchi au problème physique) :

  • rendre ton espace circulaire, si on dit que ton espace est x[0,L]x\in[0,L], tu peux ensuite le prolonger par périodicité et dire que le point -\ell est le même que le point LL-\ell (et donc que le point 00 et le point LL sont identiques) ;
  • dire que tout est plat en dehors du domaine et donc avoir xh\partial_xh nul à l’extérieur.

Ensuite, pour évaluer xh\partial_x h en un point qui n’est pas forcément un des nœuds de calculs, ce n’est pas un problème non plus. Pour trouver l’expression en fonction de la valeur des deux points qui l’encadrent, écris la valeur de hh en ce point en fonction de celle des deux nœuds qui l’encadrent comme un DL au deuxième ordre et débarrasse toi des dérivées secondes par combinaison linéaire. Je peux te guider un peu plus si tu as du mal.

+0 -0

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 xh\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). :p Encore merci ! Je te tiens au courant de mes avancées.

Salut ! J’ai un peu commencé à mettre en place ce que tu m’as dit, mais je bloque au moment de l’approximation de xh\partial_x h en xx - \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 !

Salut,

Pour les premiers développements, tu les as écrit au premier ordre, pas au second. Si tu ajoutes les termes d’ordre deux, tu vas voir qu’ils partent en même temps que les termes d’ordre 0 et donc que ton erreur est en fait mieux que o(1)o(1), elle est o(δx)o(\delta x) (ou comme on le note plus souvent O(δx2)O(\delta x^2), d’où le fait que l’on dise que cette approximation est de deuxième ordre).

Pour la suite, tu te compliques pas mal la vie… Disons que le nœud le plus proche de x=xx_\ell=x-\ell est xkx_k avec

xkx=ε]δx2,δx2]x_k-x_\ell=\varepsilon\in\left]-\dfrac{\delta x}{2},\dfrac{\delta x}{2}\right]

On a

hk=h+ε(xh)+ε22(x2h)+o(δx2)h_k = h_\ell + \varepsilon(\partial_x h)_\ell + \frac{\varepsilon^2}{2}(\partial_x^2 h)_\ell+o(\delta x^2)

Puis sur les points de part et d’autre (je te laisse réfléchir à ce qu’il faut faire si xk1x_{k-1} est en dehors du domaine) :

hk+1=h+(δx+ε)(xh)+(δx+ε)22(x2h)+o(δx2)h_{k+1} = h_\ell + (\delta x+\varepsilon)(\partial_xh)_\ell + \frac{(\delta x+\varepsilon)^2}{2}(\partial_x^2 h)_\ell+o(\delta x^2)

hk1=h+(εδx)(xh)+(εδx)22(x2h)+o(δx2)h_{k-1} = h_\ell + (\varepsilon-\delta x)(\partial_xh)_\ell + \frac{(\varepsilon-\delta x)^2}{2}(\partial_x^2 h)_\ell+o(\delta x^2)

Il n’y a plus qu’à se débarrasser des termes hh_\ell et (x2h)(\partial_x^2 h)_\ell par combinaison linéaire, je te laisse essayer ça.

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 xjpx_{j - p} était le plus proche de xjx_j - \ell. Mais en effet, ça peut être xjp1x_{j - p - 1}. Je vais redéfinir mon ε\varepsilon.

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

Edit : Je n’arrive toujours pas à faire des combinaisons linéaires pour obtenir xh(x)\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. :)

+0 -0

Y’a un moyen simple de l’obtenir ?

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

Calcule hk+1hkh_{k+1}-h_k et hk1hkh_{k-1}-h_k pour te débarrasser des hh_\ell. Multiplie l’une des deux équations que tu obtiens par le facteur qui va bien pour te débarrasser des (x2h)(\partial_x^2h)_\ell par différence. Puis pour être sûr, vérifie que les cas particuliers ε=0\varepsilon=0 et ε=12δx\varepsilon=\frac12\delta x donnent un résultat cohérent avec le développement que tu as déjà pour (xh)k(\partial_xh)_k (i.e. celui où l’on évalue la dérivé première à mi-chemin de deux nœuds connus).

+1 -0

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

ça marche pour ε=0\varepsilon = 0 et ε=±12Δx\varepsilon = \pm\frac12\Delta x.

Teguad

Hmm, ça me parait bizarre que ton expression soit valide à la fois pour ε=12δx\varepsilon=\frac12\delta x et ε=12δx\varepsilon=-\frac12\delta x. Ces cas là, bien que correspondant à la même configuration, sont en fait différent en pratique (le nœud xkx_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é 12δx\frac12\delta x ou côté 12δx-\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.

+0 -0

Au temps pour moi, on obtient pas le bon résultat pour ε=12Δx\varepsilon = -\frac12\Delta x. Voici donc la relation que j’obtiens : xh(xi,t)=12δ(h(xk+1,t)2εδ[h(xk+1,t)2h(xk,t)+h(xk1,t)]h(xk1,t))+O(δ).\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é δ=Δx\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 !!!

Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte