Co oznacza „symplektyczny” w odniesieniu do integratorów numerycznych i czy używa ich odeint SciPy?

W tym komentarzu napisałem:

... domyślny integrator SciPy, który, jak zakładam, używa tylko metod symplektycznych.

w którym mam na myśli SciPy odeint, który używa albo „metody niesztywnej (Adamsa)”, albo „metody sztywnej (BDF)”. Według źródła :

def odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0,
           ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0,
           hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
           mxords=5, printmessg=0):
    """
    Integrate a system of ordinary differential equations.

    Solve a system of ordinary differential equations using lsoda from the
    FORTRAN library odepack.

    Solves the initial value problem for stiff or non-stiff systems
    of first order ode-s::
        dy/dt = func(y, t0, ...)
    where y can be a vector.
    """

Oto przykład, w którym propaguję orbitę satelity wokół Ziemi przez trzy miesiące, aby pokazać, że zachodzi zgodnie z oczekiwaniami.

Uważam, że integratory niesymplektyczne mają niepożądaną właściwość polegającą na tym, że nie będą oszczędzać energii (ani innych ilości), a zatem są niepożądane na przykład w mechanice orbitalnej. Ale nie jestem do końca pewien, co sprawia, że integrator symplektyczny staje się symplektyczny.

Czy możliwe jest wyjaśnienie, czym jest własność (co czyni symplektycznego integratora symplektycznym) w prosty i (dość) łatwy do zrozumienia, ale nie niedokładny sposób? Pytam z punktu widzenia tego, jak integrator funkcjonuje wewnętrznie , a nie jak działa w testach.

I czy moje podejrzenie jest słuszne odeinti wykorzystuje wyłącznie integratory symplektyczne?

Zasadniczo, powinieneś mieć nadzieję, że integrator czarnej skrzynki jest symplektyczny, jeśli wymaga oddzielenia równań pozycji i pędu.

origimbo

@origimbo Thanks. Tak, i wygląda na to, że odeintjest to wrapper Pythona dla dość starych, sprawdzonych i dobrze przywoływanych kodów źródłowych, (zredagowane pytanie, referencje ODEPACK i LSODA), chociaż z pewnością przyznaję, że używam go w trybie czarnej skrzynki. Mój powiązany przykład pokazuje, że wektor stanu 6D składa się z trzech pozycji i trzech prędkości.

uhoh

Integratory ODE w ODEPACK i LSODA nie są integratorami symplektycznymi.

Brian Borchers

Oto sprawdzony przykład porównujący dwa bardzo proste solwery: Euler i Symplectic Euler: idontgetoutmuch.wordpress.com/2013/08/06/… .

idontgetoutmuch

Książka Hairera, Nørsetta i Wannera daje dobre wyjaśnienie metod symplektycznych. Spójrz w szczególności na rysunek 16.1 i liczby tutaj .

Odpowiedzi:

Zacznę od korekt. Nie, odeintnie ma żadnych integratorów symplektycznych. Nie, integracja symplektyczna nie oznacza zachowania energii.

Co oznacza symplectic i kiedy należy go używać?

Przede wszystkim, co oznacza symplectic? Symplectic oznacza, że rozwiązanie istnieje na rozmaitości symplektycznej. Kolektor symplektyczny to zestaw rozwiązań, który jest zdefiniowany przez 2-formę. Szczegóły rozmaitości symplektycznych prawdopodobnie brzmią jak matematyczny nonsens, więc zamiast tego ich istotą jest bezpośrednia zależność między dwoma zestawami zmiennych na takim różnorodności. Powodem, dla którego jest to ważne dla fizyki, jest fakt, że równania Hamiltona mają naturalnie, że rozwiązania znajdują się na rozmaitości symplektycznej w przestrzeni fazowej, przy czym naturalny podział jest składnikiem pozycji i pędu. Dla prawdziwego rozwiązania hamiltonowskiego ta ścieżka przestrzeni fazowej jest stałą energią.

