Schéma numérique et ordre de convergence

Bonsoir à tous,

Je travaille actuellement quelques méthodes numériques pour approximer les équations aux dérivées partielles (que j'abrège EDP). En particulier, je m'intéresse à l'équation de transport (d'inconnue la fonction $u$) : $\partial_tu(x, t) + a(x)\partial_x u(x, t)$ pour $(x, t) \in [0, 2]\times[0, 1]$ et $u(x, 0) = e^{-30(x-1)^2}$. Je construis un schéma numérique, je démontre sa convergence ; jusque-là tout va bien. C'est ensuite que les choses se gâtent car j'ai un résultat un peu surprenant concernant le comportement de l'erreur de consistance en fonction de taille de mon maillage. Je note $\Delta_x$ le pas de maillage en espace et $\Delta_t$ celui en temps.

Pour approximer une solution $u$, je construis un maillage $(x_0, …, x_J)$ de l'espace et $(t_0, …, t_N)$ de l'espace. Je note $v_j^n$ la solution approchée de la solution $u$ au point $(x_j, t_n)$, et je construis le schéma suivant : $$\begin{aligned}v_j^{n+1} = v_j^n + \frac{\Delta_t}{\Delta_x} a(x_j)(v_j^n - v_{j-1}^n), &\quad j\in\{0, …, J\}, n\in\{0, N\}\\ v_j^0 = e^{-30(x_j-1)}, &\quad j\in\{1, …, J\}\end{aligned}$$

Enfin, pour pouvoir calculer effectivement les quantité $v_j^n$, je fais l'approximation que $v_0^n = 0$ car de toute façon $e^{-30}$ est vraiment très petit.

Bon, ça c'était pour les définitions et les notations. Je sais que ça fait beaucoup, mais ainsi va la vie. o/ J'ai démontré que mon schéma numérique est d'ordre 1 en temps et en espace. Plus précisément, $\frac{1}{\Delta_t}(u(x_j, t_n + \Delta_t) - u(x_j, t_n)) = O(\Delta_t) + O(\Delta_x)$.

