Warunki brzegowe dla równania doradczego dyskretyzowanego metodą różnic skończonych

Próbuję znaleźć zasoby, które pomogą wyjaśnić, jak wybrać warunki brzegowe podczas korzystania z metod różnic skończonych do rozwiązywania PDE.

Książki i notatki, do których mam obecnie dostęp, mówią podobne rzeczy:

Ogólne zasady rządzące stabilnością w obecności granic są zdecydowanie zbyt skomplikowane, aby można było wprowadzić tekst wprowadzający; wymagają skomplikowanych mechanizmów matematycznych

(A. Iserles Pierwszy kurs numerycznej analizy równań różniczkowych)

Na przykład, próbując zaimplementować 2-etapową metodę skoku żab dla równania porady:

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

za pomocą 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

Rozwiązanie zachowuje się ładnie, dopóki nie osiągnie granicy, gdy nagle nagle zaczyna źle się zachowywać.

Gdzie mogę nauczyć się radzić sobie z takimi warunkami brzegowymi?

Odpowiedź Sloede jest bardzo dokładna i poprawna. Chciałem tylko dodać kilka punktów, aby ułatwić zrozumienie.

Zasadniczo każde równanie fali ma nieodłączną prędkość i kierunek fali. Dla jednowymiarowego równania falowego: prędkość fali jest stałą a, która określa nie tylko prędkość, z jaką informacja rozprzestrzenia się w dziedzinie, ale także jej kierunek. Jeśli jest > 0 , informacja idzie od lewej do prawej, a jeśli jest < 0 to na odwrót.

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

W przypadku metody skoku żaby, gdy dyskretyzujesz równania, otrzymujesz: lub: u n i =u n - 2 i +μ(u n - 1 i + 1 -u n - 1 i - 1 )gdzieμ=-aΔt/Δx. W twoim przypadkuμ>

$$\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$co przekłada się na falę skierowaną w lewo. Teraz, jeśli się nad tym zastanowić, fala, która płynie w lewo, będzie potrzebować jedynie warunku granicznego na prawej granicy, ponieważ wszystkie wartości po lewej stronie są aktualizowane przez ich prawych sąsiadów. W rzeczywistości określenie dowolnej wartości na lewej granicy jest niezgodne z charakterem problemu. W niektórych metodach, takich jak prosty pod wiatr, jest to załatwiane automatycznie, ponieważ schemat obejmuje również tylko właściwych sąsiadów. W innych metodach, takich jak żaba skokowa, musisz podać pewną „poprawną” wartość.

Zwykle odbywa się to poprzez ekstrapolację z wewnętrznej domeny, aby znaleźć brakującą wartość. W przypadku problemów wielowymiarowych i niekanonicznych wymaga to znalezienia wszystkich wektorów własnych strumienia Jakuba, aby ustalić, które części granicy faktycznie wymagają warunków brzegowych, a które części wymagają ekstrapolacji.

Ogólna odpowiedź
Twój problem polega na tym, że wcale nie ustawiasz (ani nawet nie określasz) warunków brzegowych - twój problem liczbowy jest źle zdefiniowany.

Zasadniczo istnieją dwa możliwe sposoby określenia warunków brzegowych:

1. $u_{0}$$u_{101}$
2. Zmień szablon numeryczny, aby używał tylko informacji wewnętrznych na granicy.

To, w którą stronę idziesz, zależy w dużym stopniu od fizyki twojego problemu. W przypadku problemów typu równania falowego zwykle określa się wartości własne strumienia Jakobianu, aby zdecydować, czy potrzebne są zewnętrzne warunki brzegowe, czy też należy zastosować rozwiązanie wewnętrzne (ta metoda jest powszechnie nazywana „przewijaniem w górę”).

---

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

Poniżej znajdziesz zmodyfikowaną wersję swojego kodu źródłowego:

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

Spojrzałem na to bardziej szczegółowo i wydaje się, że to (przynajmniej w podstawowych przypadkach, którymi się zajmuję) zależy od prędkości grupowej metody.

Metoda leapfrog (na przykład) to:

$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]$

Teraz musimy dowiedzieć się, jaka jest prędkość grupowa warunków brzegowych:

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

Możemy obliczyć prędkość grupy granicznej w następujący sposób:

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

więc aby znaleźć pewne prędkości grupowe, na które pozwalają granice, musimy znaleźć:

$\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}$ dałbym $\mu = 2 \sin (\frac{c \mu \pi}{2})$ dla którego rozwiązanie $c \in [-1,1]$będzie istnieć. (Dla większości wyborów$\mu$ przynajmniej)

---

Rozwiązaniem, które znalazłem w literaturze, jest podjęcie $u_0^{n+1} = u_1^{n}$ ponieważ ma on liczbę fal granicznych, która leży na zewnątrz $[-1,1]$.

Wciąż muszę przeczytać o tym więcej, zanim całkowicie to zrozumiem. Myślę, że kluczowymi słowami, których szukam, są teoria GKS.

Źródło wszystkich notatek A Iserles część III

---

Dokładniejsze obliczenia tego, co zrobiłem, można znaleźć tutaj: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Chłopaki, jestem bardzo nowy na tej stronie. Może nie jest to miejsce, o które warto zapytać, ale proszę wybacz mi, ponieważ jestem tu bardzo nowy :) Mam bardzo podobny problem, jedyną różnicą jest funkcja początkowa, która w moim przypadku jest falą cosinusową. Mój kod jest następujący: wyczyść wszystko; clc; zamknij wszystko;

M = 1000; N = 2100;

mu = 0,5;

c = [mu 0-im]; f = @ (x) 1-cos (20 * pi * x-0,025). ^ 2; u = zera (M, N); x = 0: (1 / M): 0,05; u (1: długość (x), 1) = f (x); u (1: długość (x), 2) = f (x - mu / (M)); x = przestrzeń linowa (0,1, M);

dla i = 3: N wstrzymanie;

% 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);

wykres (x, u (:, i)); oś ([0 1,5 -0,5 2]) drawnow; % pauzy na końcu

Jest już ten kod, ale z jakiegoś powodu, prawdopodobnie związany z falą cosinusową, mój kod zawodzi: / każda pomoc byłaby mile widziana :) dzięki!