Dziwna oscylacja podczas rozwiązywania równania doradczego metodą skończonej różnicy przy całkowicie zamkniętych warunkach brzegowych Neumanna (odbicie na granicach)

33

Próbuję rozwiązać równanie doradcze, ale w rozwiązaniu pojawia się dziwna oscylacja, gdy fala odbija się od granic. Jeśli ktokolwiek widział ten artefakt wcześniej, byłbym zainteresowany, aby poznać przyczynę i jak jej uniknąć!

To jest animowany gif, otwarty w osobnym oknie, aby wyświetlić animację (będzie odtwarzany tylko raz lub nie od razu zostanie buforowany!) Propagacja impulsu gaussowskiego

Zauważ, że propagacja wydaje się bardzo stabilna, dopóki fala nie zacznie odbijać się od pierwszej granicy. Jak myślisz, co może się tu dziać? Kilka dni spędziłem na sprawdzaniu kodu i nie mogę znaleźć żadnych błędów. Jest to dziwne, ponieważ wydaje się, że istnieją dwa rozwiązania propagujące: jedno pozytywne i jedno negatywne; po odbiciu od pierwszej granicy. Rozwiązania wydają się przemieszczać wzdłuż sąsiednich punktów siatki.

Szczegóły implementacji znajdują się poniżej.

Równanie doradcze,

ut=vux

gdzie to prędkość propagacji.v

Crank-Nicolson jest bezwarunkową (pdf link) stabilną dyskretyzacją dla równania doradczego, pod warunkiem, że powoli zmienia się w przestrzeni (zawiera tylko składowe niskich częstotliwości po transformacji Fouriera).u(x)

Zastosowałem dyskretyzację,

ϕjn+1ϕjnΔt=v[1β2Δx(ϕj+1nϕj1n)+β2Δx(ϕj+1n+1ϕj1n+1)]

Umieszczenie niewiadomych po prawej stronie umożliwia zapisanie tego w formie liniowej,

βrϕj1n+1+ϕjn+1βrϕj+1n+1=(1β)rϕj1n+ϕjn+(1β)rϕj+1n

gdzie (aby przyjąć średnią czasową równomiernie ważoną między punktem bieżącym a przyszłym) i .R = V Δ tβ=0.5r=vΔt2Δx

Te równania mają postać macierzy , gdzie,Aun+1=Mun

A=(1βr0βr1βrβr1βr0βr1)

M=(1(1β)r0(1β)r1(1β)r(1β)r1(1β)r0(1β)r1)

Wektory i są znanymi i nieznanymi wielkościami, dla których chcemy rozwiązać.U N + 1unun+1

Następnie stosuję zamknięte warunki brzegowe Neumanna na lewej i prawej granicy. Przez zamknięte granice rozumiem na obu interfejsach. W przypadku zamkniętych granic okazuje się, że (nie pokażę tutaj mojej pracy) musimy tylko rozwiązać powyższe równanie macierzowe. Jak wskazał @DavidKetcheson, powyższe równania macierzowe faktycznie opisują warunki brzegowe Dirichleta . W przypadku warunków brzegowych Neumannaux=0

A=(100βr1βrβr1βr001)

Aktualizacja

Zachowanie wydaje się dość niezależne od wyboru stałych, których używam, ale są to wartości wykresu, które widzisz powyżej:

  • v = 2
  • dx = 0,2
  • dt = 0,005
  • σ = 2 (gaussowski hwhm)
  • β = 0,5

Aktualizacja II

Symulacja z niezerowym współczynnikiem dyfuzji, (patrz komentarze poniżej), oscylacja zanika, ale fala już nie odbija !? Nie rozumiem dlaczego?D=1

Dyfuzja i porady

boyfarrell
źródło
co wziąłeś za ? v
Chris
v=2 w tych symulacjach. Zaktualizuję ustawienia symulacji. Dobry pomysł.
boyfarrell,
Wtedy oczekiwałbym, że warunek początkowy zostanie wskazany po prawej stronie i zniknie przez prawą granicę. Wszystko, co przychodzi na myśl, to to, że centralny schemat może dawać fałszywe oscylacje, chyba że zastosuje się go do równania doradzania-dyfuzji z liczbą komórek Pecleta poniżej 2. Być może wypróbować schemat wichury?
Chris
Czy uważasz, że w równaniu może wystąpić błąd znaku? W rzeczywistości moim ostatecznym celem jest zastosowanie tego równania rada-dyfuzja. Obecnie testuję różne przypadki ograniczające. W powyższym przykładzie współczynnik dyfuzji został ustawiony na zero. Dołączyłem nową animację powyżej. Nie rozumiem, dlaczego pik nie odzwierciedla się, gdy współczynnik dyfuzji jest niezerowy? Robi dokładnie tak, jak wspomniałeś (oprócz kierunku).
boyfarrell
Myślałem o , więc znak jest w porządku. Druga fabuła wygląda dla mnie dobrze. Dlaczego miałbyś oczekiwać, że coś się zastanowi? Może się to zdarzyć tylko wtedy, gdy zmiany jakoś się podpiszą. Spróbuj użyć schematu podmuchu dla porady zamiast schematu centralnego, wtedy powinieneś zobaczyć coś podobnego dla . v D = 0tu+vxu=0vD=0
Chris

