Auf der Suche nach Runge-Kutta 8. Ordnung in C / C ++

10

Ich möchte die Runge-Kutta-Methode 8. Ordnung (89) in einer in C ++ geschriebenen Anwendung für Himmelsmechanik / Astrodynamik unter Verwendung eines Windows-Computers verwenden. Daher frage ich mich, ob jemand eine gute Bibliothek / Implementierung kennt, die dokumentiert und kostenlos zu verwenden ist. Es ist in Ordnung, wenn es in C geschrieben ist, solange keine Kompilierungsprobleme zu erwarten sind.

Bisher habe ich diese Bibliothek (mymathlib) gefunden . Der Code scheint in Ordnung zu sein, aber ich habe keine Informationen zur Lizenzierung gefunden.

Können Sie mir helfen, indem Sie einige der Alternativen aufzeigen, die Sie vielleicht kennen und die zu meinem Problem passen würden?

EDIT:
Ich sehe, dass nicht wirklich so viele C / C ++ - Quellcodes verfügbar sind, wie ich erwartet hatte. Daher wäre auch eine Matlab / Octave-Version in Ordnung (muss noch kostenlos sein).

James C.
quelle

Antworten:

8

Sowohl die GNU Scientific Library (GSL) (C) als auch die Boost Odeint (C ++) verfügen über Runge-Kutta-Methoden 8. Ordnung.

Beide sind Open Source und sollten unter Linux und Mac direkt im Paketmanager verfügbar sein. Unter Windows ist es wahrscheinlich einfacher für Sie, Boost anstelle von GSL zu verwenden.

GSL wird unter der GPL-Lizenz und Boost Odeint unter der Boost-Lizenz veröffentlicht.

Bearbeiten: Ok, Boost Odeint hat NICHT die Runge-Kutta 89-Methode, sondern nur die 78, aber es bietet ein Rezept für die Herstellung beliebiger Runge-Kutta-Stepper.

Methoden 8. Ordnung sind jedoch ziemlich hoch und höchstwahrscheinlich übertrieben für Ihr Problem.

Prince-Dormand bezieht sich auf eine bestimmte Art von Runge-Kutta und ist nicht direkt mit der Reihenfolge verbunden, aber die häufigste ist 45. Matlabs ode45, der empfohlene ODE-Algorithmus, ist eine Prince-Dormand 45-Implementierung. Dies ist der gleiche Algorithmus wie in Boost Odeint Runge_Kutta_Dopri5 implementiert .

LKlevin
quelle
1
Danke für die Antwort. OK, das ist jetzt peinlich. Ich habe mir Boost Odeint angesehen, bevor ich hier gefragt habe, und nur "runge_kutta_fehlberg78" gefunden. Ist das das Richtige? Eigentlich kenne ich die Unterschiede zwischen Methanten in der Praxis nicht, aber ich suchte nach einem RK89 (auch Dormand-Prince genannt, wenn ich im Internet suche). Können Sie bitte Ihre Antwort zu diesem Thema kommentieren oder erweitern? Vielen Dank.
James C
Der Beitrag wurde aktualisiert, um Ihre Fragen zu beantworten. Prince-Dormand 45 wird Ihre Probleme höchstwahrscheinlich gut lösen.
LKlevin
15

Wenn Sie über lange Zeiträume Himmelsmechanik betreiben, spart die Verwendung eines klassischen Runge-Kutta-Integrators keine Energie. In diesem Fall wäre es wahrscheinlich besser, einen symplektischen Integrator zu verwenden. Boost.odeint implementiert auch ein symplektisches Runge-Kutta-Schema 4. Ordnung, das für lange Zeitintervalle besser funktioniert. Soweit ich das beurteilen kann, implementiert GSL keine symplektischen Methoden.

