Dlaczego metody Runge – Kutta wyższego rzędu nie są używane częściej?

17

Byłem ciekawy, dlaczego wysokopoziomowe (tj. Większe niż 4) metody Runge – Kutty prawie nigdy nie są omawiane / stosowane (przynajmniej o ile mi wiadomo). Rozumiem, że wymaga to dłuższego czasu obliczeniowego na krok (np. RK14 z osadzonym krokiem 12. rzędu ), ale czy są jeszcze inne wady stosowania metod Runge – Kutty wyższego rzędu (np. Problemy ze stabilnością)? Czy przy zastosowaniu do równań z wysoce oscylującymi rozwiązaniami w ekstremalnych skalach czasowych takie metody wyższego rzędu nie byłyby zwykle preferowane?

Mathews24
źródło
2
Myślę, że to bardzo subiektywne pytanie. Największym minusem, jak już zauważyłeś, jest koszt obliczeń. Ogólnie staramy się zachować równowagę między dokładnością a czasem obliczeniowym. W PDE, kiedy ludzie mówią o wyższym porządku, zazwyczaj myślą o 3. lub 4. rzędzie. Krok czasowy jest również utrzymywany w tej samej kolejności.
Vikram
3
W PDE schemat dokładności wysokiego rzędu dla zależności czasowej nie ma sensu, jeśli dokładność przestrzenna jest gorsza. W rzeczywistości dokładność zależności przestrzennej wynosi głównie około 2 lub 3 rzędu, szczególnie podczas pracy na nieustrukturyzowanych siatkach. Ludzie muszą kontrolować globalne obcinanie błędów przy jak najmniejszym koszcie, dlatego uważa Runge-Kutta z wystarczająco wysoką dokładnością w określonych przypadkach.
tqviet
@tqviet Jeśli stosujesz przybliżenia różnicowe wstecz lub centralnie do rzędu 8 dla pochodnych przestrzennych, RK8 byłoby odpowiednie, nie? Zasadniczo, czy są jakieś problemy z dokładnością lub stabilnością przy stosowaniu tak dużych przybliżeń różnic skończonych rzędu pochodnych przestrzennych?
Mathews24,
1
@ Mathews24: Nie wspomniałem o stabilności, która silnie zależy od równania. Gdy bardzo dokładne system jest stosowany do zależności przestrzennej przyjąć RK do zależności czasowej z co najmniej tego samego rzędu dokładności, ale stan stabilności może wymagać mniejszej wartości . Δt
tqviet

Odpowiedzi:

17

Istnieją tysiące artykułów i setki kodów przy użyciu metod Runge-Kutta piątego rzędu lub wyższych. Zauważ, że najczęściej używany jawny integrator w MATLAB jest ODE45, który rozwija rozwiązanie za pomocą metody Runge-Kutta 5. rzędu.

Przykłady szeroko stosowanych wysokiej klasy metod Runge-Kutta

Papier Dormand i Książę daje metodę 5-ty rzędu ponad 1700 cytowań według Google Scholar . Większość z nich to papiery wykorzystujące swoją metodę rozwiązania jakiegoś problemu. Artykuł z metody Cash-Karp ma ponad 400 cytowań . Być może najczęściej stosowaną metodą rzędu wyższego niż 5 jest metoda Prince-Dormand rzędu 8 która ma ponad 400 cytowań w Google Scholar . Mógłbym podać wiele innych przykładów; i pamiętaj, że wielu (jeśli nie większość) osób stosujących te metody nigdy nie powołuje się na dokumenty.

Zauważ też, że ekstrapolacja i odroczona korekta wysokiego rzędu są metodami Runge-Kutty .

Metody wysokiego rzędu i błąd zaokrąglania

Jeśli twoja dokładność jest ograniczona błędami zaokrąglania , powinieneś użyć metody wyższego rzędu . Wynika to z faktu, że metody wyższego rzędu wymagają mniejszej liczby kroków (i mniejszej oceny funkcji, nawet jeśli jest więcej ocen na krok), więc popełniają mniej błędów zaokrąglania. Możesz to łatwo zweryfikować samodzielnie za pomocą prostych eksperymentów; jest to dobry problem do zaliczenia pierwszego kursu analizy numerycznej.