Odpowiedzi:

28

Równanie, które rozwiązujesz, nie pozwala na właściwe rozwiązania, więc nie ma czegoś takiego jak odzwierciedlająca warunek brzegowy dla tego równania. Jeśli weźmiesz pod uwagę cechy, zdasz sobie sprawę, że możesz narzucić warunek brzegowy tylko na właściwej granicy. Próbujesz narzucić jednorodny warunek Dirichleta na lewej granicy, który jest matematycznie niepoprawny.

Dla przypomnienia: Metoda charakterystyki mówi, że roztwór ten musi być stała wzdłuż każdej linii, formy przez dowolną stałą . Zatem rozwiązanie wzdłuż lewej granicy jest określane przez rozwiązanie wcześniej w obrębie domeny problemu; nie możesz narzucić tam rozwiązania.xνt=CC

W przeciwieństwie do równania, Twój system liczbowy ma przyznać prawo rozwiązania toku. Prawidłowe tryby są nazywane trybami pasożytniczymi i obejmują bardzo wysokie częstotliwości. Zauważ, że fala skierowana w prawo to pakiet fali piłokształtnej, związany z najwyższymi częstotliwościami, które mogą być reprezentowane na twojej siatce. Ta fala jest czysto numerycznym artefaktem, stworzonym przez twoją dyskretyzację.

Dla podkreślenia: nie zapisałeś pełnego problemu wartości początkowej, który próbujesz rozwiązać. Jeśli to zrobisz, będzie jasne, że nie jest to dobrze postawiony problem matematyczny.

Cieszę się, że opublikowałeś to tutaj, ponieważ jest to piękna ilustracja tego, co może się zdarzyć, gdy dyskrecjonujesz problem, który nie jest dobrze postawiony, oraz zjawiska pasożytniczych modów. Duża +1 dla twojego pytania ode mnie.

David Ketcheson
źródło
dziękuję za dyskusję i poprawki. Nie doceniłem, że macierz nie będzie miała takich samych właściwości jak równanie różniczkowe. Więcej komentarzy do naśladowania ...
boyfarrell,
Tak, teraz widzę, jak to właściwie są warunki brzegowe Dirichleta, które zanotowałem powyżej dla korekty. To pierwszy raz, kiedy naprawdę próbuję zrozumieć proces rozwiązywania tych równań, ciągle uderzam w nierówności na drodze. Z przyjemnością opublikuję moje postępy!
boyfarrell
@David Ketcheson: Mam ten sam problem i opublikowałem swój problem w poniższym linku scicomp.stackexchange.com/questions/30329/... Czy możesz mi wyjaśnić, dlaczego według ciebie problem nie jest „dobrze postawiony matematycznie” ? Dzięki.
Herman Jaramillo
@HermanJaramillo Próbujesz narzucić wartość rozwiązania na lewej granicy, gdzie jest ona już określona przez PDE. Każdy podręcznik, który zawiera omówienie porady, będzie również wskazywał, jakie są prawidłowe warunki brzegowe i dlaczego. Dodałem drugi akapit z dodatkowym wyjaśnieniem; mam nadzieję, że to pomaga.
David Ketcheson
1
@HermanJaramillo: brak „matematycznie dobrze postawionego” oznacza w zasadzie, że masz dwa równania dla jednej wartości funkcji na granicy, warunek brzegowy, a także sam PDE. Zasadniczo te dwa równania są ze sobą sprzeczne. Mówiąc bardziej ogólnie, można to potraktować jako problem optymalizacji, w którym oba cele powinny być spełnione tak dobrze, jak to możliwe.
davidhigh
0

Wiele się nauczyłem z powyższych odpowiedzi. Chcę dołączyć tę odpowiedź, ponieważ uważam, że oferuje ona inne spojrzenie na problem.

Rozważmy równanie ze stałą prędkością fali .

uxx+1cutt=0.
c

Bez warunków początkowych i brzegowych to równanie ma rozwiązanie postaci . (puls przesuwa się w prawo)u(x,t)=f(xct)

