Randbedingungen für die Advektionsgleichung, diskretisiert durch eine Finite-Differenzen-Methode

14

Ich versuche einige Ressourcen zu finden, um zu erklären, wie man Randbedingungen wählt, wenn man Finite-Differenzen-Methoden zur Lösung von PDEs einsetzt.

Die Bücher und Notizen, auf die ich momentan Zugriff habe, sagen ähnliche Dinge aus:

Die allgemeinen Regeln für die Stabilität bei Vorhandensein von Grenzen sind für einen Einführungstext viel zu kompliziert. Sie erfordern ausgefeilte mathematische Maschinen

(A. Iserles Ein erster Kurs in der numerischen Analyse von Differentialgleichungen)

Beispiel: Wenn Sie versuchen, die 2-Schritt-Sprungmethode für die Advektionsgleichung zu implementieren:

uin+1=uin1+μ(ui+1nui1n)

mit 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

Die Lösung verhält sich gut, bis sie die Grenze erreicht, wenn sie sich sehr plötzlich schlecht verhält.

Wo kann ich lernen, mit solchen Randbedingungen umzugehen?

Simon Morris
quelle

Antworten:

12

Sloedes Antwort ist sehr gründlich und korrekt. Ich wollte nur ein paar Punkte hinzufügen, um das Erfassen zu erleichtern.

Grundsätzlich hat jede Wellengleichung eine inhärente Wellengeschwindigkeit und -richtung. Für eine eindimensionale Wellengleichung: die Wellengeschwindigkeit die Konstante a, die nicht nur die Geschwindigkeit bestimmt, mit der sich die Information in der Domäne ausbreitet, sondern auch ihre Richtung. Wenn a > 0 , gehen die Informationen von links nach rechts und wenn a < 0, ist es umgekehrt.

ut+einux=0
einein>0ein<0

Für die Schaltfroschmethode erhalten Sie, wenn Sie die Gleichungen diskretisieren: oder: u n i =u n - 2 i +μ(u n - 1 i + 1 -u n - 1 i - 1 )wobeiμ=-aΔt/Δx. In Ihrem Fall istμ>0

uichn-uichn-22Δt+einuich+1n-1-uich-1n-12Δx=0
uichn=uichn-2+μ(uich+1n-1-uich-1n-1)
μ=-einΔt/Δxμ>0was bedeutet, dass eine Welle nach links geht. Wenn Sie jetzt darüber nachdenken, benötigt eine Welle, die sich nach links bewegt, nur eine Randbedingung an der rechten Grenze, da alle Werte auf der linken Seite über ihre rechten Nachbarn aktualisiert werden. Die Angabe eines beliebigen Wertes an der linken Grenze entspricht nicht der Art des Problems. Bei bestimmten Methoden, wie dem einfachen Gegenwind, wird dies automatisch erledigt, da das Schema auch nur die richtigen Nachbarn in seiner Schablone enthält. Bei anderen Methoden wie dem Schaltfrosch müssen Sie einen "richtigen" Wert angeben.

Dies erfolgt normalerweise durch Extrapolation aus der internen Domäne, um den fehlenden Wert zu ermitteln. Bei mehrdimensionalen und nicht-kanonischen Problemen müssen dabei alle Eigenvektoren des Flusses Jacobi ermittelt werden, um zu bestimmen, welche Teile der Grenze tatsächlich Randbedingungen benötigen und welche Teile extrapoliert werden müssen.

GradGuy
quelle
Was würde es physikalisch bedeuten, diese Gleichung mit einer Randbedingung auf der linken und rechten Seite zu verwenden?
Frank
5

Allgemeine Antwort
Ihr Problem ist, dass Sie die Randbedingungen überhaupt nicht festlegen (oder sogar spezifizieren) - Ihr numerisches Problem ist schlecht definiert.

Grundsätzlich gibt es zwei Möglichkeiten, die Randbedingungen festzulegen:

  1. u0u101
  2. Ändern Sie die numerische Schablone so, dass nur Inneninformationen an der Grenze verwendet werden.

Welchen Weg Sie einschlagen, hängt stark von der Physik Ihres Problems ab. Für Wellengleichungsprobleme bestimmt man normalerweise die Eigenwerte des Flusses Jacobian, um zu entscheiden, ob äußere Randbedingungen erforderlich sind oder ob die innere Lösung verwendet werden soll (diese Methode wird allgemein als "Aufwind" bezeichnet).



uich-1nuich+1nichn+1ich=1u0nu100n+1u101n

u1nu100n