Metody dziesiątego rzędu są niezwykle przydatne w arytmetyce podwójnej precyzji. Wręcz przeciwnie, gdybyśmy mieli tylko metodę Eulera, błąd zaokrąglania byłby poważnym problemem i potrzebowalibyśmy bardzo precyzyjnych liczb zmiennoprzecinkowych dla wielu problemów, w których solwery wysokiego rzędu dobrze sobie radzą.

Metody wysokiego rzędu mogą być równie stabilne

@RichardZhang odniósł się do drugiej bariery Dahlquista, ale dotyczy to tylko metod wieloetapowych. Zadane tutaj pytanie dotyczy metod Runge-Kutta, a istnieją metody Runge-Kutta dla każdego zamówienia, które są nie tylko stabilne , ale także stabilne BZAb (właściwość stabilności przydatna dla niektórych nieliniowych problemów). Aby dowiedzieć się o tych metodach, zobacz na przykład tekst Hairer & Wanner.

Metody wysokiego rzędu w mechanice niebieskiej

Ty pytasz

Czy przy zastosowaniu do równań z wysoce oscylującymi rozwiązaniami w ekstremalnych skalach czasowych takie metody wyższego rzędu nie byłyby zwykle preferowane?

Masz rację! Najlepszym tego przykładem jest mechanika niebieska. Nie jestem ekspertem w tej dziedzinie. Ale na przykład ten artykuł porównuje metody mechaniki niebieskiej i nawet nie uważa, że ​​porządek jest niższy niż 5. Stwierdza się, że metody rzędu 11 lub 12 są często najskuteczniejsze (przy metodzie Prince-Dormand rzędu 8 również często bardzo wydajny).

David Ketcheson
źródło
Ketchson: czy mógłbyś podać dowody lub wyjaśnienia dotyczące tego stwierdzenia: „ekstrapolacja wysokiego rzędu i metody odroczonej korekty są metodami Runge-Kutta”? Zwłaszcza „metody odroczonej korekty”. Dzięki.
tqviet
@David Ketcheson Czy możesz omówić, jak zmieniłaby się twoja odpowiedź, jeśli użyjesz zwalidowanych (zweryfikowanych) technik obliczeniowych, takich jak interwał zaokrąglony na zewnątrz lub arytmetyka radialna? Co powiesz na to, że zastosowano arytmetykę zaokrągloną na zewnątrz o podwójnej precyzji na zewnątrz lub arytmetykę promieniową? Co stanie się z zawijaniem i zależnością w miarę zwiększania kolejności Runge-Kutta, i dla zabawy powiedzmy, że ODE jest bardzo sztywny.
Mark L. Stone,
@ MarkL.Stone To zupełnie inny zestaw pytań. Jeśli chcesz je zadać, prześlij je jako osobne pytania. Nie jestem jednak ekspertem w tych sprawach i nie będę w stanie odpowiedzieć.
David Ketcheson
1
@tqviet Spójrz na ten dokument, aby uzyskać wyjaśnienie.
David Ketcheson
12

Tak długo, jak używasz standardowej arytmetyki zmiennoprzecinkowej podwójnej precyzji, metody bardzo wysokiego rzędu nie są potrzebne, aby uzyskać rozwiązanie z wysoką dokładnością w rozsądnej liczbie kroków. W praktyce uważam, że dokładność rozwiązania jest zwykle ograniczona do błędu względnego 1,0e-16 przez reprezentację zmiennoprzecinkową podwójnej precyzji, a nie liczbę / długość kroków, które są podejmowane z RKF45.

Jeśli przejdziesz na schemat arytmetyczny zmiennoprzecinkowy wyższy niż podwójna precyzja, warto skorzystać z metody dziesiątego rzędu.

Brian Borchers
źródło
5
Myślę, że ta odpowiedź jest myląca. Metody wysokiego rzędu prowadzą do znacznie mniejszego błędu zaokrąglania, podczas gdy metody niskiego rzędu cierpią z powodu dominującego błędu zaokrąglania, gdy wymagana dokładność jest duża lub przedział czasu jest długi; patrz moja odpowiedź poniżej.
David Ketcheson
2
Chodzi o to, że w zmiennoprzecinkowym podwójnej precyzji nie można nawet przedstawić rozwiązania z dokładnością względną większą niż 1.0e-16. W wielu praktycznych sytuacjach stary dobry RKF45 doprowadzi Cię do tego poziomu dokładności w okresie, który Cię interesuje, bez konieczności wykonywania drobnych kroków. Może nie być dobrym wyborem dla sztywnych systemów lub sytuacji, w których wymagany jest integrator symplektyczny, ale metoda Runge Kutta wyższego rzędu również nie jest świetnym rozwiązaniem w takich sytuacjach. Zgadzam się, że przez bardzo długi czas metody Runge Kutta wyższego rzędu mogą mieć sens.
Brian Borchers,
10