Geoff Oxberry
quelle
Danke für die Antwort. Würde ein symplektisches Runge-Kutta 4. Ordnung bessere Ergebnisse liefern als RKF78, wenn es mit Erdsatelliten (niedrige Umlaufbahn sowie tiefere Weltraumumlaufbahn) verwendet wird, möglicherweise über einen Zeitraum von 1-3 Umlaufbahnen?
James C
@ JamesC Ja. In einer langen Zeit ist die symplektische Methode viel besser.
eccstartup
@eccstartup - Was würdest du für eine lange Zeit hier halten? Weil es so lang sein kann wie eine Umlaufbahn eines Planeten um die Sonne oder einige Umlaufbahnen eines Wettersatelliten um die Erde usw.
James C
@ JamesC Ich habe dieses große Problem nicht beobachtet. Aber für meine Modellprobleme mit vielen berechneten Umlaufbahnen ergeben symplektische Methoden sehr perfekte Umlaufbahnen.
eccstartup
Es ist daher ein Ratschlag, eine Version der impliziten Runge-Kutta-Methode zu programmieren, die viele symplektische Methoden mit beliebig höherer Ordnung enthält.
eccstartup
4

einige Punkte zusammenfassen:

  1. Wenn es sich um eine langfristige Integration eines nicht dissapativen Modells handelt, ist ein symplektischer Integrator genau das, wonach Sie suchen.
  2. Andernfalls sind Runge-Kutta-Nystrom-Methoden, da es sich um eine Bewegungsgleichung handelt, effizienter als eine Transformation in ein System erster Ordnung. Aufgrund von DP gibt es RKN-Methoden hoher Ordnung. Es gibt einige Implementierungen, wie hier in Julia, die dokumentiert sind, und hier ist eine MATLAB- Implementierung .
  3. Runge-Kutta-Methoden hoher Ordnung werden nur benötigt, wenn Sie eine hochgenaue Lösung wünschen. Wenn die Toleranzen niedriger sind, ist ein RK 5. Ordnung wahrscheinlich schneller (bei gleichem Fehler). Das Beste, was Sie tun müssen, wenn Sie dies häufig lösen müssen, ist, eine Reihe verschiedener Methoden zu testen. In dieser Reihe von Benchmarks für 3-Körper-Probleme sehen wir, dass (für denselben Fehler) RK-Methoden hoher Ordnung nur eine marginale Verbesserung der Geschwindigkeit darstellen, obwohl Sie als Fehler -> 0 sehen können, dass die Verbesserung gegenüber Dormand bereits> 5x beträgt -Prinz 45 ( DP5), wenn Sie sich 4 Stellen Genauigkeit ansehen (Toleranzen sind dafür allerdings viel niedriger. Toleranzen sind nur ein Baseballstadion in jedem Problem). Wenn Sie die Toleranzen noch weiter senken, wächst die Verbesserung gegenüber einer RK-Methode höherer Ordnung. Möglicherweise müssen Sie jedoch Zahlen mit höherer Genauigkeit verwenden.
  4. Der 7/8-Algorithmus der Dormand-Prince-Ordnung hat ein anderes Tableau 8. Ordnung als die DP853-Methode von Hairer's dop853und DifferentialEquations.jl DP8(die gleich sind). Die letztere 853-Methode kann nicht in der Standard-Tableau-Version einer Runge-Kutta-Methode implementiert werden, da ihr Fehlerschätzer nicht dem Standard entspricht. Diese Methode ist jedoch viel effizienter und ich würde nicht einmal die Verwendung der älteren Fehlberg 7 / 8- oder DP 7/8-Methoden empfehlen.
  5. Für RK-Methoden hoher Ordnung sind die Verner "Efficient" -Methoden der Goldstandard. Das zeigt sich in den Benchmarks, die ich verlinkt habe. Sie können diese selbst in Boost codieren oder eines der beiden Pakete verwenden, die sie implementieren, wenn Sie dies einfacher möchten (Mathematica oder DifferentialEquations.jl).
Chris Rackauckas
quelle
2

Ich möchte hinzufügen, dass das, was Geoff Oxberry für die langfristige Integration vorschlägt (unter Verwendung symplektischer Integratoren), zwar wahr ist, in einigen Fällen jedoch nicht funktioniert. Insbesondere wenn Sie dissipative Kräfte haben, bewahrt Ihr System die Energie nicht mehr und daher können Sie in diesem Fall nicht auf symplektische Integratoren zurückgreifen. Die Person, die die Frage stellte, sprach von niedrigen Erdumlaufbahnen, und solche Umlaufbahnen weisen einen großen atmosphärischen Widerstand auf, dh eine dissipative Kraft, die die Verwendung solcher symplektischer Integratoren ausschließt.

