Was bedeutet "symplektisch" in Bezug auf numerische Integratoren und werden sie von SciPy Odeint verwendet?

25

In diesem Kommentar schrieb ich:

... Standard-SciPy-Integrator, von dem ich annehme, dass er nur symplektische Methoden verwendet.

in dem ich mich auf SciPy's beziehe odeint, das entweder eine "nicht steife (Adams) Methode" oder eine "steife (BDF) Methode" verwendet. Nach der Quelle :

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

Hier ist ein Beispiel, in dem ich einen Satelliten drei Monate lang um die Erde bewege, um zu zeigen, dass er erwartungsgemäß funktioniert.

Ich glaube, dass nicht-symplektische Integratoren die unerwünschte Eigenschaft haben, dass sie dazu neigen, keine Energie (oder andere Größen) zu sparen, und daher beispielsweise in der Orbitalmechanik unerwünscht sind. Aber ich bin nicht ganz sicher, was einen symplektischen Integrator symplektisch macht.

Ist es möglich, die Eigenschaft (die einen symplektischen Integrator symplektisch macht) auf einfache und (ziemlich) leicht verständliche, aber nicht ungenaue Weise zu erklären? Ich frage unter dem Gesichtspunkt, wie der Integrator intern funktioniert und nicht wie er beim Testen abschneidet.

Und ist mein Verdacht richtig, dass odeintnur symplektische Integratoren eingesetzt werden?

uhoh
quelle
4
Als Faustregel sollten Sie nur hoffen, dass ein Black-Box-Integrator symplektisch ist, wenn Sie Positions- und Impulsgleichungen trennen müssen.
Origimbo
@origimbo Danke. Diese funktionieren, und es sieht so aus, als wäre odeintes ein Python-Wrappoer für ziemlich alte, etablierte und gut referenzierte Quellcodes (bearbeitete Frage, Verweise auf ODEPACK und LSODA), obwohl ich es mit Sicherheit zugebe, im Black-Box-Modus zu arbeiten. Mein verknüpftes Beispiel zeigt, dass der 6D-Zustandsvektor aus drei Positionen und drei Geschwindigkeiten besteht.
Uhoh
11
Die ODE-Integratoren in ODEPACK und LSODA sind keine symplektischen Integratoren.
Brian Borchers
2
Hier ist ein Beispiel, das zwei sehr einfache Löser vergleicht: Euler und Symplectic Euler: idontgetoutmuch.wordpress.com/2013/08/06/… .
Idontgetoutmuch
2
Das Buch von Hairer, Nørsett und Wanner gibt eine gute Erklärung für symplektische Methoden. Betrachten Sie insbesondere Abbildung 16.1 und die Abbildungen hier .
JM

Antworten:

47

Lassen Sie mich mit Korrekturen beginnen. Nein, odeinthat keine symplektischen Integratoren. Nein, symplektische Integration bedeutet keine Energieeinsparung.

Was bedeutet symplektisch und wann sollten Sie es anwenden?

Was bedeutet symplektisch? Symplektisch bedeutet, dass die Lösung auf einer symplektischen Mannigfaltigkeit besteht. Eine symplektische Mannigfaltigkeit ist eine Lösungsmenge, die durch eine 2-Form definiert ist. Die Details symplektischer Mannigfaltigkeiten klingen wahrscheinlich wie mathematischer Unsinn, und stattdessen besteht der Kern der Sache darin, dass eine direkte Beziehung zwischen zwei Mengen von Variablen auf einer solchen Mannigfaltigkeit besteht. Der Grund, warum dies für die Physik wichtig ist, besteht darin, dass Hamilton'sche Gleichungen natürlich besagen, dass sich die Lösungen auf einer symplektischen Mannigfaltigkeit im Phasenraum befinden, wobei die natürliche Aufspaltung die Positions- und Impulskomponenten sind. Für die wahre Hamiltonsche Lösung ist dieser Phasenraumpfad eine konstante Energie.