Aby dodać do doskonałej odpowiedzi Briana Borchera, wiele rzeczywistych aplikacji dopuszcza bardzo sztywne ODE lub DAE. Intuicyjnie problemy te doświadczają nierównomiernych, nagłych zmian w czasie, więc lepiej je modelować przy użyciu wielomianów niskiego rzędu rozmieszczonych drobno na krótkich krokach, w przeciwieństwie do wielomianów wysokiego rzędu rozciągniętych na długie rozmiary stopni. Również stabilność często wymaga użycia metod niejawnych , dla których kara obliczeniowa metod wyższego rzędu jest znacznie bardziej rygorystyczna.

Bardziej rygorystycznie, metody wyższego rzędu są mniej stabilne niż metody niższego rzędu w przypadku sztywnych problemów. Mamy na przykład bariery Dahlquista dla liniowych metod wieloetapowych.

Twierdzenie (druga bariera Dahlquista). Konieczna jest metoda wielostopniowa stabilna na poziomie Ar2). Spośród wszystkich wieloetapowych metod rzędu 2 reguła trapezowa ma najmniejszą stałą błędu.

Podobne (ale o wiele bardziej skomplikowane) stwierdzenia dotyczące stabilności L we wzorach RK. We wszystkich przypadkach wzrost kolejności często nie zawsze prowadzi do dokładniejszych rozwiązań. Poniżej znajduje się fragment przełomowego artykułu Prothero i Robinsona z 1974 r .:

Stosując jednoetapowe metody stabilne A do rozwiązywania dużych układów sztywnych nieliniowych równań różniczkowych stwierdziliśmy, że
(a) niektóre metody stabilne A dają wysoce niestabilne rozwiązania i
(b) dokładność rozwiązań uzyskanych, gdy równania są sztywność często wydaje się być niezwiązana z kolejnością zastosowanej metody.

Jeszcze bardziej rygorystyczne podejście do tego tematu można znaleźć w klasycznym tekście Hairer & Wannera, „Rozwiązywanie zwykłych równań różniczkowych II: Sztywne i różnicowe - problemy algebraiczne”, 1991.

W praktyce równania sztywne prawie zawsze rozwiązuje się za pomocą reguły trapezoidalnej lub wzoru TR-BDF2 (funkcje ode23t i ode23tb w MATLAB). Oba są niejawnymi metodami drugiego rzędu. Oczywiście tam, gdzie stabilność nie jest problemem (tj. W równaniach nietrwałych), mamy swobodę wyboru spośród wielu opcji; RK45 jest najczęstszym wyborem.

Richard Zhang
źródło
Bardzo interesujące. Czy istnieje (intuicyjne) wyjaśnienie, dlaczego kolejność musi być mniejsza lub równa 2, aby była stabilną metodą A wieloetapową? I żeby wyjaśnić, kiedy mówicie, że można sformułować podobne stwierdzenia dla formuł RK, to czy jest to znowu rząd 2?
Mathews24,
Ale w przypadku metod Runge-Kutta istnieją stabilne metody A dowolnego porządku.
David Ketcheson,
@DavidKetcheson Tak, ale nie są silnie stabilne na A (tj. Stabilne na L). Mają wiele problemów przy rozwiązywaniu problemów DAE, np. Symulują proste obwody tranzystorowe. Rzeczywiście, TR jest niesławny z powodu wywoływania sztucznego dzwonienia w SPICE, co motywowało rozwój TR-BDF2.
Richard Zhang,
@DavidKetcheson W celach informacyjnych patrz doi.org/10.1090/S0025-5718-1974-0331793-2 . Pojęcie stabilności A nie jest wystarczająco silne dla DAE, a metody stabilizacji A wysokiego rzędu często dają dziwne wyniki, gdy są stosowane do rozwiązywania DAE.
Richard Zhang,
Jasne, ale pytanie nie dotyczy DAE ani metod wieloetapowych.
David Ketcheson
9

