Condições de contorno para a equação de advecção discretizada por um método de diferenças finitas

Estou tentando encontrar alguns recursos para ajudar a explicar como escolher condições de contorno ao usar métodos de diferenças finitas para resolver PDEs.

Os livros e notas aos quais tenho acesso atualmente dizem coisas semelhantes:

As regras gerais que governam a estabilidade na presença de limites são muito complicadas para um texto introdutório; eles exigem maquinaria matemática sofisticada

(A. Iserles, um primeiro curso de análise numérica de equações diferenciais)

Por exemplo, ao tentar implementar o método de salto em duas etapas para a equação de advecção:

$u_i^{n+1} = u_i^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

usando MATLAB

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);
    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

A solução se comporta bem até atingir o limite, quando de repente começa a se comportar mal.

Onde posso aprender a lidar com condições de contorno como essa?

A resposta de Sloede é muito completa e correta. Eu só queria acrescentar alguns pontos para facilitar a compreensão.

Basicamente, qualquer equação de onda tem uma velocidade e direção inerentes. Para uma equação de onda unidimensional: a velocidade da onda é a constante a que determina não apenas a velocidade na qual a informação está se propagando no domínio, mas também sua direção. Se um > 0 , as informações vão da esquerda para a direita e se um < 0 é o contrário.

$$u_t + a u_x = 0$$ $a$$a>0$$a<0$

Para o método do sapo com salto, quando você discretiza as equações, obtém: ou: u n i =u n - 2 i +μ(u n - 1 i + 1 -u n - 1 i - 1 ) emqueμ=-aΔt/Δx. No seu caso,μ>

$$\frac{u_i^{n} - u_i^{n-2}}{2 \Delta t} + a \frac{u_{i+1}^{n-1} - u_{i-1}^{n-1}}{2 \Delta x} = 0$$ $$u_i^{n} = u_i^{n-2} + \mu (u_{i+1}^{n-1} - u_{i-1}^{n-1})$$ $\mu = -a \Delta t/\Delta x$$\mu > 0$que se traduz em uma onda indo para a esquerda. Agora, se você pensar bem, uma onda que está viajando para a esquerda precisará apenas de uma condição de limite no limite direito, pois todos os valores à esquerda estão sendo atualizados por meio de seus vizinhos direitos. De fato, especificar qualquer valor no limite esquerdo é inconsistente com a natureza do problema. Em certos métodos, como simples contra o vento, isso é resolvido automaticamente, pois o esquema também envolve apenas os vizinhos corretos no estêncil. Em outros métodos, como o sapo de salto, você precisa especificar algum valor "correto".

Isso geralmente é feito por extrapolação do domínio interno para encontrar o valor que falta. Para problemas multidimensionais e não canônicos, isso envolve encontrar todos os vetores próprios do fluxo jacobiano para determinar quais partes da fronteira realmente precisam de condições de fronteira e quais partes requerem extrapolação.

Resposta geral
Seu problema é que você não define (nem sequer especifica) as condições de contorno - seu problema numérico está mal definido.

Geralmente, existem duas maneiras possíveis de especificar as condições de contorno:

1. $u_{0}$$u_{101}$
2. Altere o estêncil numérico para que ele use apenas informações internas no limite.

O caminho a seguir depende muito da física do seu problema. Para problemas do tipo equação de onda, geralmente determina-se os valores próprios do fluxo jacobiano para decidir se são necessárias condições de contorno externas ou se a solução interior deve ser usada (esse método é geralmente chamado de 'upwinding').

---

$u_{i-1}^n$$u_{i+1}^n$$i$$n+1$$i = 1$$u_{0}^n$$u_{100}^{n+1}$$u_{101}^n$

$u_{1}^n$$u_{100}^n$

Você pode encontrar uma versão modificada do seu código-fonte abaixo:

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    %u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);

    % Apply the numerical stencil to all interior points
    for j = 2:M-1
        u(j,i) = u(j,i-2) + mu*(u(j+1,i-1) - u(j-1,i-1));
    end

    % Set the boundary values by interpolating linearly from the interior
    u(1,i) = 2*u(2,i) - u(3,i);
    u(M,i) = 2*u(M-1,i) - u(M-2,i);

    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

Portanto, analisei isso com mais detalhes e parece que isso (pelo menos nos casos básicos com os quais estou lidando) depende da velocidade do grupo do método.

O método leapfrog (por exemplo) é:

$u_{i}^{n+1} = u_{i}^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

$u_k^n = e^{i( \zeta k \Delta x + \omega(\zeta) n \Delta t)}$

$e^{2 i \omega \Delta t} = 1 + \mu e^{i\omega \Delta t} ( e^{i \zeta \Delta x} - e^{-i \zeta \Delta x} )$

$\sin(\omega \Delta t) = \mu \sin (\zeta \Delta x)$

$\frac{d\omega}{d\zeta} = \frac{\cos (\zeta \Delta x)}{\sqrt{1-\mu^2 sin^2(\zeta \Delta x)}} \in [-1,1]$

Agora precisamos descobrir a velocidade do grupo das condições de contorno:

$u_1^{n+2} = u_1^n + \mu u_2^{n+1}$

Podemos calcular a velocidade do grupo de fronteira da seguinte maneira:

$2 i \sin(\omega \Delta t) = \mu e^{i \zeta \Delta x}$

Então, para encontrar algumas velocidades de grupo que os limites permitem, precisamos encontrar:

$\omega = c \zeta$

$\cos(\zeta \Delta x) = 0, \mu \sin(\zeta \Delta x) = 2 \sin (\zeta c \Delta t)$

$\zeta = \frac{\pi}{2 \Delta x}$$\mu = 2 \sin (\frac{c \mu \pi}{2})$$c \in [-1,1]$$\mu$

---

$u_0^{n+1} = u_1^{n}$$[-1,1]$

Ainda tenho muito a ler sobre isso antes de entender completamente. Penso que as palavras-chave que procuro são a teoria GKS.

Fonte para tudo isso notas de A Iserles Part III

---

Um cálculo mais claro do que eu fiz pode ser encontrado aqui: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Pessoal, sou muito novo neste site. Talvez este não seja o lugar para perguntar, mas, por favor, me perdoe, porque sou muito novo aqui :) Estou tendo um problema extremamente semelhante, a única diferença é a função inicial que, no meu caso, é uma onda cosseno. Meu código é este: limpe tudo; clc; feche tudo;

M = 1000; N = 2100;

mu = 0,5;

c = [mu 0-mu]; f = @ (x) 1- cos (20 * pi * x-0,025). ^ 2; u = zeros (M, N); x = 0: (1 / M): 0,05; u (1: comprimento (x), 1) = f (x); u (1: comprimento (x), 2) = f (x - mu / (M)); x = espaço interno (0,1, M);

para i = 3: N adie;

% Apply the numerical stencil to all interior points
for j = 2:M-1
    u(j,i) = u(j,i-2) - mu*(u(j+1,i-1) - u(j-1,i-1));
end

% Set the boundary values by interpolating linearly from the interior
u(M,i) =  2*u(M-1,i-1) - u(M-2,i-1);

plot (x, u (:, i)); desenho do eixo ([0 1,5 -0,5 2]); % final de pausa

Já existe esse código aqui, mas por algum motivo, provavelmente relacionado à onda cosseno, meu código falha: / qualquer ajuda seria apreciada :) obrigado!