Jeśli narzucimy warunek początkowy , wówczas rozwiązaniem równania w przedziale jest . Nie ma granicy i to jest rozwiązanie.u(x,t0)=p(x)x(,)p[xc(tt0)]

Załóżmy teraz, że definiujemy ograniczoną domenę , gdy wszystkie komputery mają ograniczoną pamięć. Następnie musimy określić wartości w punkcie i , w przeciwnym razie obliczeniowo utkniemy.x[a,b]ab

Jednym ze sposobów zdefiniowania warunków brzegowych jest użycie Dirichleta na lewej granicy i warunku spójnego z propagowanym rozwiązaniem. Oznacza to, że możemy zdefiniować (założyć czas początkowy ) t0

u(a,t0)=0,u(b,t0)=p[bc(tt0)].

Powoduje to pulsowanie w prawo, aż zniknie na prawej krawędzi.

naciśnij tutaj, aby uzyskać animację na Dirichlet na lewej granicy

Nadal słyszę hałas, którego nie rozumiem (ktoś mógłby tu pomóc, proszę?)

Inną opcją jest narzucenie okresowych warunków brzegowych. To znaczy zamiast narzucać warunek Dirichleta po lewej, możemy narzucić pakiet fal, który jest zgodny z granicą po lewej. To jest:

u(a,t0)=p[ac(tt0)],u(b,t0)=p[bc(tt0)].

Jednakże o i , i ponieważ trzeba umieścić dane wewnątrz przedziału dodamy jeden „okres” o długości a następnie znaleźć , aby warunek po lewej stronie wynosiłby (to samo!) i mieliśmy impuls wychodzący na w prawo i wchodząc w lewo.ac(tt0)<at>t0c>0[a,b]baac(tt0)+ba=a+b(tt0)u(a,t0)=p[bc(tt0]

Ten link pokazuje, co nazwałbym okresowymi warunkami brzegowymi.

Wykonałem animacje w pythonie, a schemat jest schematem Crank-Nicholson, jak wskazano w pytaniu tutaj.

Nadal walczę ze wzorem szumu po tym, jak fala uderza w pierwszą (prawą) granicę.

Herman Jaramillo
źródło
1
Nie widziałem animacji w moim telefonie komórkowym, ale uważam, że twój wzór szumu jest spowodowany brakiem dokładności liczbowej. Absorpcja działa tylko dlatego, że narzucasz dokładne rozwiązanie na granicę. Po dokładnym rozwiązaniu, rozwiązanie liczbowe osiągające prawą granicę nieznacznie różni się częstotliwością i fazą. To znów prowadzi do odbicia, a tym samym zakłóceń.
davidhigh
@davidhigh: Dziękujemy za informacje. Sprawdzę to. Przykro mi, animacja nie działała w twoim telefonie. To też nie działało w moim telefonie (Samsung). Może to być brakujące oprogramowanie w telefonach. Powinien działać na komputerze. Dzięki jeszcze raz.
Herman Jaramillo
Nie ma za co. Sama animacja nie jest tak ważna, chciałem tylko zobaczyć, jak duże są błędy. Przy okazji, narzucając dokładne rozwiązanie na granicę, unikasz „zła” z założenia. Oznacza to, że nadal masz dwa równania dla jednej wartości na granicy, ale stosując dokładny wynik, zmuszasz je do zachowania spójności. Ale działa to tylko w przybliżeniu, ponieważ wynik liczbowy nie jest całkowicie dokładny.
davidhigh
I jeszcze jeden komentarz: ogólnym analitycznym rozwiązaniem równania falowego jest superpozycja impulsu poruszającego się w lewo i jednego poruszającego się w prawo. W twoim przypadku rozważasz tylko puls poruszający się w prawo, więc zastosowałeś już warunki początkowe - w przeciwieństwie do tego, co podajesz w tekście.
davidhigh
@davidhigh: Myślałem trochę o twoich spostrzeżeniach na temat hałasu po tym, jak puls przekroczy granicę. Uważam, że masz rację i istnieje różnica między dokładnym rozwiązaniem analitycznym a propagowanym liczbowo impulsem. Na granicy ta niewielka różnica generuje widoczny tam mały szum. Przedstawiony w tej dyskusji system doradztwa CN jest dyspersyjny i uważam, że chociaż dyspersja nie jest zauważana, zanim impuls nie przekroczy granicy, może wyzwalać niewielkie zaburzenie (różnica między rozwiązaniami analitycznymi i numerycznymi) na granicy. Dzięki jeszcze raz.
Herman Jaramillo