Table des matières

Étude de la stabilité et calcul automatique de CFL

Nous nous intéressons ici à résoudre un problème de transport du type :

$$ u_t + u_x = 0 $$

avec $x\in[0,1]$ et $t>0$, ainsi que la condition initiale $u_0$ donnée. Les méthodes que nous allons présentées s'étendent à des problèmes de transport non linéaire : $u_t + f_x(u) = 0$, mais la seule généralisation qui nous intéressera sera une équation de transport à vitesse constante $a$ : $u_t + au_x = 0$.

La discrétisation en espace se concentrera sur des méthodes WENO, ou des différences finies centrées d'ordre 2 ; la discrétisation en temps se concentrera sur des méthodes de type Runge-Kutta. Par la suite, pour effectuer un lien avec l'équation de Vlasov-Poisson, nous nous intéresserons à un problème généralisé du type :

$$ u_t = Lu + N(u) $$

où $L$ représentera la partie linéaire de notre problème et $N$ la partie non linéaire. Ce problème nous permettra d'étudier la résolution en temps avec des schémas exponentiels.

1. Discrétisations en espace

Pour étudier les schémas en espace nous allons utiliser l'analyse de von Neumann, pour cela introduisons quelques fonctions.

1.1 Méthodes d'ordre 1 et 2

Dans un premier temps nous allons présenter la méthode d'analyse sur le schéma upwind ou les différences finies centrées d'ordre 2.

Nous cherchons à discrétiser un problème du type :

$$ \partial_t u + a\partial_x u = 0 $$

Après une semi-discrétisation en espace, le schéma upwind s'écrit :

$$ \frac{\mathrm{d}u_j}{\mathrm{d}t} = -\frac{a^+}{\Delta x}\left(u_{j}-u_{j-1}\right)-\frac{a^-}{\Delta x}\left(u_{j+1}-u_{j}\right) $$

où $u_j(t)\approx u(t,x_j)$, $j=0,\dots,N$, $a^+=\max(a,0)$ et $a^-=\min(a,0)$.

Après une semi-discrétisation en espace, le schéma des différences finies centrées d'ordre 2 s'écrit :

$$ \frac{\mathrm{d}u_j}{\mathrm{d}t} = -\frac{a}{2\Delta x}\left(u_{j+1}-u_{j-1}\right) $$

où $u_j(t)\approx u(t,x_j)$, $j=0,\dots,N$.

La courbe pour le schéma upwind est orienté vers une partie réelle positive car il s'agit de l'approximation de $+\partial_xu$ or les méthodes Runge-Kutta que nous allons étudier par la suite s'applique à un problème écrit sous la forme $\dot{u} = L(t,u)$ donc il sera nécessaire d'approximer $-\partial_xu$.

On peut maintenant mettre en place les fonctions nécessaires pour la résultion numérique d'un problème de transport à vitesse constante (positive) :

$$ u_t + u_x = 0 $$

avec des conditions aux bords périodiques et une condition initiale $u(t=0,x) = u_0(x)$ donnée.

Maintenant, on effectue une mesure d'ordre, sur le même problème de transport avec comme condition initiale $u_0(x) = \cos(2\pi x)$.

Maintenant mettons en place ces outils pour le schéma WENO et ses variantes.

1.2 Présentation de la méthode WENO5

La méthode WENO est une famille de schémas aux différences finies (il existe une interprétation en terme de volumes finis) d'ordre potentiellement élevé, essentiellement non oscillante à l'aide de différentes approximations pondérées. De manière générale, la méthode WENO peut être utilisée pour discrétiser la dérivée spatiale dans un problème de transport du type :

$$ \partial_t u + \partial_x f(u) = 0 $$

Après une semi-discrétisation en espace, la méthode WENO consiste à approximer les flux numériques tels que :

$$ \frac{\mathrm{d}u_j}{\mathrm{d}t} = -\frac{1}{\Delta x}\left( \hat{f}_{j_\frac{1}{2}} - \hat{f}_{j-\frac{1}{2}} \right) $$

où $u_j(t) \approx u(t,x_j),\ j=0,\dots,N$. La méthode WENO, qui est une méthode up-wind, nécessite de splitter le flux entre sa partie positive et négative :

$$ f(u) = f^+(u) + f^-(u) $$

tel que $\frac{\mathrm{d}f^+}{\mathrm{d}u} \geq 0$ et $\frac{\mathrm{d}f^-}{\mathrm{d}u} \leq 0$. Nous noterons par la suite $f^\pm_j$ l\ approximation de $f^\pm(u(x_j))$ où $x_j$ est un point du maillage.

La méthode WENO consiste ensuite à estimer des poids en fonctions d'indicateurs de continuité (indicators of smoothness) $\beta_i^\pm$ (parfois noté $IS_i^\pm$), qui sont des combinaisons des carrés des approximations des dérivées successives en $x_j$ (sur des stencils potentiellement décentrés). Les coefficients $\beta_i^\pm$ permettent d'estimer des poids $w_i^\pm$, pondérant différentes estimations de la solution toujours en $x_j$ sur les différents stencils.

Pour le schéma standard WENO d'ordre 5 (aussi appelé WENO-JS pour [G.-S. Jiang et C.-W. Shu (1996)]) il y a 3 estimations sur 3 stencils différents :