Ein symplektischer Integrator ist ein Integrator, dessen Lösung auf einer symplektischen Mannigfaltigkeit beruht. Aufgrund von Diskretisierungsfehlern wird beim Lösen eines Hamilton-Systems nicht genau die richtige Flugbahn auf dem Verteiler angezeigt. Stattdessen selbst , die Trajektorie gestörtes für die Ordnung n von der tatsächlichen Flugbahn. Dann gibt es eine lineare Drift aufgrund eines numerischen Fehlers dieser Trajektorie über die Zeit. Normale Integratoren neigen zu einer quadratischen (oder größeren) Drift und haben keine guten globalen Garantien für diesen Phasenraumpfad (nur lokal).O(Δtn)n

Dies bedeutet in der Regel, dass symplektische Integratoren die Langzeitmuster aufgrund der mangelnden Drift und der fast garantierten Periodizität besser erfassen als normale Integratoren. Dieses Notizbuch zeigt diese Eigenschaften beim Kepler-Problem gut an . Das erste Bild zeigt, wovon ich mit der periodischen Natur der Lösung spreche.

Bildbeschreibung hier eingeben

Dies wurde mit dem symplektischen Integrator 6. Ordnung von Kahan und Li von DifferentialEquations.jl gelöst . Sie können sehen, dass die Energie nicht genau konserviert ist, aber ihre Variation hängt davon ab, wie weit der gestörte Lösungsverteiler vom tatsächlichen Verteiler entfernt ist. Da sich die numerische Lösung selbst auf einer symplektischen Mannigfaltigkeit befindet, ist sie in der Regel fast genau periodisch (mit einer gewissen linearen numerischen Abweichung, die Sie sehen können), was sie für die Langzeitintegration sehr nützlich macht. Wenn Sie dasselbe mit RK4 tun, können Sie eine Katastrophe bekommen:

Bildbeschreibung hier eingeben

Sie können sehen, dass das Problem darin besteht, dass die numerische Lösung keine echte Periodizität aufweist und daher dazu neigt, über die Zeit hinaus zu driften.

Dies unterstreicht den wahren Grund, sich für symplektische Integratoren zu entscheiden: Symplektische Integratoren eignen sich gut für Langzeitintegrationen bei Problemen mit der symplektischen Eigenschaft (Hamiltonsche Systeme) . Gehen wir also ein paar Dinge durch. Beachten Sie, dass Sie auch bei einem symplektischen Problem nicht immer symplektische Integratoren benötigen. In diesem Fall kann eine adaptive Runge-Kutta-Methode 5. Ordnung gut funktionieren. Hier ist Tsit5:

Bildbeschreibung hier eingeben

Beachten Sie zwei Dinge. Erstens wird die Genauigkeit so gut, dass Sie die tatsächliche Drift im Phasenraumdiagramm nicht sehen können. Auf der rechten Seite sehen Sie jedoch, dass es diese Energiedrift gibt. Wenn Sie also eine ausreichend lange Integration durchführen, funktioniert diese Methode nicht so gut wie die Lösungsmethode mit den periodischen Eigenschaften. Aber das wirft die Frage auf, wie es sich in Bezug auf Effizienz auswirkt, wenn man nur extrem genau integriert. Nun, das ist ein bisschen weniger sicher. In DiffEqBenchmarks.jl finden Sie einige Benchmarks, die diese Frage untersuchen. Zum Beispiel dieses NotebookBetrachtet man den Energiefehler gegen die Laufzeit eines Hamiltonschen Gleichungssystems anhand eines Vierfach-Boson-Modells, so zeigt sich, dass es auch bei längeren Integrationszeiten effizienter ist, nur ein RK höherer Ordnung oder ein Runge-Kutta-Nystrom zu verwenden ( RKN) Methode. Dies ist sinnvoll, da die Integratoren zur Erfüllung der symplektischen Eigenschaft auf einen gewissen Wirkungsgrad verzichten und einen festen Zeitschritt einhalten müssen (es gibt einige Untersuchungen, die Fortschritte bei letzterem machen, aber es ist nicht sehr weit).

Beachten Sie außerdem, dass Sie bei beiden Notebooks auch einfach eine Standardmethode anwenden und diese bei jedem Schritt (oder bei einigen Schritten) wieder auf die Lösungsübersicht projizieren können. Dies ist, was die Beispiele mit dem ManifoldProjection-Rückruf " DifferentialEquations.jl" tun. Sie sehen, dass die Schutzgesetze eingehalten werden, jedoch mit zusätzlichen Kosten für die Lösung eines impliziten Systems in jedem Schritt. Sie können auch einen vollständig impliziten ODE-Löser oder singuläre Massenmatrizen verwenden, um Erhaltungsgleichungen hinzuzufügen. Das Endergebnis ist jedoch, dass diese Methoden als Kompromiss rechenintensiver sind.

