Warum werden Runge-Kutta-Methoden höherer Ordnung nicht häufiger angewendet?

17

Ich war nur neugierig, warum Runge-Kutta-Methoden höherer Ordnung (dh größer als 4) (zumindest meines Wissens) so gut wie nie diskutiert / angewendet werden. Ich verstehe, dass es eine längere Rechenzeit pro Schritt erfordert (z. B. RK14 mit eingebettetem Schritt 12. Ordnung ), aber gibt es noch andere Nachteile bei der Verwendung von Runge-Kutta-Methoden höherer Ordnung (z. B. Stabilitätsprobleme)? Wären solche Methoden höherer Ordnung bei Gleichungen mit stark oszillierenden Lösungen auf extremen Zeitskalen nicht in der Regel vorzuziehen?

Mathews24
quelle
2
Ich denke, das ist eine sehr subjektive Frage. Der größte Nachteil, den Sie bereits festgestellt haben, sind die Berechnungskosten. Im Allgemeinen versuchen wir, zwischen Genauigkeit und Rechenzeit auszugleichen. Wenn Menschen in PDEs von höherer Ordnung sprechen, denken sie im Allgemeinen über die 3. oder 4. Ordnung nach. Und auch der Zeitsprung bleibt in der gleichen Reihenfolge.
Vikram
3
Bei PDE ist ein Präzisionsschema hoher Ordnung für die zeitliche Abhängigkeit nicht sinnvoll, wenn die räumliche Genauigkeit schlechter ist. Tatsächlich liegt die Genauigkeit der räumlichen Abhängigkeit meistens in der 2. oder 3. Ordnung, insbesondere bei der Arbeit an unstrukturierten Netzen. Die Menschen müssen die globale Fehlerkürzung mit den geringsten Kosten kontrollieren, daher wird Runge-Kutta in bestimmten Fällen mit einer ausreichend hohen Genauigkeitsreihenfolge betrachtet.
Dienstag,
@tqviet Wenn für die räumlichen Ableitungen Rückwärts- oder Zentraldifferenznäherungen bis zur Ordnung 8 verwendet würden, wäre RK8 geeignet, nein? Gibt es im Allgemeinen Genauigkeits- oder Stabilitätsprobleme bei der Verwendung solcher Näherungen der räumlichen Ableitungen mit endlichen Differenzen hoher Ordnung?
Mathews24
1
@ Mathews24: Ich erwähnte nicht die Stabilität, die stark von der Gleichung abhängt. Wenn eine hochgenaue Regelung räumliche Abhängigkeit angewandt wird, nehmen wir an RK zeitliche Abhängigkeit mit mindestens der gleichen Reihenfolge an Genauigkeit, aber die Stabilitätsbedingung kann einen kleineren Wert von erfordern . Δt
Dienstag,

Antworten:

17

Es gibt Tausende von Artikeln und Hunderte von Codes, die Runge-Kutta-Methoden fünfter Ordnung oder höher verwenden. Beachten Sie, dass der in MATLAB am häufigsten verwendete explizite Integrator ODE45 ist, der die Lösung mithilfe einer Runge-Kutta-Methode 5. Ordnung erweitert.

Beispiele für weit verbreitete Runge-Kutta-Methoden höherer Ordnung

Die Arbeit von Dormand & Prince, die eine Methode 5. Ordnung angibt, hat laut Google Scholar über 1700 Zitate . Die meisten davon sind Papiere, die ihre Methode verwenden, um ein Problem zu lösen. Das Cash-Karp-Methodenpapier hat über 400 Zitate . Die wahrscheinlich am häufigsten verwendete Methode höherer Ordnung als 5 ist die Methode 8. Ordnung von Prince-Dormand mit über 400 Zitaten in Google Scholar . Ich könnte viele andere Beispiele nennen; und denken Sie daran, dass viele (wenn nicht die meisten) Leute, die diese Methoden anwenden, die Papiere nie zitieren.

Beachten Sie auch, dass Extrapolationsmethoden höherer Ordnung und verzögerte Korrekturmethoden Runge-Kutta-Methoden sind .

Methoden höherer Ordnung und Rundungsfehler

Wenn Ihre Genauigkeit durch Rundungsfehler eingeschränkt ist , sollten Sie eine Methode höherer Ordnung verwenden . Dies liegt daran, dass Methoden höherer Ordnung weniger Schritte erfordern (und weniger Funktionsauswertungen, obwohl es mehr Auswertungen pro Schritt gibt), sodass weniger Rundungsfehler auftreten. Sie können dies mit einfachen Experimenten leicht selbst überprüfen. Es ist eine gute Hausaufgabe für einen ersten Kurs in numerischer Analyse.