$$ \begin{aligned} \beta_0^+ &= \frac{13}{12}(\underbrace{f^+_{j-2} - 2f^+_{j-1} + f^+_{j} }_{\Delta x^2(f''_j + \mathcal{O}(\Delta x))}))^2 + \frac{1}{4}( \underbrace{f^+_{j-2} - 4f^+_{j-1} + 3f^+_{j}}_{2\Delta x( f'_j + \mathcal{O}(\Delta x^2))} )^2 \\ \beta_1^+ &= \frac{13}{12}( \underbrace{f^+_{j-1} - 2f^+_{j} + f^+_{j+1}}_{\Delta x^2(f''_j + \mathcal{O}(\Delta x^2))} )^2 + \frac{1}{4}( \underbrace{f^+_{j-1} - f^+_{j+1}}_{2\Delta x f'_j + \mathcal{O}(\Delta x^2))})^2 \\ \beta_2^+ &= \frac{13}{12}( \underbrace{f^+_{j} - 2f^+_{j+1} + f^+_{j+2}}_{\Delta x^2(f''_j + \mathcal{O}(\Delta x))} )^2 + \frac{1}{4}(\underbrace{3f^+_{j} - 4f^+_{j+1} + f^+_{j+2}}_{-2\Delta x( f'_j + \mathcal{O}(\Delta x^2))})^2 \\ \end{aligned} $$

et de manière similaire :

$$ \begin{aligned} \beta_0^- &= \frac{13}{12}(f^-_{j+1} - 2f^-_{j+2} + f^-_{j+3})^2 + \frac{1}{4}(3f^-_{j+1} - 4f^-_{j+2} + f^-_{j+3})^2 \\ \beta_1^- &= \frac{13}{12}(f^-_{j} - 2f^-_{j+1} + f^-_{j+2})^2 + \frac{1}{4}( f^-_{j} - f^-_{j+2})^2 \\ \beta_2^- &= \frac{13}{12}(f^-_{j-1} - 2f^-_{j} + f^-_{j+1})^2 + \frac{1}{4}( f^-_{j-1} - 4f^-_{j} + 3f^-_{j+1})^2 \\ \end{aligned} $$

Ensuite nous calculons les poids non normalisés :

$$ \alpha_i^\pm = \frac{\gamma_i}{(\varepsilon + \beta_i^\pm)^2},\quad i=0,1,2 $$

où $\varepsilon$ est un paramètre numérique pour assurer la non nullité du dénominateur, il sera pris à $10^{-6}$ ; et vaec $\gamma_0=\frac{1}{10}$, $\gamma_1=\frac{6}{10}$ et $\gamma_2=\frac{3}{10}$. La normalisation des poids s'effectue comme suit :

$$ w_i^\pm = \frac{\alpha_i^\pm}{\sum_m \alpha_m^\pm},\quad i=0,1,2 $$

Nous pouvons ensuite calculer les flux numériques pour WENO5, donnés par :

$$ \begin{aligned} \hat{f}_{j+\frac{1}{2}}^+ =\ & w_0^+\left( \frac{2}{6}f^+_{j-2} - \frac{7}{6}f^+_{j-1} + \frac{11}{6}f^+_{j} \right) + w_1^+\left( -\frac{1}{6}f^+_{j-1} + \frac{5}{6}f^+_{j} + \frac{2}{6}f^+_{j+1} \right) \\ + & w_2^+\left( \frac{2}{6}f^+_{j} + \frac{5}{6}f^+_{j+1} - \frac{1}{6}f^+_{j+2} \right) \end{aligned} $$

et

$$ \begin{aligned} \hat{f}_{j+\frac{1}{2}}^- =\ & w_2^-\left( -\frac{1}{6}f^-_{j-1} + \frac{5}{6}f^-_{j} + \frac{2}{6}f^-_{j+1} \right) + w_1^-\left( \frac{2}{6}f^-_{j} + \frac{5}{6}f^-_{j+1} - \frac{1}{6}f^-_{j+2} \right) \\ + & w_0^-\left( \frac{11}{6}f^-_{j+1} - \frac{7}{6}f^-_{j+2} + \frac{2}{6}f^-_{j+3} \right) \end{aligned} $$

La méthode WENO5 prend la forme finale :

$$ (\partial_xf(u))(x_j) \approx \frac{1}{\Delta x}\left[ \left(\hat{f}_{j+\frac{1}{2}}^+ - \hat{f}_{j-\frac{1}{2}}^+ \right) + \left(\hat{f}_{j+\frac{1}{2}}^- - \hat{f}_{j-\frac{1}{2}}^- \right) \right] $$

Pour l'étude de la stabilité nous ne prendrons que le cas $f^+(u) = u$ et $f^-(u) =0$

1.2.1 Linéarisation du schéma WENO5

Le schéma WENO n'étant pas linéaire, son étude est en dehors du cadre classique de l'analyse de von Neumann. L'étude de la stabilité du schéma WENO avec différents schémas en temps, initié dans [R. Wang and R. J. Spiteri (2007)] et repris dans [M. Motamed and C. B. Macdonald and S. J. Ruuth (2010)] s'effectue à partir d'une linéarisation du schéma WENO, le calcul complet y est présenté, nous utiliserons ici que le résultat, issu d'un développement limité des poids non linéaire de la méthode WENO. Le développement limité nous donne une linéarisation des poids en :

$$ \alpha_i^{\pm} = \gamma_i + \mathcal{O}(\Delta x^2) $$

les poids $w_i^{\pm}$ s'écrivent alors :

$$ w_i^{\pm} = \gamma_i + \epsilon_i^{\pm} $$

avec $\epsilon_i = \mathcal{O}(\Delta x^2)$. Le coefficient d'amplification du schéma WENO linarisé s'écrit :

$$ \lambda^{\ell W}(\phi) = \tilde{z}(\phi) + M(\{\epsilon_i^{\pm}\}_i,\phi) $$

où $\tilde{z}$ correspond à la partie linéarie de WENO5 et $M$, dépendant de la famille des $\{\epsilon_i\}_i$, la partie non linéaire. Wang et Spiteri montrent qu'il est possible de majorer $M$ en remarquant que $M(\{\epsilon_i^{\pm}\}_i,\phi) = \mathcal{O}(\max_i(\epsilon_i))$. La fonction $\tilde{z}$ étant linéaire, il est possible d'effectuer son analyse de stabilité à l'aide de l'analyse de von Neumann.

Le flux de WENO5 non linéarisé fait intervenir les $\epsilon_i$ dont la contribution est difficilement quantifiable. On réécrit $\lambda^{\ell W} = \tilde{z} + \mathcal{O}(\Delta x^2)$, et on néglige tous les termes d'ordre $\Delta x^2$ et suppérieur.

On retrouve ainsi la courbe en haricot décrite dans [M. Motamed and C. B. Macdonald and S. J. Ruuth (2010)] ou [D. Ketcheson and C. B. Macdonald and S. J. Ruuth (2013)].

On peut maintenant s'intéresser aux non linéarités $M(\{\epsilon_i^\pm\}_i)$

En réalité la valeur des différentes composantes de $M$ nous importent peu. Ce qu'il faut remarquer est la forme général de chaque composante, prenons comme exemple la composante $\epsilon_0$ :

$$ \begin{aligned} \Re \left| \epsilon_0\left(\frac{11}{6} -\frac{7}{6}e^{-i\phi} + \frac{1}{3}e^{-2i\phi}\right)\right| &\leq c_0^{\Re}|\epsilon_0| \\ \Im \left| \epsilon_0\left(\frac{11}{6} -\frac{7}{6}e^{-i\phi} + \frac{1}{3}e^{-2i\phi}\right)\right| &\leq c_0^{\Im}|\epsilon_0| \end{aligned} $$

avec $c_0^{\Re|\Im}$ des constantes positives. Dans [Wang R. and Spiteri R. J. (2007)], une estimation de $c_0^{\Re}=\frac{10}{3}$ est prise comme exemple. On majore ainsi toutes les composantes en $\epsilon_i$ de $M$, on se retrouve alors avec :

$$ \begin{aligned} \left|\Re\left(M(\{\epsilon_i\}_{i=0,\dots,5},\phi)\right)\right| &\leq \Gamma^{\Re}\max_{0\leq m \leq 5}|\epsilon_m| \\ \left|\Im\left(M(\{\epsilon_i\}_{i=0,\dots,5},\phi)\right)\right| &\leq \Gamma^{\Im}\max_{0\leq m \leq 5}|\epsilon_m| \end{aligned} $$

avec $\Gamma^{\Re|\Im}$ des constantes postives ne dépendant que de la taille du stencil considéré. Cela permet de justifier $M = \mathcal{O}(\Delta x^2)$.

On peut vérifier quelques propriétés de ce schéma WENO linéarisé. Ce schéma correspond à l'approximation des poids $w_i^\pm = \gamma_i$, c'est-à-dire que les indicators of smoothness valent tous 1, ce qui se justifie avec une fonction suffisamment régulière.

On implémente dans un premier temps nos schémas en espace :

Puis le schéma en temps, on se contentera d'un schéma de type Runge-Kutta d'ordre 3 à 3 étages : RK(3,3) (avec la fonction rk33). Ce schéma est d'ordre inférieur en temps, l'erreur sattuera donc rapidement en $\Delta t^3$.

On teste nos schémas avec un problème de transport :

$$ \begin{cases} u_t + u_x = 0 \\ u(t=0,x) = u_0(x) = \cos(x) \end{cases} $$

On peut maintenant calculer l'ordre de la méthode sur ce problème de transport, en calculant la solution pour différents $\Delta x$ (on choisit différentes tailles de vecteur en entré), jusqu'au temps $T_f=1$ avec un pas de temps $\Delta t = \frac{1}{400}$ qui correspond environ à $\Delta t \leq \frac{1}{2}\Delta x_{\text{min}}$, pour minimiser l'erreur en temps.

On remarque bien une pente d'ordre 5 avec le schéma WENO linéarisé.

La linéarisation du schéma WENO suppose que la à dériver est suffisamment régulière. On peut donc chercher à étudier son comportement face à une discontinuité, en transportant à vitesse 1 une fonction test.

Le cas test de Shu est donné par la condtion initiale suivante :

$$ u_0(x) = \begin{cases} \frac{1}{6}\left( G(x,\beta,z-\delta) + G(x,\beta,z+\delta) + 4G(x,\beta,z) \right) & -0.8 \leq x \leq -0.6 \\ 1 & -0.4 \leq x \leq -0.2 \\ 1 - \left| 10(x-0.1) \right| & 0 \leq x \leq 0.2 \\ \frac{1}{6}\left( F(x,\alpha,a-\delta) + F(x,\alpha,a+\delta) + 4F(x,\alpha,a) \right) & 0.4 \leq x \leq 0.6 \\ 0 & \text{sinon} \end{cases} $$

avec : $$G(x,\beta,z) = e^{-\beta(x-z)^2}$$ $$F(x,\alpha,a) = \sqrt{\max(1-\alpha^2(x-a)^2,0)}$$ et les constantes sont données par $a = 0.5$, $z=-0.7$, $\delta = 0.005$, $\alpha=10$ et $\beta = \frac{\ln(2)}{36\delta^2}$

Ce cas test permet de tester la réponse d'un schéma face à différents profiles de fonctions. Il y a ainsi une fonction continue, lisse, mais relativement fine, une discontinuité, une discontinuité de gradient constant par morceaux, et une discontinuité de gradient avec un fort gradient.

La simulation s'éffectue jusqu'au temps $T_f=1$ sur le domaine périodique $[-1,1]$, avec $N_x = 100$ le nombre de points de discrétisation, et la condition CFL est donnée par $\Delta t = 0.1\Delta x$, toujours avec le schéma Runge-Kutta d'ordre 3 RK(3,3).

On remarque de manière classique que le schéma upwind diffuse énormément, et que le schéma CD2 oscille face à une discontinuité et que celles-ci se propagent sur tout le domaine. Le schéma WENO5 linéarisé oscille lorsque la variation de gradient est importante ou que le gradient présente une dicontinuité, et diffuse peu. Le schéma WENO5 complet devient plus visqueux à l'aide des poids qui permettent d'effectuer une meilleure interpolation qu'une interpolation polynomiale de degré élevé induisant localement des oscillations, mais il diffuse légèrement plus que la version linéarisé.

On peut étudier la réponse du schéma WENO linéarisé à un cosinus. La méthode servant à approximer une dérivée, on peut vérifier ce comportement, connaissant la solution exacte.

Le résultat semble satisfaisant, regardons maintenant la différence entre les deux. L'erreur commise par WENO est de $10^{-8}$, avec un $\Delta x = \frac{2\pi}{100}$ (soit $\Delta x^{6}$ ce qui est attendu avec un schéma d'ordre 5).

On trouve bien une différence pas très loin de l'ordre de $\Delta x^5$ ce qui nous permet, selon la définition de l'ordre, de considérer le schéma comme d'ordre 5.

L'étude se faisant automatiquement avec sympy, il est simple d'effectuer cette étude de coefficients de Fourier avec le schéma WENO standard non linéarisé. L'étude de l'analyse de von Neumann aussi mais il est impossible de justifier mathématiquement les résultats obtenus.

On peut maintenant tracer la différence entre les coefficients de Fourier des 2 résultats, entre WENO et WENO linéarisé.

1.3 Méthodes WENO modifiées

Il existe différentes méthodes de calcul des poids du schéma WENO d'ordre 5, qui ont été listées et comparées ici par Manuel A. Diaz. Nous étudirons ici leur différents coefficients d'amplification.

La différence réside uniquement dans l'estimation des poids, les indicateurs de continuité $\beta_i$ ainsi que les coefficients $\gamma_i$ restent identique.

Méthode Calcul des poids $w_i$ Parmètres
WENO-JS $$\begin{aligned}\alpha_i &\gets \frac{\gamma_i}{(\epsilon + \beta_i)^2} \\ w_i &\gets \frac{\alpha_i}{\sum \alpha_k} \end{aligned}$$ $\epsilon = 10^{-6}$
WENO-M $$\begin{aligned}\alpha_i &\gets \frac{\gamma_i}{(\epsilon + \beta_i)^2} \\ w_i &\gets \frac{\alpha_i}{\sum \alpha_k} \\ g_i &\gets w_i\left(\frac{\gamma_i + \gamma_i^2 - 3w_i\gamma_i + w_i^2}{\gamma_i^2+w_i(1-2\gamma_i)}\right) \\ w_i &\gets \frac{g_i}{\sum g_k} \end{aligned}$$ $\epsilon = 10^{-6}$
WENO-Z $$\begin{aligned}\alpha_i &\gets \gamma_i\left(1+\frac{\tau_5}{\epsilon + \beta_i}\right) \\ w_i &\gets \frac{\alpha_i}{\sum \alpha_k} \end{aligned}$$ $\epsilon = 10^{-40}$, $\tau_5 = \beta_0-\beta_2 $

La méthode WENO-JS est la méthode WENO classique de Jiang et Shu et sert juste pour la comparaison avec WENO-M (WENO avec une fonction de mappage) et WENO-Z. L'objectif de ces méthodes est de minimiser la perte d'ordre à l'approche d'une discontinuité, donc d'approter moins de viscosité numérique.

1.3.1 WENO-M

La méthode WENO-M fonctionnne à l'aide d'une fonction dite de mappage : $g_i$. Cette fonction est définie par :

$$ g_i : w \mapsto w\left( \frac{\gamma_i + \gamma_i^2 - 3w\gamma_i + w^2}{\gamma_i^2 + w(1-2\gamma_i)} \right) $$

Les fonctions $g_i$ n'ont besoin d'être définie que sur l'intervalle $[0,1]$ puisque les poids d'une méthode WENO vérifient $\sum_k w_k=1$ et $w_k\leq0\,\forall k$. Les points fixes de la fonction de mappage $g_i$ sont les 2 valeurs triviales 0 et 1, ainsi que $\gamma_i$ (la valeur linéarisé du poids). La fonction $g_i$ sera évaluée en $w_i = \gamma_i + \mathcal{O}(\Delta x^2)$, ce qui nous intéresse est donc son évaluation en des points proches de $\gamma_i$. On remarque que la dérivée de $g_i$ est nulle (donc un plateau) en $\gamma_i$ et que $g_i(\gamma_i) = \gamma_i$, cela permet de forcer une certaine linéarisation des poids, sans pour autant faire apparaître des oscillations.

Pour une étude plus approfondie de cette méthode avec l'analyse de von Neumann, on peut remarquer que la linéarisation des poids donne toujours $w_i = \gamma_i + \mathcal{O}(\Delta x^2)$. L'analyse qui en suit sera identique. On regardera donc son comportement sur une équation de transport avec comme condition initiale le test de Shu.

On vérifie l'ordre de la méthode, qui est bien d'ordre 5.

1.3.2 WENO-Z

La méthode WENO-Z consiste à minimiser la baisse d'ordre à l'approche d'une discontinuité. Dans le cadre de la simulation de plasmas, nous savons que nos fonctions de distribution seront continue, l'ordre élevé est nécessaire pour capturer correctement les forts gradients.

Pour l'analyse de stabilité, la linéarisation des poids donne aussi $w_i = \gamma_i + \mathcal{O}(\Delta x^2)$, l'analyse de von Neumann du schéma linéarisé reste donc la même.

On représente tous les schémas WENO d'ordre 5, ainsi que la version linéarisé pour les comparer. On effectue ce test sur le cas test de Shu, dont la condition initiale est donnée par :

$$ u_0(x) = \begin{cases} \frac{1}{6}\left( G(x,\beta,z-\delta) + G(x,\beta,z+\delta) + 4G(x,\beta,z) \right) & -0.8 \leq x \leq -0.6 \\ 1 & -0.4 \leq x \leq -0.2 \\ 1 - \left| 10(x-0.1) \right| & 0 \leq x \leq 0.2 \\ \frac{1}{6}\left( F(x,\alpha,a-\delta) + F(x,\alpha,a+\delta) + 4F(x,\alpha,a) \right) & 0.4 \leq x \leq 0.6 \\ 0 & \text{sinon} \end{cases} $$

avec : $$G(x,\beta,z) = e^{-\beta(x-z)^2}$$ $$F(x,\alpha,a) = \sqrt{\max(1-\alpha^2(x-a)^2,0)}$$ et les constantes sont données par $a = 0.5$, $z=-0.7$, $\delta = 0.005$, $\alpha=10$ et $\beta = \frac{\ln(2)}{36\delta^2}$

La simulation s'effectue jusqu'au temps $T_f=1$.

On remarque que la version linéarisé réagit très bien à la discontinuité du gradient mais oscille au niveau de la discontinuité. Le schéma WENO-Z est bien celui minimisant la perte d'ordre, donc la viscosité, sans osciller au niveau de la dicontinuité.

Avant de faire des tests en temps long, nous allons vérifier l'ordre de la méthode WENO-Z.

Maintenant nous allons effectuer un test en temps long, $T_f=100$ avec le cas test de Shu, donné par la condtion initiale suivante :

$$ u_0(x) = \begin{cases} \frac{1}{6}\left( G(x,\beta,z-\delta) + G(x,\beta,z+\delta) + 4G(x,\beta,z) \right) & -0.8 \leq x \leq -0.6 \\ 1 & -0.4 \leq x \leq -0.2 \\ 1 - \left| 10(x-0.1) \right| & 0 \leq x \leq 0.2 \\ \frac{1}{6}\left( F(x,\alpha,a-\delta) + F(x,\alpha,a+\delta) + 4F(x,\alpha,a) \right) & 0.4 \leq x \leq 0.6 \\ 0 & \text{sinon} \end{cases} $$

avec : $$G(x,\beta,z) = e^{-\beta(x-z)^2}$$ $$F(x,\alpha,a) = \sqrt{\max(1-\alpha^2(x-a)^2,0)}$$ et les constantes sont données par $a = 0.5$, $z=-0.7$, $\delta = 0.005$, $\alpha=10$ et $\beta = \frac{\ln(2)}{36\delta^2}$

Le deuxième test s'effectue sur une fonction chapeau, en temps court celle-ci ne présente pas de réel intérêt, en temps plus long on y voit apparaître les défauts des différents schémas dans des conditions plus favorable d'un créneau.

$$ u(t=0,x) = u_0(x) = \begin{cases} 4x-1 & \text{si } x\in[\frac{1}{4},\frac{1}{2}] \\ -4x+3 & \text{si } x\in[\frac{1}{2},\frac{3}{4}] \\ 0 & \text{sinon} \end{cases} $$

1.4 B-WENO

[Banks J. W. et al (2019)] est proposé une modification du schéma WENO plus intéressante pour l'étude de l'équation de Vlasov. Nous étudierons ici la présentation effectuée dans la seconde publication. Dans la continuation des schémas d'ordre élevé en espace, nous nous intéresserons au schéma BWENO d'ordre 6 présenté.

L'idée de ce schéma reste similaire à celle du WENO dans l'esprit du calcul de poids non linéaire s'adaptant aux irrégularités de la fonction. Le nombre de poids est en revanche limité à 2, et ceux-ci ne sont plus des polynôme en $\left(f_{i+k}^4\right)_{k\in [\![ -3,3 ]\!]}$ mais seulement de degré 2 : $\left(f_{i+k}^2\right)_{k\in [\![ -3,3 ]\!]}$. Cette diminution du nombre de multiplication peut rendre ce schéma compétitif au niveau du temps de calcul. Les solutions attendues dans la simulation de plasmas avec l'équation de Vlasov-Poisson sont régulières, l'ordre élevé des schémas sert à capturer les forts gradients de la solution qui peuvent apparaître lors de la filamentation.

Le flux $\hat{u}^{(B)}_{j-\frac{1}{2}}$ du schéma BWENO d'ordre $p$ s'obtient par l'addition pondéré de 2 interpolations $L$ et $R$, d'ordre $p-1$. Nous ne nous intéresserons ici qu'à la méthode d'ordre 6, nous nous abstiendrons donc d'indiquer l'ordre.

$$ \hat{u}^{(B)}_{j-\frac{1}{2}} = w^{(L)}_{j-\frac{1}{2}}u^{(L)}_{j-\frac{1}{2}} + w^{(R)}_{j-\frac{1}{2}}u^{(R)}_{j-\frac{1}{2}} $$

L'écriture du schéma commence par l'estimation des indicateurs de continuité :

$$ \begin{aligned} \beta(u_i) = (\Delta^1_4u_i)\left( (\Delta^1_4u_i) + (\Delta^2_4u_i) + \frac{1}{3}(\Delta^3_2u_i) + \frac{1}{12}(\Delta^4_2u_i) \right) \\ + (\Delta^2_4u_i)\left( \frac{4}{3}(\Delta^2_4u_i) + \frac{5}{4}(\Delta^3_2u_i) + \frac{2}{5}(\Delta^4_2u_i) \right) \\ + (\Delta^3_2u_i)\left( \frac{83}{60}(\Delta^3_2u_i) + \frac{23}{18}(\Delta^4_2u_i) \right) \\ + \frac{437}{315}(\Delta^4_2u_i)^2 \end{aligned} $$

où les $(\Delta^d_pu_i)$ sont des approximations d'ordre $p$ de la dérivée d'ordre $d$ :

$$ \begin{aligned} (\Delta^1_4u_i) &= \frac{1}{12}\left( -u_{j+2} + 8u_{j+1} - 8u_{j-1} + u_{j-2} \right) \\ (\Delta^2_4u_i) &= \frac{1}{12}\left( -u_{j+2} + 16u_{j+1} - 30u_{j} + 16u_{j-1} - u_{j-2} \right) \\ (\Delta^3_2u_i) &= \frac{1}{2}\left( u_{j+2} - 2u_{j+1} + 2u_{j-1} - u_{j-2} \right) \\ (\Delta^4_2u_i) &= \left( u_{j+2} - 4u_{j+1} + 6u_{j} - 4u_{j-1} + u_{j-2} \right) \\ \end{aligned} $$

Nous prendrons $\beta^{(L)}_{j-\frac{1}{2}} = \beta(u_{j-1})$ et $\beta^{(R)}_{j-\frac{1}{2}} = \beta(u_j)$ comme indicateurs de continuité respectivement à gauche et à droite.

Les 2 poids non normalisés sont :

$$ a^{(\Xi)}_{j-\frac{1}{2}} = \frac{d}{\epsilon + \beta^{(\Xi)}_{j-\frac{1}{2}}} $$

où $\Xi = L,R$, $d=\frac{1}{2}$ et $\epsilon$ un paramètre pour éviter la nullité du dénominateur. Ce qui nous donne les poids :

$$ \tilde{w}^{(\Xi)}_{j-\frac{1}{2}} = \frac{a^{(\Xi)}_{j-\frac{1}{2}}}{a^{(L)}_{j-\frac{1}{2}}+a^{(R)}_{j-\frac{1}{2}}} $$

Les poids sont ensuites classé selon la direction, ici pour une vitesse positive (dans l'étude du schéma nous prendrons $v=1$) :

$$ \begin{cases} w^{(L)} = \max\left(\tilde{w}^{(L)}_{j-\frac{1}{2}},\tilde{w}^{(R)}_{j-\frac{1}{2}}\right) \\ w^{(R)} = \min\left(\tilde{w}^{(L)}_{j-\frac{1}{2}},\tilde{w}^{(R)}_{j-\frac{1}{2}}\right) \end{cases} $$

Les flux à gauche et à droite sont maintenant calculés comme suit :

$$ \begin{cases} \hat{u}^{(L)}_{j-\frac{1}{2}} = \frac{1}{60}\left( -3u_{j+1}+27u_{j}+47u_{j-1} -13u_{j-2} + 2u_{j-3} \right) \\ \hat{u}^{(R)}_{j-\frac{1}{2}} = \frac{1}{60}\left( 2u_{j+2}-13u_{j+1}+47u_{j}+27u_{j-1}-3u_{j-2} \right) \\ \end{cases} $$

On peut enfin calculer une approximation de $\partial_x u$ à l'aide du flux de BWENO par :

$$ \partial_x u \approx \frac{1}{\Delta x}\left( \hat{u}^{(B)}_{j+\frac{1}{2}} - \hat{u}^{(B)}_{j-\frac{1}{2}} \right) $$

Pour l'étude du schéma BWENO, l'article [Banks J. W. et al (2019)] propose une linéarisation du schéma en prenant les poids $w^L$, $w^R$ égaux à différents couples de nombres positifs vérifiant que $w^L+w^R=1$, ainsi est tracé le coefficient d'amplification pour $w^L=1,\frac{3}{4},\frac{1}{2}$. En effet différents éléments empêche facilement d'écrire les poids sous la forme $w_i = \gamma_i + \mathcal{O}(\Delta x^2)$ comme les précédentes modifications du schéma WENO.

En linéarisant les poids ($w^L=1,\frac{3}{4},\frac{1}{2}$), on remarque un comportement similaire à celui du schéma WENO linéarisé (qui coïncide parfaitement pour $w^L=1$, $w^R=0$), avec un comportement dégénéré dans le cas centré $w^L = w^R = \frac{1}{2}$ (transport pur sans diffusion). La partie imaginaire, même dans le cas non linéarisé, la méthode BWENO reproduit le même comportement que la méthode WENO linéarisée.

On obtient un schéma qui diffuse moins que le schéma WENO-Z, il trouve ainsi son intérêt pour approximer au mieux des pics. Cela se fait au détriment du caractère oscillant.

On peut vérifier l'ordre du schéma. Pour cela on calcule l'erreur sur une advection à vitesse 1 :

$$ \begin{cases} u_t + u_x = 0 \\ u^0(x) = u(t=0,x) = \cos(2\pi x) \end{cases} $$

On peut maintenant calculer l'ordre de la méthode sur ce problème de transport, en calculant la solution pour différents $\Delta x$ (on choisit différentes tailles de vecteur en entré), jusqu'au temps $T_f=1$ avec un pas de temps $\Delta t = \frac{1}{400}$ qui correspond environ à $\Delta t \leq \frac{1}{2}\Delta x_{\text{min}}$, pour minimiser l'erreur en temps.

[Banks J. W. et al (2019)] annonce un ordre 6, mais ce n'est pas flagrant numériquement. Est-ce juste une définition différente de l'ordre ?

1.5 WENO3

Après avoir présenté les méthodes WENO dites classiques que sont la méthode WENO d'ordre 5 et ses dérivées, attardons nous un peu sur la méthode WENO3, étrangement moins étudiée. La méthode est beaucoup plus succinte mais se présente de manière similaire avec un calcul d'indicateurs de continuité, puis de poids que l'on normalise, et enfin le calcul du flux.

Les indicateurs de continuités $\beta_0$ et $\beta_1$ se calculent comme suit :

$$ \begin{aligned} \beta_0^+ &= \left(-f_{i-1}^+ + f_{i}^+\right)^2\\ \beta_1^+ &= \left(-f_{i}^+ + f_{i+1}^+\right)^2 \end{aligned} $$$$ \begin{aligned} \beta_0^- &= \left(-f_{i+2}^- + f_{i+1}^-\right)^2\\ \beta_1^- &= \left(-f_{i+1}^- + f_{i}^-\right)^2 \end{aligned} $$

Ensuite les poids :

$$ \alpha_i^\pm = \frac{\gamma_i}{(\epsilon+\beta_i^\pm)^2} $$

avec $\gamma_0 = \frac{1}{3}$ et $\gamma_1=\frac{2}{3}$ ; poids que l'on normalisent :

$$ w_i^\pm = \frac{w_i}{w_0+w_1} $$

Enfin le flux, définit par :

$$ f_{i+\frac{1}{2}}^+ = w_0^+\left(-\frac{1}{2}f_{i-1}^+ + \frac{3}{2}f_i^+\right) + w_1^+\left(\frac{1}{2}f_i^++\frac{1}{2}f_{i+1}^+\right) $$$$ f_{i+\frac{1}{2}}^- = w_0^-\left(-\frac{1}{2}f_{i+2}^- + \frac{3}{2}f_{i+1}^-\right) + w_1^-\left(\frac{1}{2}f_{i+1}^-+\frac{1}{2}f_{i}^-\right) $$

La méthode WENO3 prend la forme finale :

$$ (\partial_xf(u))(x_j) \approx \frac{1}{\Delta x}\left[ \left(f_{j+\frac{1}{2}}^+ - f_{j-\frac{1}{2}}^+ \right) + \left(f_{j+\frac{1}{2}}^- - f_{j-\frac{1}{2}}^- \right) \right] $$

Pour l'étude de la stabilité nous ne prendrons que le cas $f^+(u) = u$ et $f^-(u) =0$

Comme indiqué dans la littérature, la méthode WENO3 diffuse plus que son équivalent d'ordre 5, cela se remarque par une plus grande partie imaginaire.

Voilà qui conclu l'étude des schémas en espace.

2. Discrétisation en temps

2.1 Méthodes Runge-Kutta explicites

Les schémas Runge-Kutta sont souvent utilisés dans le domaine des EDP pour résoudre la discrétisation en temps, conforme à cette tradition c'est ce que nous allons effectuer. Nous cherchons quel est l'ordre minimale $n$ ainsi que le nombre d'étages $s$ pour stabiliser le couple RK($s$,$n$)-WENO5, ou celui qui permettra d'obtenir la plus grande CFL avec le minimum de coût numérique.

On s'intéresse ici à la fonction de stabilité d'un schéma RK. On sait que le polynôme caractéristique d'un schéma RK($n$,$n$) est une troncature de la série entière de l'exponentielle (auquel peut s'ajouter des termes de degré plus élevé si le nombre d'étage $s$ est supérieur à l'ordre de la méthode $n$), donc de la forme :

$$ p_{(n,n)}(z) = \sum_{k=0}^n \frac{z^k}{k!} $$

Il est ensuite nécessaire de tracer la courbe d'équation $|p_{(n,n)}(z)| = 1$, pour cela Miguel m'a proposé de résoudre dans un premier temps l'équation :

$$ p_{(n,n)}(z) = e^{i\theta} $$

à $\theta$ fixé, on a ainsi $n$ solutions $(b_0(\theta),\dots,b_n(\theta))$. Puis ensuite faire varier $\theta \in [0,2\pi[$. On obtient ainsi $n$ courbes paramétriques que l'on peut tracer sans trop de problèmes.

Ces fonctions de stabilités ne sont que celles d'une méthode RK($n$,$n$), avec $n<5$ pour pouvoir résoudre de façon exacte les racines du polynôme. Il est possible de tracer directement un contour avec matplotlib à l'aide de la fonction plt.contour. On peut ainsi obtenir, en un temps relativement court, un ensemble de points. En plus de tracer cet ensemble de points, il est également possible de récupérer la liste des points tracés, à l'aide de l'attribut allsegs de la classe matplotlib.contour.QuadContourSet.

Attention : matplotlib permet de récupérer la liste des points tracés, si le domaine d'étude est trop resteint, seuls les points se trouvant dans ce domaine seront listés. Comme son nom l'indique l'attribut allsegs renvoie des listes de points qui forment les segements tracés, par conséquent il est possible de récupérer chaque sous ensemble connexe (ce qui permet d'exclure l'étude des non-connexités des méthodes Runge-Kutta avec plus d'étages.

L'astuce trouvé pour obtenir le domaine principal et qui sera utilisé par la suite, est de liste les sous domaine allsegs et de les trier par nombre de points. Le domaine principal qui permettra de stabiliser un schéma type WENO, ou méthode centré par exemple, est (quasi-systématiquement) celui ayant le plus de points.

Cette méthode n'a rien de mathématique, juste une constatation empirique.

Joackim m'a fait remarquer que les bras contiennent une infinité de bulles de plus en plus petites (on voit ceci pour la petite bulle orange à l'extrémité du bras qui apparait seulement lorsque on raffine le domaine). On peut même conjecturer que la forme des bras sont des paraboles. On voit la naissance de la deuxième paire de bras avec les grosses bulles oranges.

L'étude de ces non-connexités des domaines de stabilités ne semblent pas présenter d'intérêt pour ce que l'on fait. De plus ces domaines ne correspondent pas à la stabilité d'une méthode connue. En effet il n'y a pas de méthode Runge-Kutta d'ordre $n$ à $n$ étages connues pour $n\geq 5$.

L'étude de ces méthodes RK($n$,$n$) avec $n>4$ n'est pas pertinante puisqu'il n'existe pas de méthode Runge-Kutta d'ordre 5 à 5 étages, le problème est similaire pour les méthodes d'ordre supérieur. À partir de maintenant nous allons étudier ces méthodes à partir de leur écriture sous forme de schéma ou de tableau de Butcher. Par la suite, le tableau de Butcher sera plus pratique pour tirer plus d'information comme le schéma de Lawson induit.

2.1.1 Obtention de la fonction de stabilité depuis un schéma numérique

Les calculs effectués pour l'instant ne s'appliquent que pour un RK($n$,$n$) à $n$ étages. Considérons maintenant un RK($s$,$n$) à $s$ étages ($s\geq n$ car la méthode est explicite). Il est nécessaire de calculer sa fonction de stabilité (qui est un polynôme dans le cadre des méthodes explicites).

L'obtention de la fonction de stabilité s'effectue avec une fonction $L$ linéaire, ainsi pour un problème du type :

$$ \frac{\mathrm{d}u}{\mathrm{d}t} = L(t,u) $$

on effectue la substitution :

Nous étudions ainsi un schéma RK(4,3) :

$$ \begin{aligned} u^{(1)} &= u^n + \frac{1}{2}\Delta t L(t^n,u^n) \\ u^{(2)} &= u^{(1)} + \frac{1}{2}\Delta t L(t^n+\frac{1}{2}\Delta t,u^{(1)}) \\ u^{(3)} &= \frac{2}{3}u^n + \frac{1}{3}u^{(2)} + \frac{1}{6}\Delta t L(t^n+\Delta t,u^{(2)}) \\ u^{n+1} &= u^{(3)} + \frac{1}{2}\Delta t L(t^n+\frac{1}{2}\Delta t,u^{(3)}) \end{aligned} $$

La méthode sera dite de Shu-Osher si elle n'effectue qu'un seul appel à la fonction $L$ par étage, on verra plus tard que pour un tableau de Butcher donné, il est facilement possible de shuosheriser la méthode.

2.1.2 Obtention de la fonction de stabilité depuis un tableau de Butcher

Il est souvent plus simple de ne donner qu'un tableau de Butcher d'une méthode, et cela permet en réalité d'extraire plus d'informations. Dans un premier temps, l'implémentation, se limitera aux méthodes eRK (explicit Runge-Kutta method), par la suite cette implémentation sera généralisé au cadre des méthodes DIRK (Diagonal Implicit Runge-Kutta method), le cadre général étant plus embêtant pour l'obtention algorithmique de la fonction de stabilité.

La fonction poly_butcher fonctionne pour un tableau de Butcher écrit comme suit :

$$ \begin{array}{c|c} \begin{matrix} a_{11} & \cdots & a_{1s} \\ \vdots & \ddots & \vdots \\ a_{s1} & \cdots & a_{ss} \end{matrix} & \begin{matrix} b_1 \\ \vdots \\ b_s \end{matrix} \\ \hline \begin{matrix}c_1 & \cdots & c_s \end{matrix} & \\ \end{array} $$

Le schéma qui en résulte est de la forme :

$$ \begin{aligned} u^{(i)} &= u^n + \Delta t \sum_j a_{ij} L(t^n+b_j\Delta t , u^{(j)}) \\ u^{n+1} &= u^n + \Delta t \sum_i c_i L(t^n+b_i\Delta t , u^{(i)}) \end{aligned} $$

Ce schéma est explicite si et seulement si la matrice $A$ est triangulaire strictement inférieure (pour un schéma DIRK la matrice $A$ est triangulaire inférieure, le caractère implicite est donnée par la diagonale). On ne précisera pas ici les propriétés liant $A$, $b$ et $c$ pour déterminer l'ordre de la méthode, son caractère SSP ou autre.

Dans l'obtention du polynôme caractéristique du schéma, le vecteur $b$ n'est pas nécessaire puisqu'on linéarise l'opération $L(t^n+b_j\Delta t , u^{(j)})$ par $\lambda u^{(j)}$.

Les 2 méthodes Runge-Kutta d'ordre 3 à 3 étages ont la même fonction de stabilité, par conséquent les CFL obtenues seront les mêmes et toutes les autres propriétés issues de cette fonction aussi (order-star). En revanche le schéma reste différent et l'erreur aussi (mais la convergence de l'erreur est bien la même, seule la constante diffère).

On peut souhaiter effectuer l'opération inverse, à partir d'une fonction de stabilité polynômiale (maximisant par exemple la stabilité sur l'axe imaginaire) retrouver le tableau de Butcher (et donc le schéma). Pour cela on calcule la fonction de stabilité d'une matrice $A=(a_{ij})_{ij}$ et un vecteur $c = (c_j)_j$, et on cherche à obtenir les relations d'ordre.

Maintenant il faut savoir quoi faire de ces relations d'ordre. On peut par exemple contraindre la forme de la matrice $A$ (par exemple diagonale inférieure, comme dans RK NSSP(5,3)), la positivité de certains coefficients, d'autres relations pour assurer un caractère SSP ou non, etc.

Il est normal que le résultat donné ici ne donne qu'une relation faisant intervenir $a_{10}$, $a_{20}$ et $a_{21}$ puisque pour le moment la fonction poly_butcher ne fonctionne qu'avec une méthode RK explicite. L'étude de DIRK2 par exemple nécessiterait une amélioration de la fonction.

On peut vouloir aussi écrire le schéma issu d'un tableau de Butcher. Par défaut le schéma obtenu directement à l'aide d'un tableau de Butcher n'est pas optimal d'un point de vue numérique dans le sens où il sous entend plus que $s$ (le nombre d'étages de la méthode) évaluations de la fonction $L$. Une amélioration du schéma est possible à l'aide d'un jeu de substituion des étapes précédentes.

2.1.3 Récapitulatif des domaines des différentes méthodes en temps

On peut résumer tous les schémas en temps présentés par la figure ci-dessous.

2.1.4 Recherche du meilleur domaine de stabilité

Nicolas est intéressé par effectuer une étude un peu inverse, c'est-à-dire trouver le meilleur coefficient $\alpha$ tel que le polynôme caractéristique :

$$ \alpha z^4 + \frac{z^3}{6} + \frac{z^2}{2} + z + 1 $$

ait le plus grand domaine de stabilité (en particulier sur l'axe imaginaire, pour exprimer la stabilité d'un schéma non diffusif). C'est à dire ajouter un étage sur une méthode RK3 existante, telle que celle-ci obtienne une meilleure stabilité.

Ce qui nous intéresse ici est de maximiser le domaine de stabilité sur l'axe imaginaire. Il est très compliqué de chercher à maximiser le domaine de stabilité d'une méthode RK($s$,$n$) pour WENO. En plus de cela il est nécessaire de trouver par la suite une solution aux relations d'ordre dans la fonction de stabilité.

La maximisation du domaine sur l'axe imaginaire permet de stabiliser avec la plus grande CFL possible un transport pur. Le coefficient $\alpha$ maximisant le domaine sur l'axe imaginaire avec un RK2 à 3 étages est $\alpha = \frac{1}{4}$.

Preuve : considérons le polynôme caractéristique $R_{\alpha}$ d'un schéma RK2 à 3 étages avec un coefficient $\alpha$ :

$$ R_{\alpha}(z) = 1 + z + \frac{z^2}{2} + \alpha z^3 $$

On cherche à évaluer le polynôme en $z=i\beta$, correspondant à l'intersection de l'ensemble $\{ z\ /\ |R_\alpha(z)|=1 \}$ avec l'axe imaginaire.

$$ R_\alpha(i\beta) = 1 + i\beta - \frac{\beta^2}{2} - i\beta^3\alpha $$

on souhaite un module égale à 1 :

$$ 1 = |R_\alpha(i\beta)|^2 = \left(1-\frac{\beta^2}{2}\right)^2 + \left(\beta - \beta^3\alpha \right)^2 $$

ce qui peut se simplifier par :

$$ \alpha^2\beta^6 + (\frac{1}{4}-2\alpha)\beta^4 = 0 $$

En supposant $\beta \neq 0$ :

$$ \alpha^2\beta^2 -2\alpha + \frac{1}{4} = 0 $$

soit :

$$ \beta^2 = \frac{2\alpha - \frac{1}{4}}{\alpha^2} $$

on obtient ainsi une fonction $\beta^2:\alpha\mapsto \frac{2\alpha - \frac{1}{4}}{\alpha^2}$ dont on souhaite trouver le maximum (ou plus exactement le coefficient $\alpha$ où est atteint le maximum, ce qui explique que l'on ne s'intéresse qu'au coefficient $\beta^2$ et non explicitement $\beta$).

$$ \left(\beta^2\right)'(\alpha) = \frac{-2\alpha^2 + \frac{\alpha}{2}}{\alpha^4} $$

la dérivée s'annule en $\alpha = \frac{1}{4}$ □

La maximisation du domaine sur l'axe imaginaire ne permet pas rigoureusement de maximiser la CFL du couple de schéma RK$^\alpha N$-WENO mais permet de se donner une idée ; prenons par exemple le cas de RK2 à 3 étages, comme dit précédemment le coefficient $\alpha$ maximisant le domaine de stabilité sur l'axe imaginaire est $\alpha = \frac{1}{4}$, représenté par le courbe orange sur la figure, or il est possible d'obtenir une CFL plus importante avec le coefficient $\alpha=\frac{1}{5}$, pour lequel on représente le coefficient d'amplification en rose de WENO$\times$1.4 (sans aller chercher $\alpha=\frac{1}{6}$ qui correspond à RK3).

2.1.5 Order star

Il est possible assez facilement, pour étudier certaines propriétés de nos schémas en temps, de tracer l'order star (n'ayant vu aucun article francophone en parler j'utlise le terme en anglais). L'order star est définie comme : $\{\mathcal{A}_+,\mathcal{A}_0,\mathcal{A}_-\}$ avec :

où $R$ est le polynôme caractéristique, ou approximation de Padé de la fonction à étudier.

Pour des raisons techniques, il est plus simple de colorier l'ensemble des points du plan dont la valeur d'une fonction est comprise entre 2 bornes finies, par conséquent, en bleu, nous alons afficher $ \mathcal{A}_-$.

2.2 Méthode diagonal implicit Runge-Kutta

Nous n'avons étudier jusqu'à présent que des schémas de type Runge-Kutta explicite (ou eRK), et cela se traduit par un tableau de Butcher triangulaire strictement inférieur. Le cas général d'un tableau de Butcher plein est envisageable, le schéma est alors implicite, ce qui permet souvent d'améliorer sa stabilité au détriment du coût de calcul (à cause de l'inversion d'un système). Le compromis de DIRK est de ne rendre que la diagonale implicite. Nous allons étudier ici ce genre de schéma.

Commençons par le premier schéma DIRK présenté dans [Alexander R. (1976)], dont voici le tableau de Butcher :

$$ \begin{array}{c c | c} \frac{1}{2}+\frac{1}{2\sqrt{3}} & 0 & \frac{1}{2} + \frac{1}{2\sqrt{3}} \\ -\frac{1}{\sqrt{3}} & \frac{1}{2}+\frac{1}{2\sqrt{3}} & \frac{1}{2} - \frac{1}{2\sqrt{3}} \\ \hline \frac{1}{2} & \frac{1}{2} & \end{array} $$

est une méthode DIRK(2,3). Le polynôme caractéristique des méthodes DIRK est une fonction rationnelle.

Par la suite, on ne présentera que les domaines de stabilité.

Il est intéressant de remarquer que :

  1. On trace l'opposé du coefficient d'amplification que nous avions précédemment car avec une méthode RK on résout $u_x = -u_x$, or on traçait le coefficient de $+u_x$.
  2. Il est important, surtout pour des méthodes implicites, de ne pas tracer que la bordure du domaine de stabilité car (on le verra par la suite) il est difficile de prédire de quel côté de la bordure est le domaine.

Posons $\alpha = \frac{2\cos\left(\frac{\pi}{18}\right)}{\sqrt{3}}$ : $$ \begin{array}{c c c | c} \frac{1+\alpha}{2} & 0 & 0 & \frac{1+\alpha}{2} \\ -\frac{\alpha}{2} & \frac{1+\alpha}{2} & 0 & \frac{1}{2} \\ 1+\alpha & -(1+2\alpha) & \frac{1+\alpha}{2} & \frac{1-\alpha}{2} \\ \hline \frac{1}{6\alpha^2} & 1-\frac{1}{3\alpha^2} & \frac{1}{6\alpha^2} & \end{array} $$ est une méthode DIRK(3,4).

Maintenant deux méthodes strong S-stable d'ordre 2 à 2 étages. $$ \begin{array}{c c | c } \alpha & 0 & \alpha \\ 1-\alpha & \alpha & 1 \\ \hline 1-\alpha & \alpha & \end{array} $$ avec $\alpha = 1\pm\frac{1}{2}\sqrt{2}$. Les propritétés de stabilités de la méthode changent significativement en fonction de la valeur du coefficient $\alpha$.

Maintenant une méthode strong S-stable d'ordre 3 à 3 étages :

$$ \begin{array}{c c c | c} \alpha & 0 & 0 & \alpha \\ \tau_2-\alpha & \alpha & 0 & \tau_2 \\ b_1 & b_2 & \alpha & 1 \\ \hline b_1 & b_2 & \alpha & \end{array} $$

avec $\alpha$ la racine de $x^3-3x+\frac{3}{2}x-\frac{1}{6}=0$ vivant dans $[\frac{1}{6},\frac{1}{2}]$, $\alpha \approx 0.43586652$

Un autre test issu de numipedia, où Ketcheson présente la méthode SSPIRK(3,3). Il est à noté que je ne me suis jamais intéressé au caractère SSP d'un schéma, donc je n'ai pas lu ou vérifié cet aspect sur les précédentes méthodes DIRK présentées.

Le tableau de Butcher de la méthode SSPIRK(3,3) est :

$$ \begin{array}{c c c | c} -\frac{\sqrt{2}}{4}+\frac{1}{2} & 0 & 0 & -\frac{\sqrt{2}}{4}+\frac{1}{2} \\ \frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4}+\frac{1}{2} & 0 & \frac{1}{2} \\ \frac{\sqrt{2}}{4} & \frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4}+\frac{1}{2} & \frac{\sqrt{2}}{4}+\frac{1}{2} \\ \hline \frac{1}{3} & \frac{1}{3} & \frac{1}{3} & \end{array} $$

L'objectif ici est de travailler avec un couple d'une méthode Runge-Kutta en temps et une méthode WENO en espace. Les méthodes implicites (ici uniquement sur la diagonale) nécessitent d'inverser la méthode WENO, ceci a fait l'objet d'une publication : [S. Gottlieb, J. Mullen et S. Ruuth]), il s'agit d'un problème compliqué et ces résultats ne semblent pas suffisamment convainquant pour poursuivre l'étude.

3. Étude de couples RK($s$,$n$)-WENO5

On regarde maintenant la stabilité de WENO avec les différents schémas RK($s$,$n$). Pour cela on trace le domaine de stabilité du schéma RK choisi, et on compare ce domaine avec celui du coefficient d'amplification du schéma WENO5. Le schéma est linéairement stable s'il existe un coefficient $\sigma$, appelé CFL, permettant de faire rentrer la courbe du schéma WENO dans le domaine de stabilité de RK($s$,$n$).

3.1 Automatisation du calcul de CFL

Il est maintenant intéressant d'étudier quel est le rapport d'homothétie liant le mieux RK($s$,$n$) à WENO, ce rapport d'homothétie est $\sigma$ le nombre de CFL. Dans l'article [Wang R. and Spiteri R. J. (2007)] une approximation grossière est faite en supposant que RK($s$,$n$), pour $n\geq 3$ contient un rectangle contenant la courbe de WENO linéarisé. Une estimation est ainsi faite d'un $\sigma_0 \leq \sigma$. Il est possible dans un premier temps d'estimer numériquement un coefficient $\sigma_{\texttt{num}}$ tel que :

$$ \sigma_0 \leq \sigma_{\texttt{num}} \leq \sigma $$

travail effectué dans [Motamed M. and Macdonald C. B. (2010)]. Il est aussi possible de travailler sur la version complète du schéma WENO, en prenant en compte dans l'analyse de von Neumann de la partie non linéaire.

Pour estimer une CFL entre un couple RK($s$,$n$)-WENO(lin) :

  1. On discrétise la frontière du domaine de stabilité de la méthode RK($s$,$n$), chose déjà effectuée par Python, on obtient ainsi un tableau de $\{r(\theta_k)\}_k$.
  2. On discrétise le coefficient d'amplification en fonction de $\phi$ en discrétisant l'intervalle $[0,2\pi[$, on obtient ainsi un tableau des coefficients d'amplification $\{w(\phi_i)\}_i$. Il est à noter qu'il n'y a pas de raison que la discrétisation $\theta_k$ coïncide avec les $\phi_i$.
  3. On trie ces 2 tableaux de données par ordre d'argument, cela sera plus simple pour les parcourir de manière similaire.
  4. Pour un couple $(\varphi_i,w(\phi_i)$ donné issu de la discrétisation, on recherche l'élément $\rho(\phi)$ défini comme : $$\rho(\phi_i) = \arg\!\min_{r(\theta_k)}\left(arg(r(\theta_k)-w(\phi_i)\right)$$
  5. On calcule le facteur d'homothétie locale : $$\sigma_i = \left| \frac{\rho(\phi_i)}{w(\phi_i)}\right|$$
  6. On définit la CFL comme : $$\sigma = \min_i \sigma_i$$

L'étape 4 de l'alorithme peut se représenter comme suit, ici avec la méthode RK(4,4).

Pour chaque point $w(\phi_i)$ du coefficient d'amplification de WENO linéarisé, on souhaite trouver l'élément sur la bordure du domaine de stabilité de la méthode Runge-Kutta avec l'argument le plus proche de $\arg(w(\phi_i))$. Ici par exemple on trouve deux éléments $r(\theta_k)$ et $r(\theta_{k+1})$ avec un argument semblable, on choisit celui avec l'argument le plus proche.

Si on souhaite avec un résultat plus précis, il est nécessaire d'effectuer une interpolation entre les 2 éléments de la bordure du domaine de stabilité de RK(4,4) que l'on trouve. On se limite ici à une interpolation constante par morceau, la précision ne sera approtée que par le rafinement du maillage angulaire.

3.1.1 CFL RK($s$,$n$) - WENO linéarisé

On commence par expérimenter notre méthode algorithmique sur le cas WENO linéarisé et retrouver les CFL proposées dans [Motamed M. and Macdonald C. B. (2010)] ou [Thibaut L. et al.].

On confirme bien avec ce résultat l'aspect instable du couple EE-WENO5 démontré par dans [R. Wang and R. J. Spiteri (2007)].

On peut aussi calculer des CFL avec des couples RK($s$,$n$)-WENO3.

L'obtention d'une CFL non nulle avec la méthode d'Euler explicite demande peut-être quelques vérification. La discrétisation a été augmenté juste pour ce cas pour vérifier cette prétendue stabilité, mais le script semble indiquer une certaine stabilité. Une étude approfondie à l'aide du vecteur directeur en $(0,0)$ (calculée à partir d'un développement limité des schémas) est peut-être nécessaire.

3.1.2 CFL RK($s$,$n$) - up-wind

On peut aussi effectuer ce travail pour le schéma up-wind. La méthode est exactement la même, mais à partir du cercle de valeurs propres de up-wind.

Étrangement, même en touchant aux paramètres numériques, on ne trouve pas 1. Cela doit provenir de la discrétisation en amont de upwind dans la première partie de ce document.

3.1.3 CFL RK($s$,$n$) - CD2

Pour le schéma CD2, la méthode précédente ne fonctionne pas. La méthode se concentre sur l'étude sur l'axe imaginaire des différents intégrateurs en temps et de trouver la valeur $y_\text{max}$ sur l'axe imaginaire telle que : $$ \forall y\in[-y_\text{max},y_\text{max}], |p_{RK(s,n)}(iy)| \leq 1 $$ où $p_{RK(s,n)}$ est la fonction de stabilité de la méthode de Runge-Kutta considérée.

3.1.4 Récapitulatif des CFL obtenues

On récapitule les CFL obtenues pour le schéma WENO5 en les comparant à la littérature :

- Euler RK(3,3) RK(5,3) RK(4,4)
[Thibaut L. et al.] 0.00 1.44 2.14 1.73
WENO5 linéarisé 0.00 1.430 2.558 1.730

Je n'arrive pas à expliquer la différence sur le schéma RK(5,3) entre [Thibaut L. et al.] et ma méthode de calcul de CFL avec WENO linéarisé.

Ces CFL s'appliquent à une équation de transport 1D du type :

$$ u_t + au_x = 0 $$

La CFL d'un couple de schéma RK($s$,$n$)-WENO5 pour résoudre une telle équation est de la forme :

$$ \Delta t \leq \frac{\sigma \Delta x}{a} $$

avec $\sigma$ donné par le tableau précédent.

3.2 Constante d'erreur, couple idéal

Il est possible de remarquer une constante d'erreur lors que l'on trace les courbes d'ordre entre une méthode RK($s$,$n$) et une méthode d'intégration en espace. Nous nous proposons dans cette section de l'étudier numériquement. Pour cela dans un premier temps nous allons écrire une fonction qui prend une méthode Runge-Kutta (son tableau de Butcher) et retourne une fonction python appelable pour résoudre (sur une itération) une équation du type $\dot{u} = L(t,u)$, dans la pratique on prendra $L:t,u\mapsto -a\partial_xu$ avec $a\in\mathbb{R}$ la vitesse d'advection.

On effectue maintenant un test sur une équation de transport :

$$ u_t + u_x = 0 $$

avec $u(t=0,x)=u_0(x) = \cos(2\pi x)$, avec 20 points de discrétisation en $x$ sur le segement $[0,1]$, $\Delta t = \Delta x$, jusqu'au temps $T_f=1.2$. Comme discrétisation en temps nous utiliserons les méthodes données par rk33, rk33bis et rknssp33, en espace on utilisera la méthode WENO5.

Sur ce cas test les résultats sont confondus. Cette condition initiale est trop régulière, nous allons donc effectuer un test, en temps plus long avec comme condition initiale un test de Shu. Nous savons que l'erreur est en $\mathcal{O}(\Delta x^p)+\mathcal{O}(\Delta t^q)$ avec $p$ l'ordre du schéma en espace (ici $p=5$) et $q$ l'ordre du schéma en temps ($q=3$ sur ce premier exemple). Nous souhaitons maximiser l'erreur en temps pour que celle-ci soit la plus visible, nous nous placerons un peu en dessous de la CFL, pour un couple RK(3,3)-WENO5 celle-ci est de $\sigma=1.43$. Nous allons aussi effectuer une simulation en temps plus long $T_f=10$.

Maintenant que l'on voit cette constante d'erreur sur un cas test présentant de fortes discontinuités (de la fonction et de sa dérivée), essayons de mesurer l'ordre sur ce cas.

On peut remarquer cette constante d'erreur entre les schémas pour les plus grands $\Delta t$. Ces simulations s'effectuant en temps long sur le transport d'une condition initiale non régulière il est normal de ne pas retrouver l'ordre 3 de ces méthodes. Pour de plus petits pas de temps on remarques que l'erreur des 3 méthodes convergent vers la même valeur, avec la même pente.

Le calcul de la CFL de couple RK($s$,$n$)-WENO permet de connaître les paramètres assurants la stabilité du système. On souhaite ici ce mettre dans un régime avec le plus grand $\Delta t$ et $\Delta x$ possible en assurant la stabilité et la plus faible erreur de convergence possible. On voit que le schéma RK(3,3) bis correspond mieux à nos attentes.

4. Équation avec un terme linéaire et un non linéaire

Regardons le problème :

$$ \dot{u} = Lu + N(u) $$

où $N$ est une fonction, non linéaire, de $u$. On peut vouloir construire une formule de Duhamel à partir de cette équation, en remarquant que :

$$ \partial_t\left(e^{-Lt}u\right) = e^{-Lt}N(u(t)) $$

Pour construire un schéma en temps à partir de cette équation il faut intégrer sur $[t^n,t^{n+1})$ :

$$ \int_{t^n}^{t^{n+1}} \partial_t\left(e^{-Lt}u\right)\,\mathrm{d}t = \int_{t^n}^{t^{n+1}} e^{-Lt}N(u(t))\,\mathrm{d}t $$

On évalue de manière exacte le terme de gauche :

$$ e^{-Lt^n}\left( e^{-L\Delta t}u(t^{n+1}) - u(t^n) \right) = \int_{t^n}^{t^{n+1}} e^{-Lt}N(u(t))\,\mathrm{d}t $$

soit :

$$ e^{-L\Delta t}u(t^{n+1}) - u(t^n) = \int_{t^n}^{t^{n+1}} e^{-L(t-t^n)}N(u(t))\,\mathrm{d}t $$

En effectuant le changement de variable $\tau = t-t^n$ on retrouve une intégrale sur le seul pas de temps $\Delta t$, ce qui est rassurant quant à la stabilié numérique en temps long (à propos des possibles erreurs d'overflow ou d'underflow produites par la fonction exponentielle).

On souhaite construire un schéma en temps, donc on réécrit le problème sous la forme :

$$ u(t^{n+1}) = e^{L\Delta t}u(t^n) + \int_{0}^{\Delta t} e^{-L(s-\Delta t)}N(u(t^n+s))\,\mathrm{d}s $$

Nous allons ici présenter deux classes de schéams dit exponentiels, pour résoudre cette intégrale.

4.1 Schéma de Lawson

On effectue le changement de variable $v(t) = e^{-Lt}u(t)$, ce qui permet de réécrire l'équation $\dot{u} = Lu+N(u)$ en :

$$ \dot{v} = e^{-Lt}N(e^{Lt}v(t)) $$

en posant $\tilde{N}:(t,v)\mapsto e^{-Lt}N(e^{Lt}v(t))$, on obtient :

$$ \dot{v} = \tilde{N}(t,v) $$

qui peut se résoudre avec une méthode de Runge-Kutta, que l'on réécrira en fonction de $u$ et $L$ pour obtenir un schéma dit de Lawson. Prenons par exemple le schéma RK SSP(3,3) dit de Shu-Osher :

$$ \begin{aligned} v^{(1)} &= v^n + \Delta t \tilde{N}(t^n,v^n) \\ v^{(2)} &= \frac{3}{4}v^n + \frac{1}{4}v^{(1)} + \frac{1}{4}\Delta t \tilde{N}(t^n+\Delta t,v^{(1)}) \\ v^{n+1} &= \frac{1}{3}v^n + \frac{2}{3}v^{(2)} + \frac{2}{3}\Delta t \tilde{N}(t^n+\frac{1}{2}\Delta t,v^{(2)}) \\ \end{aligned} $$

Ce qui nous donne en fonction de $u$ :

$$ \begin{aligned} v^{(1)} &= e^{-L t^n}u^n + \Delta t e^{-L t^n}N(u^n) \\ v^{(2)} &= \frac{3}{4}e^{-L t^n}u^n + \frac{1}{2}v^{(1)} + \frac{1}{4}\Delta t e^{-Lt^n}e^{-L\Delta t}N(e^{Lt^n}e^{L\Delta t}v^{(1)}) \\ e^{-L t^n}e^{-L\Delta t}u^{n+1} &= \frac{1}{3}e^{-L t^n}u^n + \frac{2}{3}v^{(2)} + \frac{2}{3}\Delta t e^{-Lt^n}e^{-\frac{1}{2}L\Delta t}N(e^{Lt^n}e^{\frac{1}{2}L\Delta t}v^{(2)}) \end{aligned} $$

Pour des raisons purement numériques liées à la précision machine, simplifions dès que possible par $e^{-Lt^n}$ pour s'assurer que le schéma ne dégénère pas en temps long :

$$ \begin{aligned} \tilde{u}^{(1)} &= u^n + \Delta t N(u^n) \\ \tilde{u}^{(2)} &= \frac{3}{4}u^n + \frac{1}{2}\tilde{u}^{(1)} + \frac{1}{4}\Delta t e^{-L\Delta t}N(e^{L\Delta t}\tilde{u}^{(1)}) \\ u^{n+1} &= \frac{1}{3}e^{L\Delta t}u^n + \frac{2}{3}e^{L\Delta t}\tilde{u}^{(2)} + \frac{2}{3}\Delta t e^{\frac{1}{2}L\Delta t}N(e^{\frac{1}{2}L\Delta t}\tilde{u}^{(2)}) \end{aligned} $$

Les coefficients exponentiels à l'intérieur de l'opérateur non-linéaire nous incitent à poser :

D'où le schéma suivant :

$$ \begin{aligned} u^{(1)} &= e^{L\Delta t}u^n + \Delta t e^{L\Delta t} N(u^n) \\ u^{(2)} &= \frac{3}{4}e^{\frac{1}{2}L\Delta t}u^n + \frac{1}{4}e^{-\frac{1}{2}L\Delta t}u^{(1)} + \frac{1}{4}\Delta t e^{-\frac{1}{2}L\Delta t}N(u^{(1)}) \\ u^{n+1} &= \frac{1}{3}e^{L\Delta t}u^n + \frac{2}{3}e^{\frac{1}{2}L\Delta t}u^{(2)} + \frac{2}{3}\Delta t e^{\frac{1}{2}L\Delta t}N(u^{(2)}) \end{aligned} $$

On retrouve ainsi le schéma proposé dans [Isherwood L. et al (2018)].pdf). Il est également possible de généralisé ce calcul au cas où l'opérateur non linéaire $N$ dépend aussi du temps, on retrouve alors les mêmes temps $t^n+c_i\Delta t$ que pour la méthode Runge-Kutta sous-jascente.

On souhaite déterminer le polynôme de stabilité de cette méthode de Lawson induite par la méthode RK(3,3) de Shu-Osher. Pour cela on linéarise le problème, et on pose $N:u\mapsto \mu u$, et pour conserver des notations similaires à l'étude des schémas de Runge-Kutta on utilisera $\lambda$ comme terme linéaire, au lieu de $L$.

On retrouve le polynôme caractéristique d'un schéma RK SSP(3,3) classique multiplié par $e^{\lambda\Delta t}$. Par conséquent pour un $\lambda$ imaginaire pur, le module de son exponentiel est égal à 1, donc le domaine de stabilité reste inchangé. On remarque que le cas général $\lambda\in\mathbb{C}$ ou le cas particulier $\lambda\in\mathbb{R}$ (intéressant pour la prise en compte d'un terme de collision BGK par exemple) est plus délicat car le module du polynôme est alors multiplié par $e^{\Delta t \Re(\lambda)}$, qui dépend de $\Delta t$, il ne s'agit donc pas d'une simple homothétie du domaine de stabilité.

Le cas $\lambda\in i\mathbb{R}$ est celui qui nous intéressera le plus par la suite, à cause de l'utilisation de méthode spectrale dans une direction et d'une méthode WENO dans l'autre. Ce cas ne changeant pas le domaine de stabilité, il n'est pas nécessaire d'effectuer de nouvelle étude sur la stabilité d'une telle méthode.

Il est possible d'automatiser la méthode d'obtention du schéma de Lawson (ou IFRK) pour la résolution de l'équation :

$$ \partial_t u = Lu + N(u) $$

à partir d'un tableau de Butcher de la méthode RK associée. La méthode algorithmique mise en place consiste à :

  1. Écrire le schéma RK de base en $v=e^{-Lt}u$ et $\tilde{N}(t,v)=e^{-Lt}N(e^{Lt}v)$, on obtient ainsi une méthode de résolution du problème :

    $$ \partial_t v = \tilde{N}(t,v) $$

  2. Substituer :

    • $v^n \mapsto e^{-Lt^n}u^n$
    • $v^{n+1} \mapsto e^{-L(t^n+\Delta t)}u^{n+1}$
    • $v^{(i)} \mapsto e^{-L(t^n+c_i\Delta t)}u^{(i)}$

Il est possible d'effectuer ce même travail à partir de l'écriture sous forme de schéma, mais il est plus compliqué d'un point de vue informatique de gérer des expressions plutôt que des tableaux. Il s'agit donc uniquement d'une contrainte technique.

On remarquera que le schéma ainsi obtenu n'est pas optimal dans le sens où il y a plus que $s$ appels à la fonction non linéaire $N$, chose évitée souvent dans l'écriture sous forme de schéma, comme c'est le cas pour la méthode RK3 de Shu-Osher. Il est donc possible de substituer récurisvement, grâce aux équations précédentes, les valeurs des différentes évaluations de la fonction $N$.

On peut limiter le nombre d'appel à la fonction $N$, souvent non-linéaire et coûteuse en temps de calcul, en shuosherisant une méthode de Lawson. En contre partie l'utilisation en espace mémoire est plus important puisqu'il est souvent nécessaire de stocker les valeurs des étages intermédiaires. L'idéal pour effectuer du calcul est donc une méthode sous-diagonale comme la méthode RK(5,3) précédemment présentée.

On retrouve bien le schéma de Lawson précédemment donné. On peut ainsi le calculer pour toute méthode de Runge-Kutta, par exemple pour DP5.

De manière plus générale, pour une méthode de Lawson à 4 étages, voici le schéma de Lawson associé.

Étude du ratio entre la fonction de stabilité d'une méthode de Lawson et la méthode Runge-Kutta dont elle est induite.

On peut démontrer ce résultat. Pour cela partons du problème de base : $$\dot{u} = Lu + N(u)$$ que l'on réécrit avec une formule de Duhamel comme : $$\partial_t\left(e^{-Lt}u\right) = e^{-Lt}N(u)$$ On applique le changement de variable habituel pour les schémas de Lawson : $v = e^{-Lt}u$ et on pose $\tilde{N}:(t,v)\mapsto e^{-Lt}N(e^{Lt}v)$, le problème s'écrit alors sous la forme canonique de la résolution d'une méthode Runge-Kutta : $$\dot{v} = \tilde{N}(t,v)$$

On écrit notre schéma type Runge-Kutta RK($s$,$n$) sur $v$, avec $\tilde{N}$ une fonction linéaire. Il est à noter que la variable de linéarisation de $N$ et de $\tilde{N}$ est la même puisqu'en linéarisant $\tilde{N}$ on obtient $e^{-Lt}e^{Lt}Nv = Nv$. Par définition de la fonction de stabilité, pour une fonction $\tilde{N}$ (ou $N$) linéaire, le schéma Runge-Kutta peut s'écrire : $$v^{n+1} = p_{(s,n)}(\Delta t N)v^n$$

En revenant à la variable $u$ cela se réécrit : $$ e^{-Lt^n}e^{-L\Delta t}u^{n+1} = p_{(s,n)}(\Delta tN)e^{Lt^n}u^n$$ soit : $$u^{n+1} = p_{(s,n)}(\Delta tN)e^{L\Delta t}u^n$$

La fonction de stabilité de la méthode de Lawson induite par la méthode RK($s$,$n$) est : $$p_{L(s,n)} = p_{(s,n)}e^{L\Delta t}$$ □

Dans notre cadre d'utilisation, $L$ est un imaginaire pur, par conséquent le domaine de stabilité de la méthode de Lawson est le même : $$ \mathcal{D}_{(s,n)} = \left\{z\in\mathbb{C}: |p_{(s,n)}(z)|\leq 1 \right\} $$

4.2 Schémas exponentiels

On présente ici une autre classe de schéma pour résoudre :

$$ u(t^{n+1}) = e^{L\Delta t}u(t^n) + \int_{0}^{\Delta t} e^{-L(s-\Delta t)}N(u(t^n+s))\,\mathrm{d}s $$

Pour cela, les méthodes exponentielles proposent de résoudre de manière exacte l'exponentielle, en approximant $u(t) \approx u(t^n)$ (ou en effectuant une quadrature pour l'ordre supérieur). On obtient alors :

$$ e^{-\lambda\Delta t}u^{n+1} - u^n = \frac{\mu}{\lambda} u^n (1-e^{-\lambda\Delta t}) $$

soit :

$$ u^{n+1} = e^{\lambda\Delta t}u^n + \frac{\mu}{\lambda} u^n (e^{\lambda\Delta t}-1) $$

Pour étudier la stabilité du schéma on se propose de poser :

On obtient alors l'expression suivante :

$$ u^{n+1} = \left(e^{z_1} + \frac{z_2}{z_1}(e^{z_1}-1)\right)u^n $$

On se propose maintenant d'étudier le module du coefficient d'amplification. En supposant que l'opérateur représenté par $z_1$ est peu visqueux, voir pas du tout dans le cas d'une méthode spectrale, c'est-à-dire que $\Re(z_1) \approx 0$, on écrit alors $z_1 = ib_1$.

$$ R(b_1,z_2) = \left| e^{ib_1} - i\frac{z_2}{b_1}(e^{ib_1}-1) \right| $$

En factorisant par $e^{ib_1}$, qui est de module 1, on obtient :

$$ R(b_1,z_2) = \left| 1 + i\frac{z_2}{b_1}(e^{-ib_1} - 1) \right| $$

On écrit maintenant $z_2 = a_2 + ib_2$ avec $a_2,b_2\in\mathbb{R}$, le module peut alors s'exprimer explicitement :

$$ R(b_1,a_2,b_2)^2 = \left| 1 + i\frac{a_2+ib_2}{b_1}(\cos(b_1) - i\sin(b_1) - 1) \right|^2 $$

On remarque que le domaine de stabilité évolue en fonction de la valeur de $\lambda = i\texttt{b1}$, ce qui peut paraître normal puisque l'on fait rentrer $e^{-\lambda t}$ dans la dérivée temporelle.

Nicolas a poussé l'étude des schémas Runge-Kuta exponentiel aux ordres plus élevés, et on obtient des domaines de stabilités très étranges et surtout asymétriques. C'est l’asymétrie qui nous pose le plus de problème pour l'obtention d'une CFL optimal avec un schéma spatial d'ordre élevé (style WENO).

Les schémas exponentiels peuvent aussi s'écrire sous la forme d'un tableau de Butcher légèrement modifié. (Pour me simplifier la vie je reprends les notations de Lukas) On s'intéresse à la résolution d'un problème du type :

$$ \dot{u} + Au = g(t,u) $$

On a le tableau de Butcher suivant :

$$ \begin{array}{c|c} \begin{matrix} c_1 \\ \vdots \\ c_s \end{matrix} & \begin{matrix} a_{11} & \cdots & a_{1s} \\ \vdots & \ddots & \vdots \\ a_{s1} & \cdots & a_{ss} \end{matrix} \\ \hline & \begin{matrix}b_1 & \cdots & b_s \end{matrix} \\ \end{array} $$

Attention : les notations $b$ et $c$ sont inversées par rapport à la présentation précédente des tableaux de Butcher.

Nota bene : nous ne nous intéressons qu'à des méthodes explicites (à cause de WENO) donc la matrice $(a)_{i,j}$ est triangulaire strictement inférieure.

Le schéma associé à ce tableau s'écrit comme :

$$ \begin{aligned} G_{nj} &= g(t^n+c_jh,U_{nj}) \\ U_{ni} &= \chi_i(-hA)u_n + h\sum_{j=1}^{s}a_{ij}(-hA)G_{nj} \\ u^{n+1} &= \chi(-hA)u_n + h\sum_{i=1}^{s}b_i(-hA)G_{ni} \end{aligned} $$

avec les fonctions :

Les coefficients $a_{i,j}$ et $b_i$ sont des fonctions pouvant faire intervenir les fonctions :

$h$ représente le pas de temps dans notre cas, $A$ est la fonction linéaire, dans le cas de Vlasov-Poisson après FFT en $x$ on a $A=i\kappa v$.

L'order star de expRK2 indique que pour $a=0$ (courbe bleue) on a une méthode d'ordre 2, mais l'ordre n'est plus que de 0 pour $a \neq 0$.

L'étude des schémas exponentiels Runge-Kutta se concentrera sur la recherche du paramètre $y_\text{max}$ tel que : $$i[-y_\text{max},y_\text{max}]\subset\mathcal{D}_{(s,n)}^a,\,\forall a$$