Crank-Nicolson est-il un schéma de discrétisation stable pour l'équation de réaction-diffusion-advection (convection)?

Je ne connais pas très bien les schémas de discrétisation courants pour les PDE. Je sais que Crank-Nicolson est un schéma populaire pour discrétiser l'équation de diffusion. Est également un bon choix pour le terme d'advection?

Je suis intéressé par la résolution de l'équation Réaction-Diffusion-Advection ,

$\frac{\partial u}{\partial t} + \nabla \cdot \left( \boldsymbol{v} u - D\nabla u \right) = f$

où $D$ est le coefficient de diffusion de la substance $u$ et $\boldsymbol{v}$ est la vitesse.

Pour mon application spécifique, l'équation peut être écrite sous la forme,

$\frac{\partial u}{\partial t} = \underbrace{D\frac{\partial^2 u}{\partial x^2}}_{\textrm{Diffusion}} + \underbrace{\boldsymbol{v}\frac{\partial u}{\partial x}}_{\textrm{Advection (convection)}} + \underbrace{f(x,t)}_{\textrm{Reaction}}$

Voici le schéma Crank-Nicolson que j'ai appliqué,

$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( u_{j-1}^{n} - 2u_{j}^{n} + u_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( u_{j-1}^{n+1} - 2u_{j}^{n+1} + u_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\alpha}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\alpha}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + f(x,t)$

Remarquez les termes $\alpha$ et $\beta$ . Cela permet au schéma de se déplacer entre:

• $\beta=\alpha=1/2$ Crank-Niscolson,
• $\beta=\alpha=1$ c'est totalement implicite
• $\beta=\alpha=0$ c'est totalement explicite

Les valeurs peuvent être différentes, ce qui permet au terme de diffusion d'être Crank-Nicolson et au terme d'advection d'être autre chose. Quelle est l'approche la plus stable, que recommanderiez-vous?

C'est une question bien cadrée et une chose très utile à comprendre. Korrok a raison de vous référer à l'analyse de von Neumann et au livre de LeVeque. Je peux ajouter un peu plus à cela. Je voudrais écrire une réponse détaillée, mais pour le moment je n'ai que le temps pour une courte:

Avec , vous obtenez une méthode qui est absolument stable pour des tailles de pas arbitrairement grandes, ainsi qu'une précision de second ordre. Cependant, la méthode n'est pas stable en L , donc les très hautes fréquences ne seront pas amorties, ce qui n'est pas physique.$\alpha=\beta=1/2$

Avec , vous obtenez une méthode qui est également inconditionnellement stable, mais précise uniquement au premier ordre. Cette méthode est très dissipative. Il est stable en L.$\alpha=\beta=1$

Si vous prenez , votre méthode peut être comprise comme l'application d'une méthode Runge-Kutta additive à la semi-discrétisation à différence centrée. L'analyse de stabilité et de précision de ces méthodes est considérablement plus compliquée. Un très bel article sur de telles méthodes est ici .$\alpha\ne\beta$

L'approche à recommander dépend fortement de l'ampleur de , du type de données initiales que vous traitez et de la précision que vous recherchez. Si une précision très faible est acceptable, alors est une approche très robuste. Si est modéré ou grand, le problème est dominé par la diffusion et très rigide; généralement donnera de bons résultats. Si est très petit, il peut être avantageux d'utiliser une méthode explicite et un remontage d'ordre supérieur pour les termes convectifs.α = β = 1 D α = β = 1 / 2 $D$$\alpha=\beta=1$$D$$\alpha=\beta=1/2$$D$

De manière générale, vous voudrez utiliser une méthode implicite pour les équations paraboliques (la partie diffusion) - les schémas explicites pour la PDE parabolique doivent avoir un pas de temps très court pour être stable. Inversement, pour la partie hyperbolique (advection), vous aurez besoin d'une méthode explicite car elle est moins chère et ne perturbe pas la symétrie du système linéaire que vous devez résoudre en utilisant un schéma implicite de diffusion. Dans ce cas, vous voulez éviter les différences centrées comme et passer aux différences unilatérales pour des raisons de stabilité.$(u_{j+1}-u_{j-1})/2\Delta t$$(u_j-u_{j-1})/\Delta t$

Je vous suggère de consulter le livre de Randy Leveque ou le livre de Dale Durran pour "l'analyse de stabilité de von Neumann". C'est une approche générale pour déterminer la stabilité de votre schéma de discrétisation, à condition que vous ayez des conditions aux limites périodiques. (Il y a aussi un bon article wiki ici .)

L'idée de base est de supposer que votre approximation discrète peut s'écrire une somme d'ondes planes , où est le nombre d'onde et la fréquence. Vous entassez une onde plane dans votre approximation du PDE et priez pour qu'elle ne explose pas. Nous pouvons réécrire l'onde plane en et nous voulons nous assurer que .$e^{i(k j\Delta x-\omega n\Delta t)}$$k$$\omega$$\xi^n e^{ikj\Delta x}$$|\xi|\le 1$

À titre d'illustration, considérons l'équation de diffusion ordinaire avec une différenciation totalement implicite:

$\frac{u^{n+1}_j-u^n_j}{\Delta t} = D\frac{u^{n+1}_{j-1}-2u^{n+1}_j+u^{n+1}_{j+1}}{\Delta x^2}$

Si nous substituons dans une onde plane, puis divisons par et , nous obtenons l'équation$\xi^n$$e^{ikj\Delta x}$

$\frac{\xi-1}{\Delta t} = D\frac{e^{-ik\Delta x}-2+e^{ik\Delta x}}{\Delta x^2}\xi$

Nettoyez-le un peu maintenant et nous obtenons:

$\xi = \frac{1}{1+\frac{2D\Delta t}{\Delta x^2}\left(1-\cos k\Delta x\right)}$ .

C'est toujours moins d'un, donc vous êtes en clair. Essayez d'appliquer ceci pour le schéma explicite et centré de l'équation d'advection:

$\frac{u_j^{n+1}-u_j^n}{\Delta t} = v\frac{u_{j-1}^n-u_{j+1}^n}{2\Delta x}$

et voir ce que vous. (Il y aura une partie imaginaire cette fois.) Vous trouverez que , ce qui est triste. D'où mon avertissement que vous ne l'utilisez pas. Si vous pouvez le faire, vous ne devriez pas avoir beaucoup de mal à trouver un schéma stable pour l'équation d'advection-diffusion complète.$\xi$$|\xi|^2 > 1$

Cela dit, j'utiliserais un schéma totalement implicite pour la partie diffusion. Remplacez la différenciation dans la partie advection par si et si et choisissez un pas de temps pour que . (Il s'agit de la condition Courant-Friedrichs-Lewy .) Elle n'est précise qu'au premier ordre, vous pouvez donc rechercher des schémas de discrétisation d'ordre supérieur si cela vous concerne. v > 0 u j - u j + 1 v < 0 V Δ t / Δ x ≤ $u_j-u_{j-1}$$v > 0$$u_j-u_{j+1}$$v < 0$$V\Delta t/\Delta x \le 1$