Poszukuję Runge-Kutta 8. rzędu w C / C ++

10

Chciałbym użyć metody Runge-Kutta 8. rzędu (89) w niebiańskiej aplikacji do mechaniki / astrodynamiki, napisanej w C ++, za pomocą komputera z systemem Windows. Dlatego zastanawiam się, czy ktoś zna dobrą bibliotekę / implementację, która jest udokumentowana i może być używana bezpłatnie? Jest w porządku, jeśli jest napisany w C, o ile nie ma żadnych problemów z kompilacją.

Do tej pory znalazłem tę bibliotekę (mymathlib) . Kod wydaje się w porządku, ale nie znalazłem żadnych informacji na temat licencjonowania.

Czy możesz mi pomóc, ujawniając niektóre z alternatyw, które możesz znać i które pasowałyby do mojego problemu?

EDYCJA:
Widzę, że tak naprawdę nie ma tak wielu kodów źródłowych C / C ++, jak się spodziewałem. Dlatego też wersja Matlab / Octave byłaby w porządku (nadal musi być darmowa).

c++ reference-request libraries c runge-kutta James C.
źródło

Odpowiedzi:

8

Zarówno Biblioteka Naukowa GNU (GSL) (C), jak i Boost Odeint (C ++) zawierają metody Runge-Kutta 8. rzędu.

Oba są opensource i pod Linuksem i Mac powinny być bezpośrednio dostępne w menedżerze pakietów. Pod oknami prawdopodobnie łatwiej będzie użyć wzmocnienia niż GSL.

GSL jest publikowany na licencji GPL, a Boost Odeint na licencji Boost.

Edycja: Ok, Boost Odeint NIE ma metody Runge-Kutta 89, tylko 78, ale zapewnia przepis na tworzenie dowolnych stepperów Runge-Kutta.

Metody 8. rzędu są jednak dość wysokie i najprawdopodobniej przesadzają z twoim problemem.

Prince-Dormand odnosi się do określonego rodzaju Runge-Kutta i nie jest bezpośrednio związany z zamówieniem, ale najczęściej występuje w 45. Matlabs ode45, który jest ich zalecanym algorytmem ODE, to implementacja Prince-Dormand 45. Jest to ten sam algorytm, który został zaimplementowany w Boost Odeint Runge_Kutta_Dopri5 .

LKlevin
źródło
1
Dziękuję za odpowiedź. OK, co jest teraz zawstydzające, spojrzałem na Boost Odeint, zanim jeszcze zapytałem, i znalazłem tylko „runge_kutta_fehlberg78”. Czy to właściwa rzecz? Właściwie nie znam różnic między metondami, gdy są używane w praktyce, ale szukałem RK89 (zwanego także Dormand-Prince, gdy przeszukuję Internet). Czy możesz skomentować lub rozszerzyć swoją odpowiedź na ten temat? Dziękuję Ci.
James C
Zaktualizowany post, aby odpowiedzieć na twoje pytania. Prince-Dormand 45 najprawdopodobniej ładnie rozwiąże twoje problemy.
LKlevin
15

Jeśli wykonujesz mechanizmy niebieskie w długich skalach czasowych, użycie klasycznego integratora Runge-Kutta nie oszczędza energii. W takim przypadku użycie integratora symplektycznego byłoby prawdopodobnie lepsze. Boost.odeint implementuje również symetryczny schemat Runge-Kutta czwartego rzędu, który działałby lepiej w długich odstępach czasu. O ile mi wiadomo, GSL nie wdraża żadnych metod symplektycznych.

