Resolution EDP

Bonjour Dans le cadre d’un projet je dois résoudre cette équation aux dérivées partielles. J’ai vu plusieurs méthodes sur internet mais je ne les comprends ] et elle ne sont pas a mon programme (spe) Pouvez vous m’aider svp ? $\frac{\partial \rho(x, t)}{\partial t}+u_{\max } \frac{\partial \rho(x, t)}{\partial x}-\frac{u_{\max }}{\rho_{\max }} \frac{\partial \rho^{2}(x, t)}{\partial x}=\beta$

Beta est constante pour le temps et l’espace et les coefficients sont connus.

Salut,

J’ai un peu de mal à comprendre le tag sur ton sujet : tu dois mettre en place un programme Python pour calculer cette dérivée ?

Salut,

Elle est marrante cette équation, il représente quoi physiquement le terme en $\frac{\partial\rho^2}{\partial x}$ ? C’est pas exactement courant comme problème.

La résolution numérique d’EDP est pas exactement un truc qui s’improvise au doigt mouillé, donc ça pourrait être pas mal d’en savoir plus sur le contexte dans lequel ça t’est demandé.

Bonjour

Je m’intéresse a la capacité des ronds points. Voila le pef d’ou est tirée l’equation par principe de conservation. https://sites.math.washington.edu/~morrow/mcm/4329.pdf L’équation est obtenue en approximant la vitesse v par une fonction affine de rho. Je dois la résoudre sous python puisque c’est le seul language a ma disposition et mes recherches ne donnent aucune indication sur la résolution.

Euh… Si t’as jamais eu le moindre cours sur la résolution numérique des PDE, on part de loin tout de même. Comme tu as une équation de conservation sous la forme archi-classique $\partial_t\rho+\nabla\cdot\vec\phi=\beta$ et dans une géométrie simple, je serais tenté de dire que ton pari le plus simple est d’écrire ça avec une formulation en volume finis et un intégrateur explicite en temps.

L’idée générale des méthodes numériques est qu’on n’a pas de solution analytique à ce problème, et qu’on peut pas bêtement essayer toutes les fonctions possibles parce qu’il y en a une infinité. Donc on essaye de réduire l’espace des solutions via diverses stratégies. Une solution qui marche pas mal (celle des volumes finis) est de découper l’espace (ici l’axe $x$) en petits cubes (ici des segments) dans lesquels on va considérer que la valeur $\rho$ au centre des cubes est la moyenne de $\rho$ dans le cube, et on évalue le flux (donc ici $\rho u=\rho\left(1-\dfrac\rho{\rho_m}\right)u_m$) sur les faces (ici les points) entre les cubes (ici les segments).

