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?

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 $A$$B$ (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).

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.

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 A$r\le2$. 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.

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