Konfiguracja testu porównawczego

W oprogramowaniu Julia DifferentialEquations.jl zaimplementowaliśmy wiele metod wyższego rzędu, w tym metody Feagina. Możesz to zobaczyć na naszej liście metod , a następnie istnieje mnóstwo innych, których możesz użyć jako dostarczonych tabel . Ponieważ wszystkie te metody są zebrane razem, możesz łatwo porównywać między nimi. Możesz zobaczyć testy, które mam online tutaj i że bardzo łatwo jest przetestować wiele różnych algorytmów. Jeśli chcesz poświęcić kilka minut na przetestowanie testów, wybierz go. Oto podsumowanie tego, co wychodzi.

Po pierwsze, należy zauważyć, że jeśli spojrzysz na każdy z testów, zobaczysz, że nasze DP5(zamówienie Dormand-Prince Prince 5) i DP8metody są szybsze niż kody Hairer Fortran ( dopri5i dop853), a zatem te implementacje są bardzo dobrze zoptymalizowane . Pokazują one, że, jak zauważono w innym wątku, nadmierne użycie metod Dormanda-Prince'a wynika z tego, że metody te są już napisane, a nie dlatego, że wciąż są najlepsze. Tak więc rzeczywiste porównanie najbardziej zoptymalizowanych implementacji jest między metodami Tsitorous, Verner i Feagin z DifferentialEquations.jl.

Wyniki

Ogólnie, metody rzędu wyższego niż 7 mają dodatkowy koszt obliczeniowy, który zwykle nie jest równoważony przez porządek, biorąc pod uwagę wybrane tolerancje. Jednym z powodów jest to, że wybory współczynników dla metod niższego rzędu są bardziej zoptymalizowane (mają małe „zasadnicze współczynniki błędu obcięcia”, które mają większe znaczenie, gdy nie jesteś asymetrycznie mały). Widać, że w wielu problemach, takich jak tutaj, metody Verner Efficient 6 i 7 działają wyjątkowo dobrze, ale metody takie jak Verner Efficient 8 mogą mieć niższe nachylenie. Wynika to z tego, że „korzyści” wyższego rzędu łączą się przy niższych tolerancjach, więc zawsze istnieje tolerancja, w której metody wyższego rzędu będą bardziej wydajne.

Pytanie brzmi jednak, jak nisko? W dobrze zoptymalizowanej implementacji poziom ten jest dość niski z dwóch powodów. Pierwszym powodem jest to, że metody niższego rzędu implementują coś o nazwie FSAL (pierwszy taki sam jak ostatni). Ta właściwość oznacza, że ​​metody niższego rzędu ponownie wykorzystują ocenę funkcji z poprzedniego kroku w następnym kroku, a tym samym mają efektywnie jedną ocenę mniejszą funkcji. Jeśli zostanie to właściwie zastosowane, wówczas coś w rodzaju metody 5. rzędu (Tsitorous lub Dormand-Prince) faktycznie bierze 5 ocen funkcji zamiast 6, które sugerowałyby tableau. Dotyczy to również metody Verner 6.

Drugi powód wynika z interpolacji. Jednym z powodów korzystania z metody bardzo wysokiego rzędu jest podejmowanie mniejszej liczby kroków i po prostu interpolowanie wartości pośrednich. Jednak w celu uzyskania wartości pośrednich funkcja interpolująca może wymagać większej liczby ocen funkcji niż użyto do wykonania kroku.Jeśli spojrzysz na metody Vernera, potrzeba 8 dodatkowych ocen funkcji dla metody Order 8, aby uzyskać interpolant rzędu 8. Wiele razy metody niskiego rzędu zapewniają „wolny” interpolant, na przykład większość metod piątego rzędu ma swobodną interpolację czwartego rzędu (bez dodatkowych ocen funkcji). Oznacza to, że jeśli potrzebujesz wartości pośrednich (które będą potrzebne dla dobrej fabuły, jeśli używasz metody wysokiego rzędu), istnieją dodatkowe ukryte koszty. Uwzględnij fakt, że te interpolowane wartości są naprawdę ważne w obsłudze zdarzeń i rozwiązywaniu równań różniczkowych opóźnienia, i rozumiesz, dlaczego wpływają na to dodatkowe koszty interpolacji.

A co z metodami Feagina?