Integrator symplektyczny to integrator, którego rozwiązanie opiera się na rozmaitości symplektycznej. Z powodu błędu dyskretyzacji, gdy rozwiązuje układ hamiltonowski, nie otrzymuje dokładnie prawidłowej trajektorii na kolektorze. Zamiast trajektorii, która sama w sobie jest zaburzony w celu od prawdziwej trajektorii. Następnie występuje dryf liniowy z powodu błędu numerycznego tej trajektorii w czasie. Normalni integratorzy mają tendencję do dryfowania kwadratowego (lub więcej) i nie mają żadnych dobrych globalnych gwarancji dotyczących tej ścieżki przestrzeni fazowej (tylko lokalna). $\mathcal{O}(\Delta t^n)$ $n$

Oznacza to, że integratory symplektyczne mają tendencję do wychwytywania wzorców długookresowych lepiej niż normalne integratory z powodu tego braku dryfu i prawie gwarancji okresowości. Ten notatnik wyświetla te właściwości dobrze w przypadku problemu Keplera . Pierwsze zdjęcie pokazuje, o czym mówię, z okresową naturą rozwiązania.

Zostało to rozwiązane za pomocą integratora symplektycznego 6. rzędu od Kahan i Li z DifferentialEquations.jl . Widać, że energia nie jest dokładnie zachowana, ale jej zmienność zależy od tego, jak daleko zaburzony kolektor rozwiązania znajduje się od prawdziwego kolektora. Ale ponieważ samo rozwiązanie numeryczne opiera się na rozmaitości symplektycznej, wydaje się być prawie dokładnie okresowe (z pewnym liniowym przesunięciem liczbowym, który można zobaczyć), co czyni go bardzo ładnym dla długoterminowej integracji. Jeśli zrobisz to samo z RK4, możesz dostać katastrofę:

Widać, że problem polega na tym, że w rozwiązaniu numerycznym nie ma prawdziwej okresowości, dlatego nadgodziny mają tendencję do dryfowania.

Podkreśla to prawdziwy powód, aby wybrać integratory symplektyczne: integratory symplektyczne dobrze sprawdzają się w długoterminowych integracjach z problemami, które mają własności symplektyczne (układy Hamiltona) . Przejdźmy zatem przez kilka rzeczy. Zauważ, że nie zawsze potrzebujesz integratorów symplektycznych nawet w przypadku problemu symplektycznego. W tym przypadku dobrze sprawdza się adaptacyjna metoda Runge-Kutta 5. rzędu. Oto Tsit5:

Zwróć uwagę na dwie rzeczy. Po pierwsze, uzyskuje wystarczająco dobrą dokładność, że nie można zobaczyć rzeczywistego dryfu na wykresie przestrzeni fazowej. Jednak po prawej stronie widać, że występuje przesunięcie energii, więc jeśli wykonasz wystarczająco długą integrację, ta metoda nie zadziała tak dobrze, jak metoda rozwiązania o właściwościach okresowych. Ale to nasuwa pytanie, jak to wygląda pod względem wydajności w porównaniu z integracją bardzo dokładnie? To jest nieco mniej pewne. W DiffEqBenchmarks.jl można znaleźć pewne testy porównawcze badające to pytanie. Na przykład ten notatnikpatrzy na błąd energii w czasie wykonywania w układzie równań hamiltonowskich z poczwórnego modelu Bosona i pokazuje, że jeśli chcesz naprawdę wysokiej dokładności, to nawet w przypadku dość długich czasów integracji bardziej efektywne jest użycie wysokiej klasy RK lub Runge-Kutta Nystrom ( RKN). Ma to sens, ponieważ aby zaspokoić właściwość symplektyczną, integratorzy rezygnują z pewnej wydajności i prawie muszą być ustalone w określonym przedziale czasowym (prowadzone są badania nad tym ostatnim, ale nie jest to zbyt daleko).

Ponadto zauważ z obu tych komputerów przenośnych, że możesz również zastosować standardową metodę i rzutować ją z powrotem do kolektora rozwiązań na każdym etapie (lub co kilka kroków). To właśnie robią przykłady wykorzystujące wywołanie zwrotne ManifoldProjection DifferentialEquations.jl . Widzisz, że zasady ochrony gwarancji są przestrzegane, ale z dodatkowym kosztem rozwiązania systemu niejawnego na każdym etapie. Możesz również użyć w pełni niejawnego solwera ODE lub pojedynczych macierzy masy, aby dodać równania konserwatorskie, ale efekt końcowy jest taki, że metody te są bardziej kosztowne obliczeniowo jako kompromis.