Zusammenfassend lässt sich sagen, dass die Klasse von Problemen, bei denen Sie nach einem symplektischen Integrator suchen, diejenigen sind, die eine Lösung für eine symplektische Mannigfaltigkeit (Hamilton-Systeme) haben, bei der Sie die Rechenressourcen nicht investieren möchten, um eine sehr genaue (Toleranz <1e-12) zu haben. Lösung und brauchen keine exakte Energie / etc. Erhaltung. Dies unterstreicht, dass es sich um langfristige Integrationseigenschaften handelt. Sie sollten sich also nicht einfach wie in einigen Literaturstellen beschrieben zu ihnen begeben. Aber sie sind immer noch ein sehr wichtiges Werkzeug in vielen Bereichen wie der Astrophysik, in denen Sie lange Zeitintegrationen haben, die Sie schnell genug lösen müssen, ohne eine absurde Genauigkeit zu haben.

Wo finde ich symplektische Integratoren? Welche Art von symplektischen Integratoren gibt es?

Im Allgemeinen gibt es zwei Klassen symplektischer Integratoren. Es gibt die symplektischen Runge-Kutta-Integratoren (die in den obigen Beispielen gezeigt werden) und es gibt implizite Runge-Kutta-Methoden, die die symplektische Eigenschaft haben. Wie @origimbo erwähnt, müssen die symplektischen Runge-Kutta-Integratoren mit einer partitionierten Struktur versehen werden, damit sie die Positions- und Momentum-Teile getrennt verarbeiten können. Entgegen dem Kommentar sind die impliziten Runge-Kutta-Methoden symplektisch, ohne dies zu erfordern, sondern erfordern die Lösung eines nichtlinearen Systems. Das ist nicht so schlimm, denn wenn das System nicht steif ist, kann dieses nicht lineare System durch funktionale Iteration oder Anderson-Beschleunigung gelöst werden.

Das sei gesagt, odeint hat keine Methoden von einem dieser Familien , so dass es nicht eine gute Wahl , wenn Sie für symplektischer Integratoren suchen. In Fortran gibt es auf der Website von Hairer ein kleines Set, das Sie verwenden können . Mathematica hat einige eingebaut . Die GSL-ODE-Löser haben implizite RK-Gauß- Punktintegratoren, die IIRC-Symplektiker sind, aber das ist ungefähr der einzige Grund, die GSL-Methoden zu verwenden.

Aber die umfassendste Reihe von symplektischer Integratoren finden Sie in DifferentialEquations.jl in Julia ( man erinnere sich dies für die Notebooks oben verwendet wurde). Die Liste der verfügbaren symplektischen Runge-Kutta-Methoden finden Sie auf dieser Seite. Sie werden feststellen, dass die implizite Mittelpunktmethode auch symplektisch ist (die implizite Runge-Kutta-Trapezoid-Methode wird als "fast symplektisch" angesehen, weil sie reversibel ist). Es verfügt nicht nur über die meisten Methoden, sondern ist auch Open-Source-fähig (Sie können den Code und seine Tests in einer Hochsprache anzeigen) und verfügt über zahlreiche Benchmarks . Ein gutes Einführungsnotizbuch zur Lösung physikalischer Probleme ist dieses Notizbuch. Aber natürlich wird empfohlen, dass Sie mit dem Paket durch das erste ODE-Tutorial beginnen .

Im Allgemeinen finden Sie in diesem Blogbeitrag eine detaillierte Analyse numerischer Differentialgleichungssuiten . Es ist ziemlich detailliert, aber da es viele Themen abdecken muss, ist jedes weniger detailliert. Sie können also jederzeit darum bitten, dass es in irgendeiner Weise erweitert wird.

