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

14

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:

ujan+1=ujan-1+μ(uja+1n-uja-1n)

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?

Simon Morris
źródło

Odpowiedzi:

12

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.

ut+zaux=0
zaza>0za<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μ>0

ujan-ujan-2)2)Δt+zauja+1n-1-uja-1n-12)Δx=0
ujan=ujan-2)+μ(uja+1n-1-uja-1n-1)
μ=-zaΔt/Δxμ>0co 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.

GradGuy
źródło
Fizycznie, co oznaczałoby użycie tego równania z warunkiem brzegowym po lewej i prawej stronie?
Frank
5

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. u0u101
  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ę”).



uja-1nuja+1njan+1ja=1u0nu100n+1u101n

u1nu100n

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
Michael Schlottke-Lakemper
źródło
Ładna odpowiedź i witamy w scicomp, Sloede. Jedno pytanie, zwykle widzę „przewijanie w górę” zdefiniowane jako użycie jednostronnego szablonu, w którym informacje są pobierane tylko z jednej granicy domeny. Czy chciałeś powiedzieć to w swojej odpowiedzi?
Aron Ahmadia,
1
W rzeczy samej. Przepraszam, jeśli moja odpowiedź nie była wystarczająco jasna. Zasadniczo jednak „przewrócenie w górę” oznacza po prostu uwzględnienie kierunku przepływu informacji. Nie musi to oznaczać, że całkowicie odrzucasz jedną stronę rozwiązania, oznacza to tylko, że preferujesz tę część rozwiązania, która leży w kierunku „pod wiatr”.
Michael Schlottke-Lakemper
Jeśli utworzysz N = 1000i uruchomisz kod nieco dłużej, okaże się, że nie działa on zgodnie z oczekiwaniami.
Simon Morris
Powodem tego jest to, że moje rozwiązanie „szybkiej naprawy” nie jest fizycznie poprawne, a ponadto dość wrażliwe na fałszywe oscylacje w rozwiązaniu. Nie używaj tego do rzeczywistych obliczeń naukowych!
Michael Schlottke-Lakemper
2

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:

ujan+1=ujan-1+μ(uja+1n-uja-1n)

ukn=mija(ζkΔx+ω(ζ)nΔt)

mi2)jaωΔt=1+μmijaωΔt(mijaζΔx-mi-jaζΔx)

grzech(ωΔt)=μgrzech(ζΔx)

reωreζ=sałata(ζΔx)1-μ2)sjan2)(ζΔx)[-1,1]

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

u1n+2)=u1n+μu2)n+1

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

2)jagrzech(ωΔt)=μmijaζΔx

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

ω=doζ

sałata(ζΔx)=0,μgrzech(ζΔx)=2)grzech(ζdoΔt)

ζ=π2)Δx dałbym μ=2)grzech(doμπ2)) dla którego rozwiązanie do[-1,1]będzie istnieć. (Dla większości wyborówμ przynajmniej)


Rozwiązaniem, które znalazłem w literaturze, jest podjęcie u0n+1=u1n 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

Simon Morris
źródło
-2

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!

Jan
źródło
2
Witamy w SciComp.SE! Powinieneś uczynić to nowym pytaniem. (Odpowiedzi są przeznaczone tylko na rzeczywiste odpowiedzi.) Jeśli użyjesz linku „zadaj własne pytanie” u dołu (jest ciemnożółty na jasnożółtym, co prawda trochę trudny do zauważenia, jeśli nie wiesz, że tam jest) , automatycznie połączy pytanie z tym pytaniem. (Możesz również zamieścić link do tego pytania w swoim.)
Christian Clason,