¿Es Crank-Nicolson un esquema de discretización estable para la ecuación Reacción-Difusión-Advección (convección)?

No estoy muy familiarizado con los esquemas comunes de discretización para PDEs. Sé que Crank-Nicolson es un esquema popular para discretizar la ecuación de difusión. ¿También es una buena opción para el término de advección?

Estoy interesado en resolver la ecuación Reacción-Difusión-Advección ,

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

donde $D$ es el coeficiente de difusión de la sustancia $u$ y $\boldsymbol{v}$ es la velocidad.

Para mi aplicación específica, la ecuación se puede escribir en la forma,

$\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}}$

Aquí está el esquema de Crank-Nicolson que he aplicado,

$\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)$

Observe los términos $\alpha$ y $\beta$ . Esto permite que el esquema se mueva entre:

• $\beta=\alpha=1/2$ 1/2 Crank-Niscolson,
• $\beta=\alpha=1$ es totalmente implícito
• $\beta=\alpha=0$ es completamente explícito

Los valores pueden ser diferentes, lo que permite que el término de difusión sea Crank-Nicolson y que el término de advección sea otra cosa. ¿Cuál es el enfoque más estable, qué recomendarías?

Esta es una pregunta bien enmarcada y una cosa muy útil de entender. Korrok tiene razón al referirlo al análisis de von Neumann y al libro de LeVeque. Puedo agregar un poco más a eso. Me gustaría escribir una respuesta detallada, pero por el momento solo tengo tiempo para una breve:

Con , obtienes un método que es absolutamente estable para tamaños de paso arbitrariamente grandes, así como también una precisión de segundo orden. Sin embargo, el método no es estable en L , por lo que no se amortiguarán las frecuencias muy altas, lo que no es físico.$\alpha=\beta=1/2$

Con , obtienes un método que también es incondicionalmente estable, pero solo preciso de primer orden. Este método es muy disipativo. Es estable en L.$\alpha=\beta=1$

Si toma , su método puede entenderse como la aplicación de un método aditivo Runge-Kutta a la semi-discretización de diferencia centrada. El análisis de estabilidad y precisión para tales métodos es considerablemente más complicado. Un documento muy bueno sobre tales métodos está aquí .$\alpha\ne\beta$

El enfoque a recomendar depende en gran medida de la magnitud de , el tipo de datos iniciales con los que trata y la precisión que busca. Si se acepta una precisión muy baja, entonces es un enfoque muy sólido. Si es moderado o grande, entonces el problema está dominado por la difusión y es muy rígido; típicamente dará buenos resultados. Si es muy pequeño, entonces puede ser ventajoso usar un método explícito y un orden ascendente de orden superior para los términos convectivos.$D$$\alpha=\beta=1$$D$$\alpha=\beta=1/2$$D$

En términos generales, querrá utilizar un método implícito para las ecuaciones parabólicas (la parte de difusión): los esquemas explícitos para la PDE parabólica deben tener un tiempo muy breve para ser estables. Por el contrario, para la parte hiperbólica (advección) querrá un método explícito, ya que es más barato y no interrumpe la simetría del sistema lineal que debe resolver utilizando un esquema implícito para la difusión. En ese caso, desea evitar diferencias centradas como y cambiar a diferencias unilaterales por razones de estabilidad.$(u_{j+1}-u_{j-1})/2\Delta t$$(u_j-u_{j-1})/\Delta t$

Le sugiero que mire el libro de Randy Leveque o el libro de Dale Durran para el "análisis de estabilidad von Neumann". Es un enfoque general para determinar la estabilidad de su esquema de discretización, siempre que tenga condiciones límite periódicas. (También hay un buen artículo wiki aquí ).

La idea básica es suponer que su aproximación discreta se puede escribir una suma de ondas planas , donde es el número de onda y la frecuencia. Aprietas una onda plana en tu aproximación al PDE y rezas para que no explote. Podemos reescribir la onda plana como y queremos asegurarnos de que .$e^{i(k j\Delta x-\omega n\Delta t)}$$k$$\omega$$\xi^n e^{ikj\Delta x}$$|\xi|\le 1$

A modo de ilustración, considere la ecuación de difusión ordinaria con diferenciación totalmente implícita:

$\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 sustituimos en una onda plana, luego dividimos por y , obtenemos la ecuación$\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$

Limpie esto un poco ahora y obtenemos:

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

Esto siempre es menos de uno, por lo que está en claro. Intente aplicar esto para el esquema explícito y centrado de la ecuación de advección:

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

y mira qué obtienes. (Tendrá una parte imaginaria esta vez). Encontrarás que , que son momentos tristes. De ahí mi advertencia de que no lo uses. Si puede hacer eso, entonces no debería tener muchos problemas para encontrar un esquema estable para la ecuación de advección-difusión completa.$\xi$$|\xi|^2 > 1$

Dicho esto, usaría un esquema totalmente implícito para la parte de difusión. Cambie la diferencia en la parte de advección a si y si y elija un paso de tiempo para que . (Esta es la condición de Courant-Friedrichs-Lewy ). Solo es precisa de primer orden, por lo que es posible que desee buscar esquemas de discretización de orden superior si eso le preocupa.$u_j-u_{j-1}$$v > 0$$u_j-u_{j+1}$$v < 0$$V\Delta t/\Delta x \le 1$