Como você depura o código numérico, qual poderia ser a fonte desse erro oscilatório?

Pode-se obter uma grande quantidade de insights sobre a experiência, apenas queria saber se alguém já viu algo semelhante a isso antes. O gráfico mostra a condição inicial (verde) para a equação de advecção-difusão, depois a solução na iteração 200 (azul) e depois novamente na iteração 400 (vermelha).

A solução da equação de advecção-difusão explode após algumas iterações. O número do e a condição da CFL são satisfeitos, , portanto as equações devem ser estáveis. Antecipo que tenho um erro no código numérico.$\mu\approx0.07$$C\approx 0.0015$

Fundo. A discretização é a diferença central para termos de advecção e difusão. Acredito que esta seja a primeira ordem de advecção e a segunda ordem de difusão. Eu implementei isso usando uma abordagem de volume finito (pela primeira vez) na qual os valores dos coeficientes (velocidade e coeficiente de difusão) nas faces das células são encontrados por interpolação linear a partir das médias das células. Aplico a condição de limite de Robin nas superfícies esquerda e direita e defino o fluxo nos limites como zero.

Como você depura seu código numérico? Alguém já viu algo parecido com isso antes, onde seria um bom lugar para começar a procurar?

Atualizar

• Aqui estão minhas notas pessoais de estilo "livro de laboratório" sobre a implementação de um método de volume finito para a equação de difusão de advecção, http://danieljfarrell.github.io/FVM/
• O código fonte do Python está disponível aqui, http://github.com/danieljfarrell/FVM.git

Atualizar

A solução não poderia ser mais simples! Acabei de cometer um erro de sinal no termo de difusão. É estranho, tenho certeza que não tinha postado isso, não teria encontrado o erro! Se alguém quiser compartilhar dicas sobre como depurar seu código numérico, ainda estou interessado. Eu não tenho um método, é um pouco imprevisível, continuo tentando coisas para obter pistas, mas esse processo pode levar semanas (às vezes).

Prova de que funciona ( NB: com o método de volume finito, tudo o que você precisa fazer para calcular a área é um somatório de largura altura para todas as células, se você usar um método de integração como numpy.trapz, seus resultados incluirão o código numérico erro do método trapézio). O que esta acontecendo aqui? Existem coeficientes de velocidade e difusão constantes, mas com condições de contorno fechadas . Portanto, na fronteira, vemos o equilíbrio entre o campo de velocidade empurrando para a direita e o empurrão de difusão para a esquerda.$\times$

Reuni algumas das minhas experiências em depuração de códigos numéricos aqui: deal.II FAQ: depuração . Não sei se isso teria ajudado você nesse caso em particular, mas pode ser em outros.