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

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!)

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,

$\frac{\partial u}{\partial t} = \boldsymbol{v}\frac{\partial u}{\partial x}$

gdzie to prędkość propagacji.$\boldsymbol{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ę,

$\frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right]$

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

$\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} -\beta r\phi_{j+1}^{n+1} = -(1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} + (1-\beta)r\phi_{j+1}^{n}$

gdzie (aby przyjąć średnią czasową równomiernie ważoną między punktem bieżącym a przyszłym) i .R = V Δ $\beta=0.5$$r=\boldsymbol{v}\frac{\Delta t}{2\Delta x}$

Te równania mają postać macierzy , gdzie,$A\cdot u^{n+1} = M\cdot u^n$

$\boldsymbol{A} = \left( \begin{matrix} 1 & -\beta r & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & \beta r & 1 \\ \end{matrix} \right)$

$\boldsymbol{M} = \left( \begin{matrix} 1 & (1 - \beta)r & & & 0 \\ -(1 - \beta)r & 1 & (1 - \beta)r & & \\ & \ddots & \ddots & \ddots & \\ & & -(1 - \beta)r & 1 & (1 - \beta)r \\ 0 & & &-(1 - \beta)r & 1 \\ \end{matrix} \right)$

Wektory i są znanymi i nieznanymi wielkościami, dla których chcemy rozwiązać.U N + $u^n$$u^{n+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 Neumanna$\frac{\partial u}{\partial x} = 0$

$\boldsymbol{A} = \left( \begin{matrix} 1 & 0 & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & 0 & 1 \\ \end{matrix} \right)$

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:

• $\boldsymbol{v}$ = 2
• dx = 0,2
• dt = 0,005
• $\sigma$ = 2 (gaussowski hwhm)
• $\beta$ = 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$

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-\nu t = C$$C$

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.

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 .

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

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,t_0)=p(x)$$x \in (-\infty, \infty)$$p[x-c(t-t_0)]$

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 \in [a,b]$$a$$b$

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

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:

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.$a-c(t-t_0) <a$$t>t_0$$c>0$$[a,b]$$b-a$$a-c(t-t_0)+b-a=a+b(t-t_0)$$u(a,t_0)=p[b-c(t-t_0]$

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ę.