Methoden zehnter Ordnung sind in der Arithmetik mit doppelter Genauigkeit äußerst nützlich. Im Gegenteil, wenn alles, was wir hätten, die Euler-Methode wäre, wäre der Rundungsfehler ein großes Problem, und wir würden für viele Probleme, bei denen Löser höherer Ordnung gut funktionieren, sehr hochpräzise Gleitkommazahlen benötigen.

Methoden hoher Ordnung können genauso stabil sein

EINB

Methoden höherer Ordnung in der Himmelsmechanik

Du fragst

Wären solche Methoden höherer Ordnung bei Gleichungen mit stark oszillierenden Lösungen auf extremen Zeitskalen nicht in der Regel vorzuziehen?

Du hast genau recht! Ein Paradebeispiel dafür ist die Himmelsmechanik. Ich bin kein Experte auf diesem Gebiet. Aber diese Arbeit vergleicht zum Beispiel Methoden für die Himmelsmechanik und berücksichtigt nicht einmal eine niedrigere Ordnung als 5. Sie kommt zu dem Schluss, dass Methoden der Ordnung 11 oder 12 oft am effizientesten sind (mit der Prince-Dormand-Methode der Ordnung 8 auch oft sehr effizient).

David Ketcheson
quelle
Ketchson: Können Sie bitte Beweise oder Erklärungen zu dieser Aussage liefern: "Hochwertige Extrapolations- und verzögerte Korrekturmethoden sind Runge-Kutta-Methoden"? Insbesondere die "aufgeschobenen Korrekturmethoden". Vielen Dank.
Dienstag,
@ David Ketcheson Können Sie diskutieren, wie sich Ihre Antwort bei Verwendung validierter (verifizierter) Rechentechniken wie nach außen gerundetem Intervall oder radialer Arithmetik ändern würde? Wie wäre es, wenn ein nach außen gerundetes Intervall oder eine radiale Arithmetik mit einer höheren Genauigkeit als der doppelten verwendet würden? Was passiert mit dem Wickeln und der Abhängigkeit, wenn die Runge-Kutta-Reihenfolge erhöht wird, und sagen wir zum Spaß, dass die ODE sehr steif ist?
Mark L. Stone
@ MarkL.Stone Das sind ganz andere Fragen. Wenn Sie sie fragen möchten, stellen Sie sie bitte als separate Fragen. Ich bin jedoch kein Experte in diesen Dingen und werde nicht in der Lage sein zu antworten.
David Ketcheson
1
@tqviet In diesem Artikel finden Sie eine Erklärung.
David Ketcheson
12

Solange Sie standardmäßige Gleitkomma-Arithmetik mit doppelter Genauigkeit verwenden, sind keine Methoden hoher Ordnung erforderlich, um eine Lösung mit hoher Genauigkeit in einer angemessenen Anzahl von Schritten zu erhalten. In der Praxis stelle ich fest, dass die Genauigkeit der Lösung normalerweise durch die Gleitkommadarstellung mit doppelter Genauigkeit auf einen relativen Fehler von 1,0e-16 beschränkt ist und nicht durch die Anzahl / Länge der Schritte, die mit RKF45 ausgeführt werden.

Wenn Sie zu einem Gleitkomma-Arithmetikschema mit einer höheren Genauigkeit als der doppelten wechseln, lohnt es sich möglicherweise, eine Methode 10. Ordnung zu verwenden.

Brian Borchers
quelle
5
Ich denke, diese Antwort ist irreführend. Methoden höherer Ordnung führen zu viel weniger Rundungsfehlern, wohingegen Methoden niedriger Ordnung unter Rundungsfehlern leiden, die dominieren, wenn die erforderliche Genauigkeit groß ist oder das Zeitintervall lang ist. Siehe meine Antwort unten.
David Ketcheson
2
Der Punkt ist, dass Sie in Gleitkommazahlen mit doppelter Genauigkeit nicht einmal eine Lösung mit einer relativen Genauigkeit von 1,0e-16 darstellen können. In vielen praktischen Situationen bringt Sie der gute alte RKF45 über den Zeitraum, an dem Sie interessiert sind, auf dieses Genauigkeitsniveau, ohne dass winzige Schritte erforderlich sind. Es ist möglicherweise keine gute Wahl für steife Systeme oder Situationen, in denen ein symplektischer Integrator erforderlich ist, aber eine Runge-Kutta-Methode höherer Ordnung ist auch in diesen Situationen keine gute Lösung. Ich bin damit einverstanden, dass für sehr lange Zeiträume Runge-Kutta-Methoden höherer Ordnung sinnvoll sein können.
Brian Borchers
10

