Wybór wielkości kroku za pomocą ODE w Matlabie

Cześć i dzięki za poświęcenie czasu na spojrzenie na moje pytanie. To jest zaktualizowana wersja mojego pytania, które opublikowałem wcześniej w physics.stackexchange.com

Obecnie badam wirnik ekscytonowy 2D Kondensat Bosego-Einsteina i jestem ciekawy stanu podstawowego tego układu. Matematyczna metoda dojścia do stanu podstawowego nazywa się metodą wyobrażoną czasu .

Metoda jest bardzo prosta, w której czas w mechanice kwantowej zastępuje się urojonym jeden Podstawienie to powoduje, że cząsteczki o wysokiej energii w moim układzie rozpadają się szybciej niż cząstki o niskiej energii. Ponownie normalizując liczbę cząstek na każdym etapie obliczeń, otrzymujemy układ cząstek o najniższej energii, aka. stan podstawowy.

t = - i τ

$t = -i \tau$

Omawiane równanie jest nieliniowe, nazywane nieliniowym równaniem Schrödingera , czasem równaniem Grossa-Pitaevskiego . Aby rozwiązać problem, używam Matlabs ode45, który ewoluuje system do przodu i ostatecznie osiąga stan podstawowy.

Uwaga! Nieliniowe równanie Schrödingera obejmuje laplacian i niektóre inne wyrażenia różniczkowe w przestrzeni. Wszystko to rozwiązuje się za pomocą szybkiej transformacji Fouriera. W końcu mamy tylko czas ODE. *

Mój problem i pytanie: obliczenia idą od do . Ode45 jest umieszczany w pętli for , aby nie obliczał jednocześnie gigantycznego wektora . Pierwsza runda zaczynałaby się od od45 (odefun, ), a następnie następna od $t_0$ $t_f$ $[t_0,\dots,t_f]$ $[t_0, t_0+\Delta/2, t_0 + \Delta],y,\dots$ $t_0 + \Delta$ . Tutaj krok czasu jest moim problemem. Różne wybory kroków czasowych dają mi różne rozwiązania stanu podstawowego i nie mam pojęcia, jak ustalić, który krok czasowy daje mi „najbardziej” prawidłowy stan podstawowy! $\Delta$

Moja próba: zdaję sobie sprawę, że w tym schemacie duże kroki czasowe spowodują rozpad dużej liczby cząstek przed ponowną normalizacją do pierwotnej liczby cząstek, podczas gdy małe kroki czasowe spowodują rozpad mniejszej ilości cząstek przed ponowną normalizacją. Początkowo sądzę, że małe kroki czasowe powinny dać dokładniejsze rozwiązanie, ale wydaje się, że jest odwrotnie.

Nie jestem ekspertem numerycznym, więc wybór ode45 był po prostu arbitralny. ode113 daje mi to samo. :(

Czy ktoś ma jakieś przemyślenia na ten temat. Daj mi znać, jeśli potrzebne będą dodatkowe szczegóły.

Dziękuję Ci.

Aktualizacja 1: Badam urojoną metodę czasu i ODE. Wydaje się, że jeśli krok czasowy nie jest wystarczająco mały, wszystko staje się niestabilne. To sprawia, że zastanawiam się, czy moje równania nieliniowe są sztywne, co znacznie utrudnia to, co rozumiem. Będę Cie informować na bieżąco.

Aktualizacja 2: NAPRAWIONO: Problemem była normalizacja poza ODE. Jeśli normalizacja jest utrzymywana w odefun, wówczas ODE zwraca ten sam wynik dla różnych wyborów „zewnętrznych” kroków czasowych. Mój kolega pokazał mi starsze kody i po prostu dodałem jedną linię do mojego odefun.

function y_out = odefun(t,y_in,...variables...) 

    ...
    [ Nonlinear equations evaluated ]  
    ...


    y_out = y_out + 0.1*y_in*(N0-Ntemp) ;
end

Ostatni wiersz oblicza różnicę w bieżącej liczbie cząstek (Ntemp) i liczbie cząstek, które system powinien zatrzymać (N0). Dodaje część cząstek z powrotem do wyjścia, a tym samym tworzy całkowitą stabilność liczby cząstek w układzie zamiast ich wszystkie zanikają.

Zadam również nowe pytanie dotyczące wymiarów problemu i pewnych różnic w pracy z pikosekundami lub nanosekundami jako kroków czasowych w ODE.

Dziękuję wam wszystkim. :)

Podstawowym problemem jest to, że siłą używasz metody adaptacyjnej, takiej jak ode45()robienie równych kroków. Dlaczego właśnie unikasz generowania „gigantycznego wektora”? Jeśli absolutnie potrzebujesz równych punktów, ode45()postępuj jak zwykle, a następnie użyj interpolacji.

y

$y$

Jeśli pamięć służy, powinna istnieć opcja ode45()umożliwiająca zachowanie kroków większych niż określony próg; możesz przyjrzeć się temu.