Geoff Oxberry
źródło
Dziękuję za odpowiedź. Czy symetryczna Runge-Kutta czwartego rzędu dawałaby lepsze wyniki niż RKF78, gdyby była używana z satelitami Ziemi (niska orbita, a także głębsza orbita kosmiczna), być może w okresie 1-3 orbit?
James C,
@JamesC Tak. W długim okresie metoda symplektyczna jest znacznie lepsza.
eccstartup
@eccstartup - Co byś tutaj uważał za długi okres? Ponieważ może to być jedna orbita planety wokół Słońca lub kilka orbit satelity pogodowego wokół Ziemi itp.
James C
@JamesC Nie zauważyłem tego dużego problemu. Ale dla moich problemów z modelem, przy obliczonych wielu orbitach, metody symplektyczne dają bardzo doskonałe orbity.
eccstartup
Dlatego zaleca się, aby zaprogramować własną wersję niejawnej metody Runge-Kutty, która zawiera wiele metod symplektycznych o tak wysokim porządku, jak chcesz.
eccstartup
4

podsumowując niektóre punkty:

  1. Jeśli jest to długoterminowa integracja modelu nieulegającego dezaprobacie, to integrator symplektyczny jest tym, czego szukasz.
  2. W przeciwnym razie, ponieważ jest to równanie ruchu, metody Runge-Kutta Nystrom będą bardziej wydajne niż transformacja do układu pierwszego rzędu. Istnieją metody RKN wysokiego rzędu ze względu na DP. Istnieje kilka implementacji, na przykład tutaj w Julii, są one udokumentowane i tutaj jest MATLAB .
  3. Metody Runge-Kutta wysokiej klasy są potrzebne tylko wtedy, gdy potrzebujesz rozwiązania o wysokiej dokładności. Jeśli ma niższe tolerancje, RK 5. rzędu będzie prawdopodobnie szybszy (dla tego samego błędu). Najlepszym rozwiązaniem, jeśli często trzeba to rozwiązać, jest przetestowanie wielu różnych metod. W tym zestawie testów porównawczych dotyczących problemów z trzema ciałami widzimy, że (dla tego samego błędu) metody RK wyższego rzędu są jedynie marginalną poprawą prędkości, chociaż jako błąd -> 0 widać, że poprawa już idzie do> 5x w porównaniu z Dormandem - Książę 45 ( DP5), gdy patrzysz na 4 cyfry dokładności (tolerancje są jednak znacznie niższe w tym przypadku. Tolerancje są tylko problemem dla każdego problemu). Gdy obniżysz tolerancje, jeszcze niższa poprawa z metody wysokiego rzędu RK, ale może być konieczne rozpoczęcie korzystania z liczb o wyższej precyzji.
  4. Algorytm 7/8 kolejności Dormanda-Prince'a ma inny układ tabel ósmego rzędu niż metoda DP853 metod Hairer's dop853i DifferentialEquations.jl DP8(które są takie same). Ta ostatnia metoda 853 nie może zostać zaimplementowana w standardowej wersji tableau metody Runge-Kutta, ponieważ jej estymator błędów jest niestandardowy. Ale ta metoda jest znacznie bardziej wydajna i nie poleciłbym nawet stosowania starszych metod Fehlberga 7/8 lub DP 7/8.
  5. W przypadku metod RK wysokiego rzędu metody „wydajne” Vernera są złotym standardem. Widać to w testach, które połączyłem. Możesz sam je zakodować w Boost lub użyć jednego z 2 pakietów, które je implementują, jeśli chcesz, aby były łatwiejsze (Mathematica lub DifferentialEquations.jl).
Chris Rackauckas
źródło
2

Chciałbym dodać, że choć to, co sugeruje Geoff Oxberry dla długoterminowej integracji (przy użyciu integratorów symplektycznych), jest prawdą, w niektórych przypadkach nie zadziała. Mówiąc dokładniej, jeśli masz siły rozpraszające, twój system nie zachowuje już energii, a zatem nie możesz w tym przypadku skorzystać z integratorów symplektycznych. Osoba zadająca pytanie mówiła o niskich orbitach Ziemi, a takie orbity wykazują dużą siłę oporu atmosferycznego, czyli siłę rozpraszającą, która wyklucza stosowanie takich integratorów symplektycznych.