Um die ausgezeichnete Antwort von Brian Borcher zu ergänzen, lassen viele reale Anwendungen sehr steife ODEs oder DAEs zu. Intuitiv erfahren diese Probleme nicht glatte, abrupte Änderungen über die Zeit, so dass sie besser mit Polynomen niedriger Ordnung modelliert werden können, die sich über kurze Schrittgrößen verteilen, als mit Polynomen hoher Ordnung, die sich über lange Schrittgrößen erstrecken. Außerdem erfordert die Stabilität häufig die Verwendung impliziter Methoden, für die der Rechenaufwand für Methoden höherer Ordnung viel steiler ist.

Strenger sind Methoden höherer Ordnung für steife Probleme weniger stabil als Methoden niedrigerer Ordnung. Wir haben zum Beispiel die Dahlquist-Barrieren für lineare Mehrschrittverfahren.

r2

Ähnliche (aber weitaus kompliziertere) Aussagen können für die L-Stabilität in den RK-Formeln gemacht werden. In allen Fällen führt die Erhöhung der Reihenfolge häufig nicht immer zu genaueren Lösungen. Das Folgende ist ein Auszug aus der wegweisenden Arbeit von Prothero und Robinson von 1974:

Bei der Verwendung von A-stabilen einstufigen Methoden zur Lösung großer Systeme von steifen nichtlinearen Differentialgleichungen haben wir festgestellt, dass
(a) einige A-stabile Methoden sehr instabile Lösungen ergeben und
(b) die Genauigkeit der Lösungen, die mit den Gleichungen erhalten werden steif scheint oft unabhängig von der Reihenfolge der verwendeten Methode zu sein.

Weitere Informationen zu diesem Thema finden Sie im klassischen Text von Hairer & Wanner, "Lösung gewöhnlicher Differentialgleichungen II: Steife und Differential - Algebraische Probleme", 1991.

In der Praxis werden steife Gleichungen fast immer mit der Trapezregel oder der TR-BDF2-Formel (ode23t- und ode23tb-Funktionen in MATLAB) gelöst. Beides sind implizite Methoden zweiter Ordnung. Natürlich können wir, wenn Stabilität kein Thema ist (z. B. bei nicht steifen Gleichungen), aus einer Reihe von Optionen wählen. RK45 ist die häufigste Wahl.

Richard Zhang
quelle
Sehr interessant. Gibt es eine (intuitive) Erklärung dafür, warum die Reihenfolge kleiner oder gleich 2 sein muss, damit es sich um eine A-stabile Mehrschrittmethode handelt? Und nur um zu verdeutlichen, wenn Sie sagen, dass ähnliche Aussagen für RK-Formeln gemacht werden können, ist es erneut in der Größenordnung 2?
Mathews24
Für Runge-Kutta-Methoden gibt es A-stabile Methoden in beliebiger Reihenfolge.
David Ketcheson
@DavidKetcheson Ja, aber sie sind nicht stark A-stabil (dh L-stabil). Bei der Lösung von DAEs treten viele Probleme auf, z. B. die Simulation einfacher Transistorschaltungen. In der Tat ist TR dafür berüchtigt, in SPICE künstliche Klingeleffekte zu verursachen, was die Entwicklung von TR-BDF2 motivierte.
Richard Zhang
@DavidKetcheson Zum Nachschlagen siehe doi.org/10.1090/S0025-5718-1974-0331793-2 . Der Begriff der A-Stabilität ist für DAEs nicht stark genug, und Methoden mit A-Stabilität höherer Ordnung führen bei der Lösung von DAEs häufig zu seltsamen Ergebnissen.
Richard Zhang
Sicher, aber es geht nicht um DAEs oder um mehrstufige Methoden.
David Ketcheson
9

Das Benchmark-Setup

In der Julia-Software DifferentialEquations.jl haben wir viele Methoden höherer Ordnung implementiert, einschließlich der Feagin-Methoden. Sie können es in unserer Liste der Methoden sehen , und dann gibt es Tonnen von anderen, die Sie als mitgelieferte Tableaus verwenden können . Da alle diese Methoden zusammengesetzt sind, können Sie sie problemlos miteinander vergleichen. Sie können die Benchmarks, die ich online habe, hier sehen und sehen, dass es sehr einfach ist, viele verschiedene Algorithmen zu vergleichen. Wenn Sie sich also ein paar Minuten Zeit nehmen möchten, um die Benchmarks auszuführen, sollten Sie es versuchen. Hier ist eine Zusammenfassung der Ergebnisse.