Odpowiedzią jest użycie lokalnego oszacowania błędu. Jest jeden wbudowany w ODE45, więc najłatwiej jest z niego skorzystać, ale możesz też samodzielnie napisać kod.

David Ketcheson

0.1

$0.1$

1 / time

$1/\text{time}$

\frac{α}{Δ t} \cdot (N_{t} - N_{0})

$\frac{\alpha}{\Delta t} \cdot (N_t - N_0)$

Δ t

$\Delta t$

Odpowiedzi:

Ponieważ nie opublikowałeś kodu MATLAB, nie jestem pewien, jak dzwonisz ode45. Zgaduję, że zmieniasz wektor tspan (drugi argument) przy każdym wywołaniu ode45. Pierwszą rzeczą do zrozumienia jest to, że wektor tspan nie ma (prawie) żadnego wpływu na krok czasowy wykorzystywany przez ode45. Wektor tspan pozwala po prostu przekazać do ode45 przedział czasu integracji i o której godzinie chcesz uzyskać wynik. Krok czasu wykorzystywany przez algorytm Runga-Kutta w ode45 jest dostosowywany wewnętrznie, aby osiągnąć zalecaną dokładność. Dwa parametry kontrolujące tę dokładność to RelTol i AbsTol w strukturze opcji przekazanych do ode45. Mają rozsądne wartości domyślne, a ponieważ o nich nie wspomniałeś, zakładam, że ich nie zmieniłeś.

Powiedziałem „prawie” brak wpływu na normalny czas ode45. Jeśli żądasz danych wyjściowych w odstępach czasu bardzo małych w stosunku do kroku czasowego, który w przeciwnym razie zajmowałby ode45, wówczas konieczne będzie zmniejszenie kroku czasowego, aby spełnić twoje żądanie wyjściowe. Wierzę, że tak właśnie zakłada JM. Dlaczego potrzebujesz rozwiązania przy tak wielu wyjściach? Zwykle wystarczy po prostu zażądać danych wyjściowych w wystarczającym czasie, aby wygenerować gładki wykres.

Jeśli chodzi o zmianę w rozwiązaniu, które widzisz, być może domyślne wartości RelTol i AbsTol nie są odpowiednie dla twojego problemu. Sugeruję zastąpienie pętli na ode45 pojedynczym wywołaniem, żądanie wyjścia w rozsądnej liczbie razy i eksperymentowanie z mniejszymi wartościami RelTol i AbsTol, aż do uzyskania zbieżnego rozwiązania.

Bill Greene
źródło

Dziękuję za odpowiedź. Powodem, dla którego potrzebowałem rozwiązania przy tak wielu wyjściowych czasach, jest to, że jeśli funkcja fali nie jest normalizowana regularnie, wszystko się rozkłada i mój system jest pusty. Dlatego umieściłem ode45 w pętli z małymi wektorami tspan, aby móc ponownie normalizować po każdym wektorze tspan.

Ponieważ nieliniowe równanie Schrödingera jest, cóż, nieliniowe, może mieć dużą liczbę stanów stacjonarnych, z których niektóre mogą być stabilne. W rzeczywistości fizycznej, poczynając od jednego określonego stanu, system deterministycznie ewoluuje w jeden stan końcowy. Jeśli schemat numeryczny daje różne wyniki dla różnych dyskretyzacji (przedziałów czasowych), to jest to podstawowa wada dyskretyzacji. Sprawdź swój kod.

$\psi_0$

\frac{d ψ}{d t} = F (ψ),

$\frac{\text{d}\psi}{\text{d}t} = F(\psi),$

F (ψ_{0}) = 0.

$F(\psi_0) = 0.$

G (ψ) = \int_{Ω} E (ψ)

$G(\psi) = \int_\Omega E(\psi)$

E (\cdot)

$E(\cdot)$

F (ψ) = 0

$F(\psi)=0$

E (ψ)

$E(\psi)$

E (ψ) = - | ψ |^{4}

$E(\psi) = -|\psi|^4$

Nico Schlömer
źródło

Tak. Wykreślam profil gęstości mojego rozwiązania wyjściowego i kiedy nie zmienia się on przez długi czas, zasadniczo przestał ewoluować, zakładam, że osiągnąłem stan stacjonarny. Ale nie jestem pewien, czy spojrzenie na gęstość energii może pomóc, ponieważ funkcją fali jest spinor z (+2, +1, -1, -2) składnikami spinowymi. Nie sądzę, że zintegrowanie każdego elementu powie mi energię kondensatu, ale kiedy dojdę do stanu podstawowego, gęstość energii powinna być stacjonarna, a zatem stała w czasie, jest to wskazówka dla prawidłowego rozwiązania.

Problem rozwiązany:

Normalizacja musi być częścią funkcji ocenianej w ODE. Rozbijanie ODE na wiele etapów i normalizowanie między nimi powoduje pozornie niestabilność numeryczną i daje różne wyniki w zależności od przedziałów czasowych, w których ODE jest dzielony. (Aby uzyskać więcej informacji, patrz edycja 2 w pytaniu).

źródło