Chris Rackauckas
quelle
10
Mit dieser Antwort scheine ich den Stack Exchange Jackpot erreicht zu haben! Dies ist die perfekte Antwort für mich, da ich einige davon sofort verstehe und die Teile, die ich nicht weiter lesen möchte. Ich freue mich sehr über die Zeit, die Sie gebraucht haben, um diese Antwort gründlich zu finden, und auch über andere hilfreiche und lehrreiche Links.
Uhoh
Bevor wir uns den mathematischen Details zuwenden, können wir grob sagen, dass Symplektikum die Erhaltung des Volumens bedeutet , nicht wahr?
Miguel
2
FTR, der Grund, warum das adaptive Runge-Kutta 5. Ordnung so viel besser abschneidet als RK4, ist nicht, dass es eine höhere Ordnung hat, sondern dass es geeignetere Schrittgrößen auswählt. Der Grund, warum RK4 so schlecht abschneidet, liegt hauptsächlich darin, dass die Schrittweite am Perigäum unangemessen hoch ist. der gleiche Löser mit der halben Schrittgröße würde eine viel bessere Lösung ergeben. (Nur, es würde viel Zeit verschwenden, die Umlaufbahn fein um das Apogäum
aufzulösen
1
hervorragende ausstellung. Als Nebenfrage: Das OP beginnt mit einem Verweis auf Python - gibt es empfohlene Python-Tutorials / -Pakete gemäß den verknüpften Julia-Beispielen?
Quetzalcoatl
1
Das einzige mir bekannte Python-Paket für diese Art von Integratoren ist diffeqpy , wo es nicht in der README- Datei dokumentiert ist, aber Sie können auf alle diese Methoden zugreifen und dies in Python mit diesem Paket neu schreiben.
Chris Rackauckas
14

pqH(p,q)

dqdt=+Hp
dpdt=Hq.
H
p(t),q(t)=ϕt(p(t0),q(t0))
dpdqpqWenn Sie eindimensional sind, können Sie davon ausgehen, dass der Bereich innerhalb geschlossener Kurven im Phasenraum erhalten bleibt. Dies stellt alle Arten von guten Stabilitätseigenschaften sicher, da "Kugeln" von Flugbahnen "nahe" zueinander bleiben müssen.

In Bezug auf die Numerik verhält sich ein symplektischer Integrator auf die gleiche Art und Weise, wobei auch diese Flächen- / Zweiformen erhalten bleiben. Dies bedeutet wiederum, dass es einen konservierten "numerischen Hamilton-Operator" gibt (der möglicherweise nicht mit dem exakten identisch ist). Beachten Sie, dass Stabilität nicht mit Genauigkeit gleichzusetzen ist, sodass die meisten Vorteile symplektischer Methoden bei einer Integration über einen sehr langen Zeitraum erzielt werden (z. B. kann Ihre Methode einen Satelliten schnell auf der falschen Seite der Erde platzieren, ohne dass er verfällt) es).

origimbo
quelle
Danke dafür! Ich werde jetzt Wörter verwenden, die über meiner Gehaltsstufe liegen. n-Kugeln von Trajektorien sind gefährdeter, wenn sie sich in der Nähe von Gabelungen befinden, z. B. bei 3-Körper-Simulationen. vgl. Doedel et al. 2007, Int. J. Bifurcation and Chaos, 17, No. 8 (2007) 2625–2677 Wie habe ich es gemacht? Also ieec.cat/hosted/web-libpoint/papers/…
uhoh
2
Sofern dem Leser die mathematischen Details nicht bekannt sind, ist die Erwähnung der Stabilität irreführend, da die Erhaltung des Volumens nicht bedeutet, dass die einzelnen Flugbahnen eng bleiben.
Miguel
1
@ Miguel Ich denke, dies ist eine dieser Situationen, in denen der Leser, der den mathematischen Details nicht folgt, in keiner Weise beeinträchtigt ist, aber in Bezug auf die Genauigkeit, Stabilität und Recheneffizienz des üblichen Numerikers würde ich argumentieren, dass die Stabilität betont wird Vorteile ist nützlich. Gerne nehme ich Vorschläge für ein Umschreiben entgegen, wenn Sie sich etwas Besseres einfallen lassen, das nicht absichtlich ungenau ist.
Origimbo
22
1
@ Miguel: Aber der Partikelklumpen darf sich in zwei oder mehr Teile teilen. Sein Gesamtvolumen hat nur konstant zu bleiben.
Wolfgang Bangerth