Zunächst ist zu beachten, dass bei Betrachtung der einzelnen Benchmarks festgestellt wird, dass unsere DP5(Dormand-Prince Order 5) und DP8Methoden schneller sind als die Hairer Fortran-Codes ( dopri5und dop853). Daher sind diese Implementierungen sehr gut optimiert . Diese zeigen, dass, wie in einem anderen Thread angemerkt, die Überbeanspruchung der Dormand-Prince-Methoden darauf zurückzuführen ist, dass die Methoden bereits geschrieben wurden und nicht darauf, dass sie immer noch die besten sind. Der tatsächliche Vergleich zwischen den am besten optimierten Implementierungen besteht also zwischen den Tsitorous-Methoden, den Verner-Methoden und den Feagin-Methoden von DifferentialEquations.jl.

Die Ergebnisse

Im Allgemeinen haben die Verfahren mit einer Bestellung von mehr als 7 zusätzliche Berechnungskosten, die in der Regel aufgrund der gewählten Toleranzen nicht aufgewogen werden. Ein Grund dafür ist, dass die Koeffizientenauswahl für Methoden niedrigerer Ordnung optimiert ist (sie haben kleine "Hauptabschneidungsfehlerkoeffizienten", die wichtiger sind, wenn Sie nicht asymtopisch klein sind). Sie können sehen, dass die Methoden Verner Efficient 6 und 7 bei vielen Problemen wie hier sehr gut funktionieren, aber Methoden wie die Verner Efficient 8 können eine geringere Steigung aufweisen. Dies liegt daran, dass sich die "Gewinne" höherer Ordnung bei niedrigeren Toleranzen addieren, sodass es immer eine Toleranz gibt, bei der die Methoden höherer Ordnung effizienter sind.

Die Frage ist aber dann, wie niedrig? In einer gut optimierten Implementierung wird dies aus zwei Gründen ziemlich niedrig. Der erste Grund ist, dass Methoden niedrigerer Ordnung so genannte FSAL implementieren (first same as last). Diese Eigenschaft bedeutet, dass die Methoden niedrigerer Ordnung eine Funktionsbewertung aus dem vorherigen Schritt im nächsten Schritt wiederverwenden und somit effektiv eine Funktionsbewertung weniger haben. Wenn dies richtig angewendet wird, führt eine Methode fünfter Ordnung (Tsitorous oder Dormand-Prince) tatsächlich 5 Funktionsbewertungen durch, anstatt der 6, die die Tableaus vorschlagen würden. Dies gilt auch für die Verner 6-Methode.

Der andere Grund ist auf Interpolationen zurückzuführen. Ein Grund für die Verwendung einer Methode sehr hoher Ordnung besteht darin, weniger Schritte zu unternehmen und einfach Zwischenwerte zu interpolieren. Um jedoch die Zwischenwerte zu erhalten, benötigt die Interpolationsfunktion möglicherweise mehr Funktionsauswertungen, als für den Schritt verwendet wurden. Wenn Sie sich die Verner-Methoden ansehenbenötigt die Order 8-Methode 8 zusätzliche Funktionsauswertungen, um eine Order 8-Interpolante zu erhalten. Oftmals liefern die Methoden niedriger Ordnung eine "freie" Interpolation, zum Beispiel haben die meisten Methoden 5. Ordnung eine freie Interpolation 4. Ordnung (keine zusätzlichen Funktionsbewertungen). Wenn Sie also Zwischenwerte benötigen (die Sie für ein gutes Diagramm benötigen, wenn Sie eine Methode hoher Ordnung verwenden), fallen zusätzliche versteckte Kosten an. Berücksichtigen Sie, dass diese interpolierten Werte für die Ereignisbehandlung und das Lösen von Verzögerungsdifferentialgleichungen wirklich wichtig sind, und Sie sehen, warum die zusätzlichen Interpolationskosten eine Rolle spielen.

Was ist also mit den Feagin-Methoden?

