O Crank-Nicolson é um esquema estável de discretização para a equação de reação à difusão-advecção (convecção)?

Não estou muito familiarizado com os esquemas comuns de discretização para PDEs. Eu sei que Crank-Nicolson é um esquema popular para discretizar a equação de difusão. Também é uma boa escolha para o termo de advecção?

Eu sou interessante em resolver a equação Reação-Difusão-Advecção ,

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

onde $D$ é o coeficiente de difusão da substância $u$ e $\boldsymbol{v}$ é a velocidade.

Para minha aplicação específica, a equação pode ser escrita no formulário

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

Aqui está o esquema Crank-Nicolson que apliquei,

$\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 os termos $\alpha$ e $\beta$ . Isso permite que o esquema se mova entre:

• $\beta=\alpha=1/2$ Crank-Niscolson,
• $\beta=\alpha=1$ está totalmente implícito
• $\beta=\alpha=0$ é totalmente explícito

Os valores podem ser diferentes, o que permite que o termo de difusão seja Crank-Nicolson e o termo de advecção seja outra coisa. Qual é a abordagem mais estável, o que você recomendaria?

Esta é uma pergunta bem estruturada e uma coisa muito útil para entender. Korrok está correto em encaminhá-lo à análise de von Neumann e ao livro de LeVeque. Eu posso adicionar um pouco mais a isso. Gostaria de escrever uma resposta detalhada, mas no momento só tenho tempo para uma resposta curta:

Com , você obtém um método absolutamente estável para tamanhos de etapas arbitrariamente grandes, além de precisão de segunda ordem. No entanto, o método não é estável em L , portanto, frequências muito altas não serão amortecidas, o que não é físico.$\alpha=\beta=1/2$

Com , você obtém um método que também é incondicionalmente estável, mas com precisão de primeira ordem. Este método é muito dissipativo. É estável em L.$\alpha=\beta=1$

Se você usar , seu método pode ser entendido como a aplicação de um método aditivo de Runge-Kutta à semi-discretização de diferença centralizada. A análise de estabilidade e precisão de tais métodos é consideravelmente mais complicada. Um artigo muito bom sobre esses métodos está aqui .$\alpha\ne\beta$

Qual abordagem a recomendar depende fortemente da magnitude de , do tipo de dados iniciais com que você lida e da precisão que você procura. Se uma precisão muito baixa for aceitável, então é uma abordagem muito robusta. Se é moderado ou grande, o problema é dominado pela difusão e muito rígido; normalmente dará bons resultados. Se for muito pequeno, pode ser vantajoso usar um método explícito e um upwinding de ordem superior para os termos convectivos.α = β = 1 D α = β = 1 / 2 $D$$\alpha=\beta=1$$D$$\alpha=\beta=1/2$$D$

De um modo geral, convém usar um método implícito para equações parabólicas (a parte da difusão) - esquemas explícitos para PDE parabólico precisam ter um tempo muito curto para serem estáveis. Por outro lado, para a parte hiperbólica (advecção), você desejará um método explícito, pois é mais barato e não interrompe a simetria do sistema linear que você precisa resolver usando um esquema implícito de difusão. Nesse caso, você deseja evitar diferenças centralizadas como mudar para diferenças unilaterais por razões de estabilidade.( u j - u j - 1 ) / Δ $(u_{j+1}-u_{j-1})/2\Delta t$$(u_j-u_{j-1})/\Delta t$

Eu sugiro que você olhe o livro de Randy Leveque ou o livro de Dale Durran para "análise de estabilidade de von Neumann". É uma abordagem geral para determinar a estabilidade do seu esquema de discretização, desde que você tenha condições de contorno periódicas. (Há também um bom artigo wiki aqui .)

A idéia básica é supor que sua aproximação discreta possa ser escrita como uma soma das ondas planas , onde é o número da onda e a frequência. Você coloca uma onda plana na sua aproximação ao PDE e reza para que ela não exploda. Podemos reescrever a onda plana como e queremos garantir que . k ω ξ n e i k j Δ x | ξ | ≤ $e^{i(k j\Delta x-\omega n\Delta t)}$$k$$\omega$$\xi^n e^{ikj\Delta x}$$|\xi|\le 1$

A título de ilustração, considere a equação de difusão comum com diferenciação 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}$

Se substituirmos em uma onda plana, então dividimos por e , obtemos a equaçãoe i k j Δ $\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$

Limpe isso um pouco agora e teremos:

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

Isso é sempre menos que um, então você está limpo. Tente aplicar isso no esquema explícito e centrado da equação de advecção:

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

e ver o que você começa. (Desta vez, terá uma parte imaginária.) Você verá que , que são momentos tristes. Daí minha advertência de que você não a usa. Se você puder fazer isso, não deverá ter muita dificuldade em encontrar um esquema estável para a equação completa de difusão de advecção.$\xi$$|\xi|^2 > 1$

Dito isto, eu usaria um esquema totalmente implícito para a parte da difusão. Altere a diferenciação na parte de advecção para se e se e escolha um timestep para que . (Essa é a condição de Courant-Friedrichs-Lewy .) É apenas de primeira ordem, portanto, convém procurar esquemas de discretização de ordem superior, se isso lhe interessar. 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$