Zobaczysz więc, że podejrzanie brakuje metod Feagina w testach porównawczych. Są w porządku, testy zbieżności działają na liczbach o dowolnej dokładności itp., Ale aby je dobrze wykonać, musisz poprosić o kilka absurdalnie niskich tolerancji. Na przykład w niepublikowanych testach porównawczych stwierdziłem, że Feagin14osiąga lepsze wyniki Vern9(metoda Vernera 9 rzędu) przy tolerancjach takich jak 1e-30. W przypadku aplikacji z chaotyczną dynamiką (jak w przypadku problemów Pleides lub astrofizyki 3-ciał) możesz chcieć takiej dokładności ze względu na wrażliwą zależność (błędy w układach chaotycznych szybko się łączą). Jednak większość ludzi prawdopodobnie wykonuje obliczenia na liczbach zmiennoprzecinkowych o podwójnej precyzji i nie znalazłem testu porównawczego, w którym osiągają lepsze wyniki w tej dziedzinie tolerancji.

Ponadto nie ma interpolanta zgodnego z metodami Feagina. Więc po prostu umieszczam na nich interpolację Hermite'a trzeciego rzędu, aby taki istniał (i działa zaskakująco dobrze). Jeśli jednak nie ma standardowej funkcji interpolacji, możesz wykonać rekurencyjną metodę Hermite (użyj tej interpolacji, aby uzyskać punkt środkowy, a następnie interpolacji 5. rzędu itp.), Aby uzyskać interpolację wysokiego rzędu, ale jest to bardzo kosztowne, a wynikowe interpolacja niekoniecznie ma niską zasadę błędu skracania (więc jest dobra, gdy dtjest naprawdę mała, co jest dokładnym przeciwieństwem pożądanego przypadku!). Więc jeśli kiedykolwiek potrzebujesz naprawdę dobrej interpolacji w celu dopasowania do swojej dokładności, musisz przynajmniej powrócić do czegoś takiego Vern9.

Uwaga na temat ekstrapolacji

Zauważ, że metody ekstrapolacji są po prostu algorytmami do generowania metod Runge-Kutta o dowolnym porządku. Jednak dla swojej kolejności podejmują więcej kroków niż to konieczne i mają wysokie współczynniki błędów skracania, a zatem nie są tak wydajne, jak dobrze zoptymalizowana metoda RK przy danym zamówieniu. Ale biorąc pod uwagę poprzednią analizę, oznacza to, że istnieje dziedzina o wyjątkowo niskiej tolerancji, w której metody te będą lepsze niż „znane” metody RK. Ale w każdym teście, który przeprowadziłem, wydaje mi się, że nie byłem tak niski.

Uwaga na temat stabilności

Wybór naprawdę nie ma nic wspólnego z kwestiami stabilności. W rzeczywistości, jeśli przejdziesz przez tabelę DifferentialEquations.jl (możesz tylko plot(tab)dla regionów stabilności) zobaczysz, że większość metod ma podejrzanie podobne regiony stabilności. To jest właściwie wybór. Zwykle podczas opracowywania metod autor zwykle wykonuje następujące czynności:

  1. Znajdź najniższe współczynniki błędu obcięcia zasady (czyli współczynniki dla warunków następnego zamówienia)
  2. Z zastrzeżeniem ograniczeń zamówienia
  3. I uczyń obszar stabilności zbliżonym do tego z metody Dormand-Prince Order 5.

Dlaczego ostatni warunek? Cóż, ponieważ metoda ta jest zawsze stabilna w sposobie dokonywania wyborów adaptacyjnych stopniowania kontrolowanych przez PI, więc jest to dobry słupek dla „wystarczająco dobrych” obszarów stabilności. To nie przypadek, że wszystkie regiony stabilności są zwykle podobne.

Wniosek

W każdym wyborze metody występują kompromisy. Metody RK najwyższego rzędu po prostu nie są tak wydajne przy niższych tolerancjach, zarówno dlatego, że trudniej jest zoptymalizować wybór współczynników, jak i dlatego, że liczba złożonych funkcji oceny związków (i rośnie nawet szybciej, gdy zaangażowane są interpolacje). Jednakże, jeśli tolerancja stanie się wystarczająco niska, wygrywają, ale wymagane tolerancje mogą być znacznie poniżej „standardowych” aplikacji (tj. Naprawdę dotyczą tylko systemów chaotycznych).

Chris Rackauckas
źródło