Podsumowując, klasą problemów, do których chcesz sięgnąć po integrator symplektyczny, są te, które mają rozwiązanie na rozmaitości symplektycznej (systemy Hamiltona), w której nie chcesz inwestować zasobów obliczeniowych w celu uzyskania bardzo dokładnej (tolerancji <1e-12) rozwiązanie i nie potrzebujesz dokładnej energii / itp. ochrona. To podkreśla, że chodzi przede wszystkim o długoterminowe właściwości integracyjne, więc nie powinieneś po prostu gromadzić się wokół nich, chcąc nie chcąc, jak sugeruje to literatura. Ale nadal są one bardzo ważnym narzędziem w wielu dziedzinach, takich jak astrofizyka, w których istnieje długotrwała integracja, którą należy rozwiązać wystarczająco szybko bez absurdalnej dokładności.

Gdzie znajdę integratory symplektyczne? Jakie istnieją integratory symplektyczne?

Generalnie istnieją dwie klasy integratorów symplektycznych. Istnieją symplektyczne integratory Runge-Kutta (które są pokazane w powyższych przykładach) i istnieją niejawne metody Runge-Kutta, które mają właściwość symplektyczną. Jak wspomina @origimbo, symplektyczni integratorzy Runge-Kutta wymagają, aby zapewnić im podzieloną strukturę, aby mogli osobno obsługiwać elementy położenia i pędu. Jednak w przeciwieństwie do komentarza, ukryte metody Runge-Kutty są symplektyczne, nie wymagając tego, ale zamiast tego wymagają rozwiązania układu nieliniowego. Nie jest to takie złe, ponieważ jeśli układ jest niesztywny, ten nieliniowy układ można rozwiązać za pomocą iteracji funkcjonalnej lub przyspieszenia Andersona, ale symplektyczne metody RK powinny nadal być preferowane ze względu na wydajność („

To powiedziawszy, odeint nie ma metod z żadnej z tych rodzin , więc nie jest dobrym wyborem, jeśli szukasz integratorów symplektycznych. W Fortran strona Hairera ma mały zestaw, z którego możesz skorzystać . Mathematica ma kilka wbudowanych . Solwery ODE GSL mają niejawne integratory punktów RK Gaussa, które IIRC są symplektyczne, ale to jedyny powód, dla którego warto zastosować metody GSL.

Ale najbardziej wszechstronny zestaw integratorów symplektycznych można znaleźć w pliku DifferentialEquations.jl w Julii (pamiętaj, że użyto go w powyższych notatnikach). Lista dostępnych symplektycznych metod Runge-Kutta znajduje się na tej stronie i zauważysz, że niejawna metoda punktu środkowego jest również symplektyczna (niejawna metoda trapezoidalna Runge-Kutta jest uważana za „prawie symplektyczną”, ponieważ jest odwracalna). Nie tylko ma największy zestaw metod, ale jest również open source (można zobaczyć kod i jego testy w języku wysokiego poziomu) i ma wiele testów porównawczych . Dobry notatnik wprowadzający do używania go do rozwiązywania problemów fizycznych to ten notatnik z samouczkiem. Ale oczywiście zaleca się, aby rozpocząć pracę z pakietem od pierwszego samouczka ODE .

Ogólnie można znaleźć szczegółową analizę pakietów równań różniczkowych numerycznych w tym poście na blogu . Jest dość szczegółowy, ale ponieważ musi obejmować wiele tematów, każdy robi to mniej szczegółowo niż to, więc możesz poprosić o jego rozszerzenie w jakikolwiek sposób.

Chris Rackauckas
źródło

Przy tej odpowiedzi wydaje mi się, że trafiłem w Jackpot Stack Exchange! To dla mnie idealna odpowiedź, ponieważ rozumiem niektóre z nich od razu i fragmenty, których nie chcę, aby czytać dalej. Naprawdę doceniam czas poświęcony na dokładne pozyskanie tej odpowiedzi, a także podanie innych pomocnych i pouczających linków.

uhoh

Zanim przejdziemy do szczegółów matematycznych, moglibyśmy z grubsza powiedzieć, że symplektyczny oznacza zachowanie objętości , prawda?

Miguel

FTR, powodem, dla którego adaptacyjna Runge-Kutta 5. rzędu działa o wiele lepiej niż RK4, nie jest to, że ma wyższy porządek, ale wybiera bardziej odpowiednie rozmiary stopni. Powodem, dla którego RK4 działa tak źle, jest głównie to, że wielkość kroku jest nieodpowiednio wysoka na perygeum; ten sam solver z połową wielkości kroku dałby znacznie lepsze rozwiązanie. (Po prostu marnowanie czasu na dokładne orbicie wokół apogeum byłoby dużo czasu, gdy nie jest to wcale konieczne).

lewej około

doskonała ekspozycja. Jako pytanie poboczne: OP zaczyna się od odwołania do Pythona - czy istnieją zalecane samouczki / pakiety Pythona zgodne z linkami z przykładami Julii?

Quetzalcoatl

Jedyny znany mi pakiet Pythona dla tego rodzaju integratorów to diffeqpy , gdzie nie jest udokumentowany na README, ale możesz uzyskać dostęp do wszystkich tych samych metod i ponownie napisać to w Pythonie za pomocą tego pakietu.

Chris Rackauckas,

$\mathbf{p}$ $\mathbf{q}$ $\mathcal{H(\mathbf{p},\mathbf{q})}$

\frac{re q}{re t} = + \frac{\partial H.}{\partial p}

$\frac{d\mathbf{q}}{dt}=+\frac{\partial \mathcal{H}}{\partial\mathbf{p}}$

\frac{re p}{re t} = - \frac{\partial H.}{\partial q} .

$\frac{d\mathbf{p}}{dt}=-\frac{\partial \mathcal{H}}{\partial\mathbf{q}}.$

H

$\mathcal{H}$

p (t), q (t) = ϕ_{t} (p (t_{0}), q (t_{0}))

$\mathbf{p}(t),\mathbf{q}(t)= \phi_t(\mathbf{p}(t_0),\mathbf{q}(t_0))$

d p \land d q

$d\mathbf{p}\wedge d\mathbf{q}$

p

$p$

q

$q$ są jednowymiarowe, o czym można pomyśleć, mówiąc, że obszar wewnątrz zamkniętych krzywych w przestrzeni fazowej jest zachowany. Zapewnia to wszelkiego rodzaju ładne właściwości stabilności, ponieważ „kule” trajektorii muszą pozostać „blisko”.

Pod względem numerycznym integrator symplektyczny działa w ten sam sposób, zachowując również ten obszar / dwie formy. To z kolei oznacza, że istnieje zachowany „numeryczny hamiltonian” (który może nie być [czytaj „nie jest”] taki sam jak ten dokładnie). Zauważ, że stabilność nie jest tym samym co dokładność, więc większość zalet metod symplektycznych pojawia się podczas integracji przez bardzo długi czas (np. Twoja metoda może szybko umieścić satelitę po niewłaściwej stronie Ziemi, nigdy nie pozwalając jej na rozpad na to).

origimbo
źródło

Dziękuję Ci za to! Użyję teraz słów powyżej mojej kategorii płac. Kulki n trajektorii są bardziej zagrożone, gdy znajdują się w pobliżu rozgałęzień, takich jak symulacje 3-ciałowe. por. Doedel i in. 2007, Int. J. Bifurcation and Chaos, w. 17, no. 8 (2007) 2625–2677 Jak to zrobiłem? Również ieec.cat/hosted/web-libpoint/papers/…

uhoh

O ile czytelnik nie zna szczegółów matematycznych, wzmianka o stabilności jest myląca, ponieważ zachowanie objętości nie oznacza, że poszczególne trajektorie pozostają blisko.

Miguel

@Maguel Myślę, że jest to jedna z tych sytuacji, w których czytelnik, który nie podąża za matematycznymi szczegółami, jest tak czy inaczej przeklęty, ale jeśli chodzi o trojkę dokładności, stabilności i wydajności obliczeniowej zwykłego numerycy, twierdzę, że podkreślam stabilność korzyści są przydatne. Z przyjemnością przyjmę sugestie dotyczące przepisania, jeśli możesz wymyślić coś lepszego, co nie jest celowo niedokładne.

origimbo

^{2}

$^2$

^{2}

$^2$

@Miguel: Ale kropla cząstek może się rozdzielić na dwie lub więcej części. Jego całkowita objętość musi pozostać stała.

Wolfgang Bangerth