Par exemple, si tu dis que ton domaine est $x\in[0,1)$ (avec $x=1$ étant le même point que $x=0$ comme on est sur un rond-point, et donc périodique), on va juste résoudre l’équation pour un ensemble de point $\rho_{i+1/2}$ situés en $x_{i+1/2}=(i+1/2)/N$ avec $i$ entier de $0$ à $N-1$, $N$ étant le nombre de segment de ton domaine (typiquement quelques dizaines, voire une centaine pour ton problème qui devrait être relativement lisse). Sur le segment $i$ (compris entre $x\in[x_i, x_{i+1}[=[i/N, (i+1)/N[$), de longueur $\delta x=1/N$, l’équation de conservation intégrée s’écrit

$\int_i\partial_t\rho + \int_i\nabla\cdot(\phi\hat x) = \int_i\beta$

On peut l’approximer comme expliqué ci-dessus en disant que la valeur moyenne de $\rho$ dans le segment $i$ est $\rho_{i+1/2}$ (et similaire pour $\beta$), ce qui donne

$\partial_t\rho_{i+1/2}\delta x + \phi_{i+1} - \phi_i = \beta_{i+1/2}\delta x$

Soit en réarrangeant un peu

$\partial_t\rho_{i+1/2} = -\dfrac{\phi_{i+1} - \phi_i}{\delta x} + \beta_{i+1/2}$

On connait $\beta_{i+1/2}$, et on peut calculer $\phi_i$ à partir des $\rho_{i\pm 1/2}$ :

$\phi_i = \rho_i\left(1-\dfrac{\rho_i}{\rho_m}\right)u_m \text{ avec }\rho_i\approx \dfrac{\rho_{i-1/2}+\rho_{i+1/2}}{2}$

(Cette approximation de $\rho_i$ par la moyenne des voisins fonctionne parce qu’on fait une erreur du même ordre que celle qu’on fait en supposant que $\rho_{i+1/2}$ est la moyenne de $\rho$ sur le segment, mais montrer ça dépasse le cadre de ce qu’on peut expliquer à l’arrache sur un message de forum).

On se retrouve alors avec un système d’ODE sur les $\rho_{i+1/2}$. Évidemment, comme on a envie de simplifier encore plus pour se retrouver avec des systèmes d’équations linéaires triviales à écrire et communiquer à un ordinateur, on va aussi discrétiser en temps. On approxime

$\partial_t\rho_{i+1/2}\approx \dfrac{\rho_{i+1/2}^{n+1}-\rho_{i+1/2}^{n}}{\delta t}$

où les exposants $n$ et $n+1$ ne sont pas des puissances mais dénotent des quantités évaluées à l’instant $n$ et à l’instant suivant. $\delta t=t^{n+1}-t^n$ étant le pas de temps.

Tout un pan des méthodes numériques s’intéresse alors à quel instant on va évaluer le membre de droite, $t^n$ ou $t^{n+1}$ (ou encore un mélange) ? Ici je suggère de pas se prendre la tête et d’intégrer de façon explicite, ce qui signifie que le membre de droite est évalué au moment $n$. Ça veut dire que si on a une solution à un instant $n$, c’est à dire un ensemble de valeurs $\rho_{i+1/2}^n$ connues, on peut calculer les $\rho_{i+1/2}^{n+1}$ au pas de temps suivant de façon directe :

$\rho_{i+1/2}^{n+1} =\rho_{i+1/2}^{n} + \delta t \left(-\dfrac{\phi_{i+1}^n - \phi_i^n}{\delta x} + \beta_{i+1/2}^n\right)$

Les conditions au bord périodiques permettent de calculer les flux en $i=0$ et $i=N$ comme on dit que $\rho_{N+1/2}=\rho_{1/2}$ et $\rho_{-1/2} = \rho_{N-1/2}$, et la condition initiale nous donne les $\rho_{i+1/2}$ à l’instant $n=0$. On peut alors calculer itérativement l’évolution de $\rho$ jusqu’au temps voulu.

D’accord merci beaucoup J’ai juste quelques questions. Mon domaine est plus compliqué qu’un simple segment si je faisais d’hypothèse d’un rond point carre cela me faciliterait la tache ? Ensuite dans votre dernière égalité je ne comprends pas comment accéder a la différence des phi, sachant que phi(i) dépend de rho( i) auquel je n’ai pas acces.

Mon domaine est plus compliqué qu’un simple segment si je faisais d’hypothèse d’un rond point carre cela me faciliterait la tache ?

Ça n’a pas d’importance parce que localement un cercle ressemble à une droite. Une autre façon de s’en convaincre est d’écrire la forme intégrée de l’équation de conservation sur une portion de tore (au lieu d’un segment), tous les termes géométriques dont la courbure vont se simplifier.

Pour $\rho_i$, c’est ce morceau qui est pertinent :

$\phi_i = \rho_i\left(1-\dfrac{\rho_i}{\rho_m}\right)u_m \text{ avec }\rho_i\approx \dfrac{\rho_{i-1/2}+\rho_{i+1/2}}{2}$

(Cette approximation de $\rho_i$ par la moyenne des voisins fonctionne parce qu’on fait une erreur du même ordre que celle qu’on fait en supposant que $\rho_{i+1/2}$ est la moyenne de $\rho$ sur le segment, mais montrer ça dépasse le cadre de ce qu’on peut expliquer à l’arrache sur un message de forum).

On le prend comme la moyenne des deux voisins $\rho_{i-1/2}$ et $\rho_{i+1/2}$.

Une question qui va peut être se poser est celle du choix du $\rho$ qui est effectivement advecté dans le terme $\phi = \rho u$. Là j’ai pris le parti de faire le choix le plus naïf de prendre $\phi_i = \rho_i u_i$, mais en fait tu peux avoir intérêt à prendre par exemple $\phi_i = \rho_{i-1/2} u_i$ (schéma d’advection dit upwind) ou un truc un peu plus subtil du genre TVD avec $\phi_i = u_i\left(\alpha \rho_{i-1/2}+(1-\alpha)\rho_{i+1/2}\right)$ et $\alpha$ qui dépend de $\rho$. Mais tout ça sont des questions que tu peux ignorer pour l’instant et partir sur l’option simple $\phi_i = \rho_i u_i$.

D’accord je vois merci. Mais je pense que cela est trop compliqué pour le niveau de mon projet, il est plutot axé physique. Il n’y aurait pas une bibliotheque python permettant de resoudre ceci simplement ?

Je pense que c’est justement la stratégie la plus simple, surtout que t’as pas besoin de rentrer dans les détails mathématique de la chose. Comme module Python qui pourrait faire ça, tu as FEniCSx mais :

• ça suppose tout de même de s’y connaitre un peu parce que c’est un solveur en éléments finis (ce qui est différent des volumes finis) qui demande donc de faire un choix sur les éléments, les fonctions de test (quoique j’imagine que c’est du Galerkin par défaut) et surtout de fournir ton problème en formulation faible au lieu de la formulation forte (une PDE locale quoi) à laquelle tu es habituée ;
• en plus de ça, comme c’est des éléments finis (encore !), ça suppose de faire un choix sur la stratégie de solveur direct que tu utilises ;
• si il s’avère que le solveur converge pas, c’est une boîte noire très difficile à contenter si tu es étranger aux méthodes numériques ;
• il sera pas forcément facile de valider les résultats que tu obtiens sans faire des calculs qui seront d’une complexité similaire voire plus grande que la stratégie que je propose.

Résoudre des PDE, c’est encore une fois pas un truc qui s’improvise comme ça, c’est vraiment un domaine des mathématiques qui est compliqué et pour lesquels il n’y a pas de solution toute prête qui marche à tous les coups sans devoir s’investir un minimum pour comprendre ce que fait le solveur (sauf les rares cas où on a des solutions analytiques…).

Ce qui m’intrigue un peu c’est qu’on t’a visiblement branché sur un projet pour lequel tu n’as pas les outils pour y répondre. Je pense honnêtement que si tu arrives à suivre ce que je raconte dans mon introduction rapide à la méthode des volumes finis, c’est la stratégie la plus rapide et simple pour que tu arrives à un résultat et que tu puisses explorer le problème physique.