Unten finden Sie eine geänderte Version Ihres Quellcodes:

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
quelle
Schöne Antwort und willkommen bei scicomp, Sloede. Eine Frage, die ich normalerweise als "Aufwind" bezeichne, ist die Verwendung einer einseitigen Schablone, bei der Informationen nur von einer Grenze der Domäne bezogen werden. Wollten Sie das in Ihrer Antwort sagen?
Aron Ahmadia
1
Ja in der Tat. Entschuldigung, wenn meine Antwort nicht klar genug war. Im Allgemeinen bedeutet "Aufwind" jedoch nur, dass Sie die Richtung des Informationsflusses berücksichtigen. Dies muss nicht bedeuten, dass Sie eine Seite der Lösung vollständig verwerfen, sondern bedeutet lediglich, dass Sie dem Teil der Lösung den Vorzug geben, der in der Aufwindrichtung liegt.
Michael Schlottke-Lakemper
Wenn Sie N = 1000den Code etwas länger machen und ausführen, stellen Sie fest, dass er sich nicht so verhält, wie erwartet.
Simon Morris
Der Grund dafür ist, dass meine "Schnellkorrektur" -Lösung physikalisch nicht einwandfrei ist und darüber hinaus ziemlich empfindlich auf unechte Schwingungen in der Lösung reagiert. Verwenden Sie dies nicht für tatsächliche wissenschaftliche Berechnungen!
Michael Schlottke-Lakemper
2

Ich habe mir das etwas genauer angesehen und es scheint, dass dies (zumindest in den grundlegenden Fällen, die ich behandle) von der Gruppengeschwindigkeit der Methode abhängt.

Die leapfrog Methode ist (zum Beispiel):

uichn+1=uichn-1+μ(uich+1n-uich-1n)

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

e2ichωΔt=1+μeichωΔt(eichζΔx-e-ichζΔx)

Sünde(ωΔt)=μSünde(ζΔx)

dωdζ=cos(ζΔx)1-μ2sichn2(ζΔx)[-1,1]

Nun müssen wir die Gruppengeschwindigkeit der Randbedingungen herausfinden:

u1n+2=u1n+μu2n+1

Wir können die Grenzgruppengeschwindigkeit wie folgt berechnen:

2ichSünde(ωΔt)=μeichζΔx

Um also einige Gruppengeschwindigkeiten zu finden, die die Grenzen erlauben, müssen wir Folgendes finden:

ω=cζ

cos(ζΔx)=0,μSünde(ζΔx)=2Sünde(ζcΔt)

ζ=π2Δxμ=2Sünde(cμπ2)c[-1,1]μ


u0n+1=u1n[-1,1]

Ich muss noch ziemlich viel darüber nachlesen, bevor ich es vollständig verstehe. Ich denke, die Schlüsselwörter, die ich suche, sind GKS-Theorie.

Quelle für all diese A Iserles Part III Notizen


Eine klarere Berechnung dessen, was ich getan habe, finden Sie hier: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Simon Morris
quelle
-2

Leute, ich bin sehr neu auf dieser Seite. Vielleicht ist dies nicht der richtige Ort, aber bitte verzeihen Sie mir, da ich sehr neu hier bin :) Ich habe ein sehr ähnliches Problem, der einzige Unterschied ist die Startfunktion, die in meinem Fall eine Kosinuswelle ist. Mein Code lautet: Alles löschen; clc; alle schließen;

M = 1000; N = 2100;

mu = 0,5;

c = & mgr; 0 & ndash; u; f = @ (x) 1-cos (20 · pi · x-0,025). ^ 2; u = Nullen (M, N); x = 0: (1 / M): 0,05; u (1: Länge (x), 1) = f (x); u (1: Länge (x), 2) = f (x - mu / (M)); x = Linspace (0,1, M);

für i = 3: N warte;

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

Diagramm (x, u (:, i)); Achse ([0 1,5 - 0,5 2]) gezeichnet; % pause end

Es gibt bereits diesen Code hier, aber aus irgendeinem Grund, wahrscheinlich im Zusammenhang mit der Kosinuswelle, schlägt mein Code fehl: / Jede Hilfe wäre dankbar :) Danke!

John
quelle
2
Willkommen bei SciComp.SE! Sie sollten dies eine neue Frage machen. (Antworten sind nur für tatsächliche Antworten gedacht.) Wenn Sie den Link "Stellen Sie Ihre eigenen Fragen" unten verwenden (es ist dunkelgelb auf hellgelb, zugegebenermaßen ein bisschen schwer zu erkennen, wenn Sie nicht wissen, dass es da ist) wird die Frage automatisch mit dieser verknüpft. (Sie können auch einen Link zu dieser Frage hinzufügen.)
Christian Clason