Sie werden also feststellen, dass die Feagin-Methoden in den Benchmarks verdächtig fehlen. Sie sind in Ordnung, Konvergenztests arbeiten mit willkürlichen Präzisionszahlen usw., aber um sie tatsächlich zu verbessern, müssen Sie nach absurd niedrigen Toleranzen fragen. Zum Beispiel fand ich in unveröffentlichten Benchmarks, dass die Feagin14Outperformance Vern9(die Verner Efficient Method 9. Ordnung) bei Toleranzen wie 1e-30. Für Anwendungen mit chaotischer Dynamik (wie bei den Pleides- oder 3-Körper-Astrophysik-Problemen) möchten Sie möglicherweise diese Genauigkeit aufgrund der empfindlichen Abhängigkeit (Fehler in chaotischen Systemen verstärken sich schnell). Die meisten Leute rechnen jedoch wahrscheinlich mit Gleitkommazahlen mit doppelter Genauigkeit, und ich habe keinen Benchmark gefunden, bei dem sie in diesem Bereich der Toleranz eine Outperformance erzielen.

Darüber hinaus gibt es keine Interpolation, die mit den Feagin-Methoden einhergeht. Also setze ich einfach eine Hermite-Interpolation dritter Ordnung auf sie, so dass eine existiert (und das funktioniert überraschend gut). Wenn es jedoch keine Standardinterpolationsfunktion gibt, können Sie die rekursive Hermite-Methode ausführen (verwenden Sie diese Interpolation, um den Mittelpunkt zu erhalten, und dann eine Interpolation 5. Ordnung usw.), um eine Interpolation höherer Ordnung zu erhalten. Dies ist jedoch sehr kostspielig und das Ergebnis Die Interpolation hat nicht unbedingt einen niedrigen prinzipiellen Trunkierungsfehlerterm (daher ist es nur gut, wenn er dtwirklich klein ist, was genau das Gegenteil des von uns gewünschten Falls ist!). Wenn Sie also jemals eine wirklich gute Interpolation benötigen, die zu Ihrer Genauigkeit passt, müssen Sie zumindest auf so etwas wie zurückgreifen Vern9.

Hinweis zur Extrapolation

Es ist zu beachten, dass Extrapolationsmethoden einfach Algorithmen zum Erzeugen von Runge-Kutta-Methoden beliebiger Ordnung sind. Für ihre Reihenfolge machen sie jedoch mehr Schritte als nötig und haben hohe prinzipielle Trunkierungsfehlerkoeffizienten, so dass sie bei einer gegebenen Reihenfolge nicht so effizient sind wie eine gut optimierte RK-Methode. In Anbetracht der vorherigen Analyse bedeutet dies jedoch, dass es einen Bereich mit äußerst geringer Toleranz gibt, in dem diese Methoden eine bessere Leistung erbringen als die "bekannten" RK-Methoden. Aber in jedem Benchmark, den ich gelaufen bin, scheint es, dass ich nicht so tief gekommen bin.

Hinweis zur Stabilität

Die Wahl hat wirklich nichts mit Stabilitätsproblemen zu tun. Wenn Sie die Tableaus von DifferentialEquations.jl durchgehen (dies ist nur plot(tab)für die Stabilitätsbereiche möglich), werden Sie feststellen, dass die meisten Methoden verdächtig ähnliche Stabilitätsbereiche aufweisen. Dies ist eigentlich eine Wahl. In der Regel führt der Autor beim Ableiten der Methoden die folgenden Aktionen aus:

  1. Ermitteln der niedrigsten Hauptabschneidungsfehlerkoeffizienten (d. H. Der Koeffizienten für die Terme nächster Ordnung)
  2. Vorbehaltlich der Bestellbedingungen
  3. Und machen Sie den Stabilitätsbereich dem der Dormand-Prince-Order-5-Methode nahe.

Warum die letzte Bedingung? Nun, da diese Methode bei der Auswahl der PI-gesteuerten adaptiven Schrittgröße in der Regel immer stabil ist, ist sie ein guter Balken für Stabilitätsbereiche, die "gut genug" sind. Es ist also kein Zufall, dass die Stabilitätsbereiche alle ähnlich sind.

Fazit

Es gibt Kompromisse bei jeder Wahl der Methode. Die RK-Methoden höchster Ordnung sind bei niedrigeren Toleranzen einfach nicht so effizient, weil es schwieriger ist, die Auswahl der Koeffizienten zu optimieren, und weil die Anzahl der Funktionsbewertungen zusammengesetzt ist (und bei Interpolationen sogar noch schneller wächst). Wenn die Toleranz jedoch niedrig genug wird, gewinnt sie, aber die erforderlichen Toleranzen können weit unter "Standard" -Anwendungen liegen (dh wirklich nur für chaotische Systeme anwendbar).

Chris Rackauckas
quelle