In diesem speziellen Fall (und in Fällen, in denen Sie keine symplektischen Integratoren verwenden können / keinen Zugriff darauf haben / nicht verwenden möchten) würde ich die Verwendung des Bulirsch-Stoer-Integrators empfehlen, wenn Sie Präzision und Effizienz über lange Zeiträume benötigen. Es funktioniert erfahrungsgemäß gut und wird auch von den Numerical Recipes empfohlen (Press et al., 2007).

viiv
quelle
Nein, empfehlen Sie keine numerischen Rezepte. Insbesondere in den meisten Fällen sollte Burlirsch-Stoer nicht empfohlen werden. Dies ist ein bekanntes Problem mit dem Buch. Hier finden Sie eine Reihe von Gegenargumenten von Spitzenforschern auf diesem Gebiet: uwyo.edu/buerkle/misc/wnotnr.html . Wenn Sie diesbezüglich Benchmarks wünschen, lesen Sie Hairers erstes Buch, in dem BS fast nie gut abschneidet. Eine höhere Ordnung ist nur dann effizienter, wenn die Fehler niedrig genug sind, und wir (und andere) haben Benchmarking durchgeführt, um ganz konsistent zu zeigen, dass sie nur für die Genauigkeit unter Gleitkommazahlen effizient ist.
Chris Rackauckas
Ich kann nicht zu viel für NR sprechen, da ich es hauptsächlich für ODEs verwendet habe, aber es scheint mir, dass die Beschwerden auf der Seite, auf die Sie verlinken, alt sind und von den Autoren von NR in ihrer Antwort (Ende der Seite) angesprochen wurden. aber das ist kein Thema. In Bezug auf die Langzeitintegration von Umlaufbahnen mit hoher Präzision (z. B. 13 bis 14 Stellen), die ich in meiner Antwort erwähnt habe, ist seit langem bewiesen, dass Extrapolationsmethoden gut funktionieren (siehe Kapitel über numerische Integration von Montenbruck & Gill). Neuere Veröffentlichungen verwenden es ebenfalls und es hat sich für mich und andere als zuverlässige und effiziente Methode erwiesen.
viiv
M & G testet es nur gegen dop853 und modernere RK-Methoden hoher Ordnung, wie die von Verner, sind viel effizienter. M & G scheint auch nur anhand von Funktionsbewertungen zu messen, die ein schwacher Indikator für das Timing sind. Es ist auch keine Zeit gegen Runge-Kutta-Nystrom-Methoden, die speziell für ODEs 2. Ordnung sind und um einiges effizienter sind als die RK-Methoden erster Ordnung, die auf ODEs zweiter Ordnung angewendet werden. Mit 13 bis 14 Stellen ist BS bei den meisten Problemen wahrscheinlich wettbewerbsfähig, aber es ist weit von der offensichtlichen Wahl entfernt, und ich habe kein Arbeitsgenauigkeitsdiagramm mit neueren Methoden gesehen, das damit nicht einverstanden ist.
Chris Rackauckas
M & G testet RKN gegen RK und BS und andere gegen RKN (Seiten 123-132 und 151-154) und sagt, dass sie die effizienteste der RK-Methoden sind (ohne Verner, obwohl sie ihn zitieren). BS hat sich bei 13 bis 14 Stellen als effizient erwiesen, was meine Behauptung war. Ich habe gesehen, dass es gegen dop853, ABM (12), Taylor und Standard RK8 getestet wurde und eine gute Leistung erbringt. Ich muss zugeben, dass ich es nicht erneut bei RKN getestet habe, aber nach dem, was ich von M & G sehen kann, ist es zum Beispiel nicht weit von FILG11 entfernt. Ich bin wirklich an Verners RK interessiert und werde mir Ihre obigen Links ansehen. Haben Sie ein Papier, in dem alle getestet werden?
viiv
Ich bin zurückgegangen und habe eine Reihe von Benchmarks bei DiffEqBenchmarks.jl erneut ausgeführt und bin in derodex Regel nicht gut. Zumindest für ODEs 1. Ordnung und für Toleranzen >=1e-13scheint die Extrapolation nicht gut zu funktionieren und ist normalerweise nicht einmal eng. Dies steht im Einklang mit dem obigen Anspruch.
Chris Rackauckas