Mais ensuite, j'ai quelque-chose d'un peu surprenant. Pour faire quelques tests, j'ai pris des fonctions $a$ simples (typiquement, $a$ constante égale à 1) et j'ai tracé le graphe en échelle log-log, de la quantité $||(v_j^n - u(x_j, t_n)_{1\le j\le J, 1\le n\le N}||_{\ell^\infty}$, en faisant varier $\Delta_x$ de façon géométrique (avec typiquement $\Delta_x = 3^{-k}$), pour $k = 2, …, 6$. $\Delta_t$ varie aussi, de façon automatique, de sorte à vérifier la relation $\frac{\Delta_t}{\Delta_x} = 1/3$ (ce qui est une condition de stabilité acceptable).

Ce qui me surprend beaucoup c'est que le graphe obtenu n'est pas linéaire ; à vrai dire il n'est même pas affine. Et c'est là que j'ai un gros doute. En effet, il me semble que normalement, l'erreur de consistance en échelle log-log est supposée être linéaire en le pas du maillage. En tout cas, je sais que c'est vrai lorsqu'il y a seulement une dimension d'espace et pas de temps. Que dire lorsqu'il y a une dimension d'espace et de temps ? A-t-on également linéarité de l'erreur en échelle log-log ? Est-ce le fait que mes courbes ne soient pas linéaires implique que mon implémentation du schéma est fausse ? Si non, avez vous des références bibliographiques ou des liens à me proposer, voire des idées de preuves, pour expliquer pourquoi l'erreur n'est pas linéaire dans ce cas ?

Merci beaucoup d'avance à tous. Je sais que la question peut paraître imbuvable, à cause des notations. Mais il y a quand même des maths intéressantes cachées derrière, c'est promis.

Salut,

En effet, il me semble que normalement, l'erreur de consistance en échelle log-log est supposée être linéaire en le pas du maillage.

Tant que tu es dans une tranche acceptable pour tes paramètres, oui. Si $\Delta x$ est trop petit, les erreurs de la machine vont faire que tu auras des oscillations bizarres. Au dessus d'une certaine taille, tu as saturation parce que l'erreur sera du même ordre de grandeur que la solution. On peut voir la tronche de ce que tu obtiens ?

Ce qui m'intrigue surtout, c'est ton critère de stabilité $\dfrac{\Delta t}{\Delta x}=\dfrac 13$… Déjà, ça devrait être $\dfrac{\Delta t}{\Delta x}<\dfrac 13$, non ? Et surtout, comment trouves tu cette valeur de $\dfrac 13$ ?

La condition sur $\frac{\Delta_t}{\Delta_x}$ est simplement la condition CFL (EN, j'ai rien trouvé de probant en français). Pour ceux qui auraient la flemme de se taper l'article, la condition CFL est une condition sur le rapport $\frac{\Delta_t}{\Delta_x}$ de stabilité en norme $\ell^\infty$.

Évidemment, cette quantité dépend de la fonction $a$, et avec le schéma que j'ai choisi, j'ai pu montrer que la condition CFL à considérer est $0 \le \max_{x\in[0, 2]} a(x) \frac{\Delta_t}{\Delta_x} \le 1$. Comme j'ai tracé les courbes avec la fonction constante $a = 1$ pour les tests, j'obtiens une CFL de type $0 \le \frac{\Delta_t}{\Delta_x} \le 1$. Du coup, j'ai choisi arbitrairement 1/3 ; j'aurais pu choisir n'importe quelle quantité entre 0 et 1, l'important étant que cela ne change pas l'ordre de l'erreur.

Enfin, voici l'allure d'une courbe d'erreur avec une CFL = 1/3 :

Pour bien te rassurer sur le fait que la CFL 1/3 fonctionne et qu'il n'est pas nécessaire d'exiger $\frac{\Delta_t}{\Delta_x} < 1/3$, j'ai aussi tracé l'erreur avec $\frac{\Delta_t}{\Delta_x} = 1/5$. Comme tu peux le voir, les deux courbes se ressemblent beaucoup :

On remarque bien que pour un pas de maillage assez grand, l'erreur tend à devenir linéaire, ce qui semble aller dans ton sens lorsque tu affirmes que le maillage doit être quand même assez précis. En outre, manifestement, l'ordre de l'erreur semble cohérent ; qu'en pensez-vous ?

P.-S. — Je n'ai pas insisté beaucoup sur les notions de conditions CFL car ce n'est pas le point qui me pose problème ici. Mais si vous êtes intéressé par la question et vous voulez en discuter, je suis prêt à créer un autre sujet pour cela.

Ah oui, j'ai oublié de préciser qui est $k$, désolé. J'ai défini $\Delta_x = 2^{-k-1}$ pour $k =1, …, 7$. Donc en gros, $k$ est le nombre de pas dans mon maillage en espace. J'ai choisi des puissances de 2, mais j'aurais pu mettre des puissances de n'importe quel entier. L'avantage de choisir des puissances et que ça me fait une suite géométrique en espace (et donc en temps aussi puisque $\Delta_t$ et $\Delta_x$ sont liés par la CFL.

Mais sauf erreur de ma part, k n'est donc pas un log, ici. Au contraire, c'est avoir une jolie courbe en passant au log-log que j'ai mis une suite géométrique sur $\Delta_x$. En même temps, lorsque je trace en échelle linéaire en k - log en erreur, je trouve une droite parfaite.

Du coup, il y a sans nul doute quelque chose qui m'échappe…

Hum, oui tout bêtement. Bien vu, et merci pour ton aide.

Pour des éventuels lecteurs à venir : les courbes que j'obtenais au début n'étaient pas des droites car j'avais représenté $k$ en échelle logarithmique, ce qui évidemment n'est pas l'approche la plus élégante.

Problème résolu, donc.

Ajout — En fait, quand tu simplifies $\log E$, tu n'as pas besoin de toutes ces étapes intermédiaires. On a directement $\log E = \log\Delta_x = \log(2^{-k-1}) = (-k-1)\log 2$.

M'en suis rendu compte juste après l'avoir écrit. Dans ma tête j'ai pensé à redéfinir l'exponentiation à partir du log en base 10 pour pouvoir sortir le $k$. Sur ce, je vais continuer mon bouquin "de l'art de se compliquer la vie pour les nuls".

Il n'empêche, tu m'as bien aidé. J'ai un peu honte de pas avoir détecté seul le problème, mais on va mettre ça sur le compte du « Je bosse sur les courbes depuis 2 jours non-stop ».

Juste une remarque, comme ton équation est une équation d'advection, je pense qu'il pourrait être judicieux de mettre en place un schéma upstream plutôt que se contenter de la dérivée spatiale à gauche. Tu gardes le même ordre de précision (avec la même condition de stabilité, nombre de Courant inférieur à 1), les calculs sont pas plus lourds mais tu gagnes en sens physique (et donc en précision). Tu peux gagner un ordre de précision avec beam-warming en prenant en compte la dérivée seconde à droite où à gauche.

Malheureusement, je n'ai pas le choix de la méthode, ce sur quoi je bosse s'insère dans un truc plus général et je choisis ce schéma pour des raisons purement scolaire. Je sais, c'est pas cool.

Ça me tue que quelqu’un comme toi, qui visiblement est assez à l’aise en math, puisse louper un truc aussi gros. Bon évidemment c’est plus facile vu de l’extérieur parce qu’on a pas la tête dans le guidon mais ça me fait tiquer quand même parce qu’il n’y a pas besoin de passer formellement par un calcul comme @dri1 l’a fait pour voir que $k$ est un log. Ça saute au yeux.

Bref, ce message ne sert à rien. Je fais juste mon petit vieux.

Ouais ouais, je sais. Mais bon, ça peut arriver à tout le monde. Je suis à peu près certain que 80 % des programmeurs expérimentés ont déjà passé un temps disons… anormalement long pour corriger une simple boucle itérative. C'est un peu pareil, là.