Como discretizar a equação de advecção usando o método Crank-Nicolson?

8

A equação de advecção precisa ser discretizada para ser usada no método Crank-Nicolson. Alguém pode me mostrar como fazer isso?

pandoragami
fonte
1
Bem-vindo ao SciComp! O escopo da sua pergunta se encaixa muito bem neste site. Para obter boas respostas, no entanto, você deve ser mais específico. Indique o que em particular você não entende. Seu código parece bem estruturado e documentado, mas, para responder sua pergunta, talvez seja necessário um trecho de código menor.
Jan
Isso pode tornar sua depuração mais simples se você usar um caso de teste em que seu vetor de entrada tenha, digamos, 5 elementos, e percorrer o código com um depurador como gdb e ddd. Fazer isso pode ajudar a restringir a fonte do erro. Eu acho que a maioria das perguntas sobre depuração de código não funciona muito bem aqui, porque a maioria do trabalho envolve descobrir onde o bug está em primeiro lugar. Depois de encontrá-lo, a explicação é freqüentemente (mas nem sempre) direta. Você pode executar um teste de unidade para descobrir se há algum erro no Tridiagonal que possa estar causando esse comportamento?
Geoff Oxberry
Veja este exemplo na entrada da Wikipedia sobre o método Crank-Nicolson. Se você definir e como zero, isso se tornará um problema de advecção simples. Resta incorporar as condições de contorno ...D xkDx
Jan
Desculpe pessoal, mas o método Crank-Nicolson é totalmente inadequado para um problema de advecção. Infelizmente, a estabilidade e a precisão da aproximação diferencial local não garantem consistência. Você obterá a solução de um problema diferente. Em alguns casos triviais, você pode ter sorte, mas a equação de advecção é geralmente implacável. Verifique os livros de análise numérica. O esquema Crank-Nicolson nunca é usado. As pessoas ainda têm de adotá-lo em alguns problemas de engenharia, por causa da estabilidade, mas eles têm de controlar a solução à medida que evolui e tem heurística para Erro de teste

Respostas:

19

Começando com a equação de advecção é de forma conservadora,

ut=(vu)x+s(x,t)

O método Crank-Nicolson consiste em uma diferença centrada no tempo médio.

ujn+1ujnΔt=v[1β2Δx(uj+1nuj1n)+β2Δx(uj+1n+1uj1n+1)]+s(x,t)

Em relação à notação, os subscritos são para pontos no espaço e os sobrescritos são para pontos no tempo.

Os pontos em estão no futuro: são desconhecidos. Agora precisamos reorganizar a equação acima para que todos os conhecidos estejam no rhs e os desconhecidos no lhs.n+1

Fazendo a substituição,

r=v2ΔtΔx

dá,

βrϕj1n+1+ϕjn+1+βrϕj+1n+1=(1β)rϕj1n+ϕjn(1β)rϕj+1n

Esta é a equação de advecção discretizada usando o método Crank-Nicolson. Você pode escrevê-lo como uma equação matricial,

β=1/2

(1βr0βr1βrβr1βr0βr1)(u1n+1u2n+1uJ1n+1uJn+1)=(1(1β)r0(1β)r1(1β)r(1β)r1(1β)r0(1β)r1)(u1nu2nuJ1nuJn)
Definir fornecerá uma integração trapezoidal no tempo, portanto, para Crank-Nicolson, é isso que você deseja.β=1/2

Algumas palavras de aviso. Esta é a solução básica que você queria, mas precisará incluir algum tipo de condição de contorno para um problema bem colocado. Além disso, Crank-Nicolson não é necessariamente o melhor método para a equação de advecção. É de segunda ordem precisa e incondicionalmente estável , o que é fantástico. No entanto, ele gerará (como em todos os estênceis de diferença centralizada) uma oscilação espúria se você tiver soluções pontiagudas muito acentuadas ou condições iniciais.

Eu escrevi o código a seguir para você em Python, ele deve começar. O código resolve a equação de advecção para uma curva gaussiana inicial movendo-se para a direita com velocidade constante.

Curva Gaussiana se movendo para a direita com velocidade constante

from __future__ import division
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
import numpy as np
import pylab

def make_advection_matrices(z, r):
    """Return matrices A and M for advection equations"""
    ones = np.ones(len(z))
    A = spdiags( [-beta*r, ones, beta*r], (-1,0,1), len(z), len(z) )
    M = spdiags( [(1-beta) * r, ones, -(1-beta) * r], (-1,0,1), len(z), len(z) )
    return A.tocsr(), M.tocsr()

def plot_iteration(z, u, iteration):
    """Plot the solver progress"""
    pylab.plot(z, u, label="Iteration %d" % iteration)

# Set up basic constants
beta = 0.5
J = 200 # total number of mesh points
z = np.linspace(-10,10,J) # vertices
dz = abs(z[1]-z[0]) # space step
dt = 0.2    # time step
v = 2 * np.ones(len(z)) # velocity field (constant)
r = v / 2 * dt / dz

# Initial conditions (peak function)
gaussian = lambda z, height, position, hwhm: height * np.exp(-np.log(2) * ((z - position)/hwhm)**2)
u_init = gaussian(z, 1, -3, 2)

A, M = make_advection_matrices(z, r)
u = u_init
for i in range(10):
    u = spsolve(A, M * u)
    plot_iteration(z, u, i)

pylab.legend()
pylab.show()
boyfarrell
fonte
3
Na verdade, também vi sua pergunta anterior, mas era tão geral que não consegui responder (quando você postou uma página de código). Na minha experiência, quando você faz uma boa pergunta neste site, as pessoas são muito úteis. Boa Viagem!
boyfarrell
Estou apenas brincando.
Pandoragami 28/05
@boyfarrel Existe alguma chance de você ter uma versão C ++ / C disso? Tudo bem se não. Não uso muito o Matlab e não sinto vontade de aprender. Até o fortran seria melhor.
Pandoragami 28/05
2
<comentários inapropriados excluídos>
Paul