W tym konkretnym przypadku (oraz w przypadkach, w których nie możesz używać / nie masz dostępu do / nie chcesz używać integratorów symplektycznych), zaleciłbym użycie integratora Bulirsch-Stoer, jeśli potrzebujesz precyzji i wydajności w dłuższych terminach. Działa dobrze z doświadczenia i jest również zalecany przez Przepisy numeryczne (Press i in., 2007).

viiv
źródło
Nie, nie polecaj przepisów numerycznych. Szczególnie w większości przypadków nie zaleca się Burlirsch-Stoer. Jest to dobrze znany problem z książką. Zobacz kilka zmian od najlepszych badaczy w tej dziedzinie tutaj: uwyo.edu/buerkle/misc/wnotnr.html . Jeśli chcesz znaleźć punkty odniesienia w tej kwestii, zobacz pierwszą książkę Hairer, w której zobaczysz, że BS prawie nigdy nie radzi sobie dobrze. Wyższe zamówienie jest bardziej wydajne tylko wtedy, gdy błędy są wystarczająco niskie, a my (i inni) przeprowadziliśmy testy porównawcze, aby dość konsekwentnie wykazać, że jest on skuteczny tylko w przypadku precyzji sub-zmiennoprzecinkowej.
Chris Rackauckas,
Nie mogę za dużo mówić za NR, ponieważ używałem go głównie do ODE, ale wydaje mi się, że skargi na stronie, do której prowadzi link, są stare i zostały rozpatrzone przez autorów NR w odpowiedzi (koniec strony), ale to nie jest temat. Jeśli chodzi o długoterminową integrację orbit z dużą precyzją (powiedzmy, 13-14 cyfr), o czym wspomniałem w mojej odpowiedzi, udowodniono, że od dawna dobrze działają metody ekstrapolacji (patrz rozdział Montenbruck & Gill o integracji numerycznej). Korzystają z niego także nowsze gazety, które udowodniły mi i innym niezawodną i wydajną metodę.
viiv,
M&G testuje go tylko pod kątem dop853, a bardziej nowoczesne metody RK wyższego rzędu, takie jak Verner, są znacznie bardziej wydajne. Wydaje się, że M&G mierzy tylko za pomocą ocen funkcji, które są słabym wskaźnikiem czasu. Nie ma też czasu w stosunku do metod Runge-Kutta Nystrom, które są specjalnie dla ODE drugiego rzędu i są bardziej wydajne niż metody RK pierwszego rzędu stosowane całkiem sporo. Na 13-14 cyfrach BS jest prawdopodobnie konkurencyjny w większości problemów, ale nie jest to oczywisty wybór i nie widziałem schematu precyzji pracy z najnowszymi metodami, które się z tym nie zgadzają.
Chris Rackauckas,
M&G testują RKN pod kątem RK, a BS i inne pod kątem RKN (strony 123–132 i 151–154) i twierdzą, że są najbardziej wydajnymi metodami RK (nie wliczając Vernera, nawet jeśli go cytują). BS okazało się wydajne na 13-14 cyfrach, co było moim twierdzeniem, widziałem, że testowano go przeciwko dop853, ABM (12), Taylor i standardowemu RK8 i działa dobrze. Muszę przyznać, że nie widziałem go przetestowanego ponownie RKN, ale z tego, co widzę w M&G, nie jest na przykład daleko od FILG11. Naprawdę interesuje mnie RK Vernera i przejrzę twoje linki powyżej. Czy masz artykuł, który testuje je wszystkie, aby je zobaczyć?
viiv,
Wróciłem i ponownie przeprowadziłem wiele testów porównawczych na DiffEqBenchmarks.jl i odexnie wydaje mi się, żeby było dobrze. Tak więc przynajmniej dla ODE pierwszego rzędu i dla tolerancji >=1e-13ekstrapolacja nie wydaje się dobrze i zwykle nie jest nawet bliska. Jest to zgodne z powyższym twierdzeniem.
Chris Rackauckas,