całkowanie numeryczne równań Newtona

Jest jakaś dobra metoda do obliczania równań różniczkowych, powiedzmy rzędu 16 ?

Chodzi o równania Newtona, czyli: y'' = f(r),
np. r'' = -GM/r^{2 * r}0^ dla symulacji orbit planet.

rząd równania różniczkowego, oznaczaj jak krotność pochodnych występuje w równaniu. Równania newtona dla planet, mają tylko drugą pochodną po czasie:

Nie rozumiem skąd ci się wzięło 16.

Do symulacji orbit wystarczą ci algorytmy: Verleta albo LeapFrog.

Chodzi mi o rząd metody 16.

Ta szkolna metoda Runge-Kutty: RK4 ma rząd tylko 4, czyli jest strasznie słaba - mało dokładna oraz czasochłonna.
http://pl.wikipedia.org/wiki/Algorytm_Rungego-Kutty

Gdybym całkował tym, powiedzmy orbitę Księżyca, wówczas należałoby stosować
bardzo krótki krok czasowy w obliczeniach, dla uzyskania dokładności nawet na poziomie 1e-8.

Ale krok metody nie może być dowolnie mały, ponieważ wówczas liczba obliczeń byłaby duża,
a błędy się kumulują, więc dokładność znowu spada dla zbyt krótkiego kroku.

Dla tej metody, przy obliczaniu orbity Księżyca, minimum wynosi około: dt = 10 minut = 600s,
co daje liczbę kroków na orbitę: n = 27 dni / 10 min = 4000.

I powiedzmy że chcemy symulować orbitę w okresie 100 lat, ponad 1300 orbit, co daje ok. 5 milionów kroków =~ 1e7.

Obliczamy z precyzją 1e-15, zatem błąd samych zaokrągleń będzie tu na poziomie: 1e7 * 1e-15 = 1e-8;
i przyjmując minimalny błąd 1e-15 na krok, ale sama metoda ma na pewno ma większy błąd na krok,
ponieważ tam jest minimum 10 operacji na liczbach double.

Zatem minimalny błąd będzie tu pewnie rzędu 1e-6.

Weźmy teraz cały Układ Słoneczny - kilkanaście ciał, powiedzmy 10^2 = 100 obliczeń samych sił,
czyli otrzymamy pewnie 100-1000 razy większy błąd!

Zatem, w symulacji kilkunastu ciał, metoda 4-go rzędu daje nam dokładność rzędu 1e-3 = 1/1000 zaledwie!

Metoda rzędu 8 miałaby błąd pewnie rzędu 1e-6, ponieważ krok byłby znacznie większy - w godzinach.
I takim sposobem dochodzimy, że błąd na poziomie 1e-15, czyli double, uzyskamy dopiero dla metody rzędu 16;
wtedy obliczamy z krokiem w dniach, czyli zarazem bardzo szybko - można nawet tysiące lat symulować.

Rząd metody 16 to by była tragedia. Rząd 4 to już jest zwykle dużo i działa to dobrze tylko w specyficznych zastosowaniach. Najbardziej uniwersalne są metody rzędu 2 i 3. LeapFrog (rzędu 2) jest IMHO bardzo dobry i prosty.

To nie jest tak, że im wyższy rząd metody, tym wyższa dokładność. Rząd metody mówi tylko o tym, że jeśli wynikiem całkowania jest wielomian stopnia < niż rząd metody całkowania, to wynik nie będzie obarczony błędami obcięcia. W przypadku, gdy tak nie jest (na ogół nie jest), to niestety nie zawsze interpolacja wielomianem wyższego stopnia jest bezwzględnie dokładniejsza niż wielomianem niższego stopnia.

Sprawa jest trochę bardziej złożona - wszystko zależy od kroku całkowania. Rząd metody dość dobrze określa ile razy maleje błąd obcięcia wraz ze zmniejszeniem kroku całkowania (dokładniej x-krotne zmniejszenie kroku powoduje spadek lokalnego błędu obcięcia x^(k+1) krotny, gdzie k rząd metody). Z tego powodu dla odpowiednio dużych kroków całkowania metody wyższych rzędów są znacznie mniej dokładne niż metody rzędów niższych. Gorzej, im wyższy rząd metody, tym krok musi być mniejszy, aby metoda miała przewagę dokładnościową nad metodą niższego rzędu. A im mniejszy krok i większy rząd metody, tym dłużej wszystko trwa, a do tego nie można kroku zmniejszać za bardzo, bo pojawia się problem błędów zaokrągleń (jak sam napisałeś).

Do tego dochodzi jeszcze jeden problem: stabilność. Wszystkie metody rzędów wyższych niż 2 nie są A-stabilne. Niektóre metody rzędu 3 (BDF) są "prawie" A-stabilne. Brak A-stabilności powoduje, że jeśli krok całkowania nie będzie mniejszy niż stała czasowa symulowanego układu, to metoda nie tłumi błędów obcięcia i w wyniku można dostać kompletne śmieci - np. rzeczywisty układ nie oscyluje, a wynik symulacji oscyluje i amplituda tych oscylacji rośnie (!). Oczywiście im wyższy rząd metody, tym łatwiej o niestabilność i tym łatwiej o otrzymanie śmieci zamiast poprawnych wyników.

Jeżeli chcesz zwiększyć dokładność, to lepiej zastosować metody całkowania ze zmiennym krokiem i ze zmiennym rzędem niż na stałe zmniejszyć krok i podnieść rząd. Są wzory na estymację lokalnego błędu obcięcia i na ich podstawie można sterować krokiem i rzędem.

Krolik napisał(a):

To nie jest tak, że im wyższy rząd metody, tym wyższa dokładność. Rząd metody mówi tylko o tym, że jeśli wynikiem całkowania jest wielomian stopnia < niż rząd metody całkowania, to wynik nie będzie obarczony błędami obcięcia. W przypadku, gdy tak nie jest (na ogół nie jest), to niestety nie zawsze interpolacja wielomianem wyższego stopnia jest bezwzględnie dokładniejsza niż wielomianem niższego stopnia.

LeapFrog i inne niskiego rzędu, są wolniejsze od RK4, ponieważ mają niższy rząd właśnie.

I to jest oczywiste: błąd jest proporcjonalny do h^n;
niech h = T/10000, T - okres orbity;
czyli mamy 10000 kroków na orbitę, a błąd samej metody wynosi:
1/10000^n, i widać że to dla n = 2 jest milion razy większe niż dla n= 4, a nie mniejsze, co tu sugerujesz.

Zatem możemy zwiększyć krok metody n = 4 zachowując błąd zgodny z metodą n=2:
1/10000^{2 = 1/k}4 => k = 100

czyli metodzie rzędu 4 wystarczy krok T/100 żeby uzyskać błąd na poziomie metody rzędu 2 z krokiem T/10000 !

Tylko 100 kroki zamiast 10 tyś., tym samym i błędy zaokrągleń będą znacznie mniejsze, a jednocześnie będzie znacznie szybciej,
ale nie aż 100 razy, lecz z 20 powiedzmy, ponieważ metoda rzędu 4 ma więcej obliczeń na krok.

20-to krotny wzrost szybkości i przy tym mniejszy błąd - przecież to czysty i gigantyczny zysk!

Dla metody 8 będzie jeszcze lepiej: ona jest tak samo lepsza od 4 jak 4 od 2.

I okazuje się że optymalnym krokiem dla liczb double jest coś w okolicach n = 16.
Wyższy rząd byłby oczywiście jeszcze lepszy, ale powyżej n=16 double mają już za małą precyzję, niestety.

W astronomii nie przypadkiem obliczają na metodach rzędu co najmniej n = 8 do 15 i wielu lat!

Krolik napisał(a):

Jeżeli chcesz zwiększyć dokładność, to lepiej zastosować metody całkowania ze zmiennym krokiem i ze zmiennym rzędem niż na stałe zmniejszyć krok i podnieść rząd. Są wzory na estymację lokalnego błędu obcięcia i na ich podstawie można sterować krokiem i rzędem.

To mogę sobie zawsze zrobić - dla metody rzędu 16 też.
Ale są znane analizy, w których wykazano że to jest szkodliwe dla długiego
czasu symulacji (przynajmniej w symulacjach tego typu, tj. grawitacja i n-body).

Sprawdziłem kilka metod i okazuje się że metody niskich rzędów nadają się tylko do prawie okrągłych orbit.

Np. dla orbity z mimośrodem e > 0.5 metoda LeapFrog kręci całą orbitą dookoła.
A dla a = 0.9 zupełnie się sypie - nawet dla bardzo małego kroku.

Rząd 8 trzyma się dobrze nawet dla e = 0.99, czyli odległość zmienia się aż od 1 do 100, i prędkość podobnie.

b = a(1-e²⁾0.5, co dla a = 1 i e = 0.99 daje: b/a = 1/7
e = 0.999 => b/a = 1/22, czyli elipsa rozciągnięta 22 razy - prawie linia prosta,
gdzie odległość zmienia się aż od 1 do 1000.

Jeżeli symulujesz bardzo "mocno eliptyczne" orbity to IMHO największy zysk dokładności / wydajności będziesz miał z zastosowania SR zmiennokrokowego.

I to jest oczywiste: błąd jest proporcjonalny do h^n;

Przy stałym n, zmiennym h i założeniu że h jest takie, że SR jest stabilny.
Zmiana n powoduje zwykle zmianę współczynnika proporcjonalności.

Możliwe, że układy n-ciał są ogólnie łatwiejsze w symulacji niż układy elektroniczne i być może tam wyższe rzędy mają uzasadnienie. Natomiast w elektronice, gdzie równania różniczkowe są bardzo często sztywne (tj. mają równocześnie ekstremalnie małe jak i ekstremalnie duże stałe czasowe), SR wyższych rzędów niż 3 są stosowane rzadko.

Z takich ciekawostek - symulacja nietłumionego układu LC za pomocą SR Geara 6 rzędu da Ci oscylacje gasnące, natomiast ten sam układ symulowany SR trapezowym lub SR Geara 23 (średnia arytmetyczna z Geara 2 i Geara 3, technicznie 2-rzędu, ale o lepszej dokładności) daje wynik poprawny i zachowuje ładunek / energię.

Podsumowując rząd metody to nie wszystko.

Krolik napisał(a):

Jeżeli symulujesz bardzo "mocno eliptyczne" orbity to IMHO największy zysk dokładności / wydajności będziesz miał z zastosowania SR zmiennokrokowego.

Oczywiście, ale im wyższy rząd tym mniej kroków wykonasz dla tego samego błędu obliczeń - o wiele mniej.

Krolik napisał(a):

Możliwe, że układy n-ciał są ogólnie łatwiejsze w symulacji niż układy elektroniczne i być może tam wyższe rzędy mają uzasadnienie. Natomiast w elektronice, gdzie równania różniczkowe są bardzo często sztywne (tj. mają równocześnie ekstremalnie małe jak i ekstremalnie duże stałe czasowe), SR wyższych rzędów niż 3 są stosowane rzadko.

Nie sądzę. Proste oscylatory nie mają osobliwości, która niestety jest w równaniach z siłami typu 1/r^2,
i właśnie dlatego orbity mocno eliptyczne są bardzo podatne na błędy metod - im wyższy rząd tym lepiej.

Krolik napisał(a):

Z takich ciekawostek - symulacja nietłumionego układu LC za pomocą SR Geara 6 rzędu da Ci oscylacje gasnące, natomiast ten sam układ symulowany SR trapezowym lub SR Geara 23 (średnia arytmetyczna z Geara 2 i Geara 3, technicznie 2-rzędu, ale o lepszej dokładności) daje wynik poprawny i zachowuje ładunek / energię.

Widocznie to są marne metody.
Najprostsza metoda symplektyczna zachowa energię, np. ta Cromera, czyli 1-go rzędu. Verlet i Frog też są symplektyczne.

Ale symplektyczne znowu fazy przesuwają, więc to też nie jest super, np. orbity się obracają dookoła wolno,
i praktycznie niezależnie od kroku, nawet dla h -> 0 jest obrót (charakterystyczny dla danej metody).

Dlatego należy stosować dokładne metody, czyli maksymalnie wysokiego rzędu... dla danej precyzji, w której obliczamy, zwykle double.

Sprawdzałem kilka metod i okazuje się że jednak rząd 16 jest za wysoki dla double.
Już rząd 8 ma błędy na poziomie 1e-16, i kroku h = 1 doba, a większy jest już raczej niewskazany
dla symulacji Układu Słonecznego - Merkury ma okres 90 dni zaledwie.

Zatem używanie rzędu 16 miałoby chyba sens dopiero w symulacji z poczwórną precyzją i bardzo długich - w milionach lat.
Wówczas ostateczny błąd byłby rzędu 1e-20 powiedzmy, a przypadku double byłoby tu pewnie minimum 1e-8, czyli już dość spore błędy.

Oczywiście, ale im wyższy rząd tym mniej kroków wykonasz dla tego samego błędu obliczeń - o wiele mniej.

Ale ja trochę nie bardzo widzę, skąd Ci się wzięło, że n > m => hⁿ < h^m, przy m > 1 (pomińmy różne stałe na razie). Weź odpowiednio duże h i będzie dokładnie na odwrót.

Poza tym wzór na błąd bezwzględny nie jest hⁿ - to jest bardzo duże uproszczenie. Jest tylko proporcjonalny do hⁿ, tzn. zmniejszając krok 2 razy zmniejszasz błąd globalny 2ⁿ krotnie (i w drugą stronę). W rzeczywistości tam są stałe, czyli prawdziwy wzór wygląda bardziej tak: błąd = c1 * (c2 * h)ⁿ, gdzie c1 i c2 zależą od rodzaju SR, jego rzędu oraz rodzaju rozwiązywanego problemu i stałych czasowych układu. Oznacza to, że im mniejszy krok tym zwykle opłaca się podnieść rząd metody. Im większy krok, tym rząd powinien być niższy. Natomiast na temat wyboru punktu, gdzie opłaca się przełączyć i jak sterować rzędem/krokiem to niejedną publikację napisano, więc nie jest to takie proste jak piszesz "im wyższy rząd tym szybciej/dokładniej".

Widocznie to są marne metody

Są to metody dostosowane do rozwiązywania układów sztywnych, tj. takich w których zjawiska zachodzą w różnych skalach czasowych równocześnie. Charakteryzują się bardzo dużą stabilnością.

Dlatego należy stosować dokładne metody, czyli maksymalnie wysokiego rzędu... dla danej precyzji, w której obliczamy, zwykle double

Wysoki rząd implikuje bardzo mały krok symulacji. Inaczej dokładność będzie gorsza niż dla metod rzędów niższych. Rząd dobiera się na podstawie wielkości kroku, a wielkość kroku na podstawie akceptowalnego błędu obcięcia oraz dopuszczalnych błędów zaokrągleń. Błąd obcięcia maleje wraz ze zmniejszaniem kroku, natomiast błędy zaokrągleń rosną i trzeba znaleźć jakiś kompromis.

I nie ma czegoś takiego jak "błąd double". Błędy zaokrągleń są zależne nie tylko od precyzji użytych liczb, ale też od tego co z nimi robisz i jak to zakodujesz (np. kolejność operacji ma znaczenie). Zwykle im wyższy rząd, tym więcej obliczeń na krok, a zarazem większy błąd zaokrągleń na krok.

<quote="975592">

Ale ja trochę nie bardzo widzę, skąd Ci się wzięło, że n > m => hⁿ < h^m, przy m > 1 (pomińmy różne stałe na razie). Weź odpowiednio duże h i będzie dokładnie na odwrót.

To się bierze wprost z teorii szeregów.
Obliczamy funkcję f(t+h), mając dane f(t), czyli warunki początkowe: y0, i t0;
w kolejnych rolę początkowych pełnią poprzednio wyliczone, czyli krokach yi, ti = t0+h*i,

f(t+h) = f(t) + f'(t)h + ...

i teraz jeśli ta funkcja istnieje w obliczanym punkcie,
czyli jest określona i skończona, wówczas to musi być szereg zbieżny, co znaczy że: f⁽ⁿ⁾ * hⁿ -> 0, z n -> oo.

Dla funkcji okresowych, oscylatory, orbity, itp. mamy tu funkcję typu: f(t) = A(t)sin(wt),
i zwykle A(t) < A_max, tz. energia jest ujemna - amplituda nie rośnie nam do nieskończoności.

Zatem obliczamy tu w zasadzie sinus:
sin(w(t+h)) = sin(wt + wh) = sin(wt) + cos(wt)w(hw) - sin(wt)w²(hw)²/2 + ...

w = 2pi/T
czyli tu będzie wyraz: (h/T)ⁿ, a nie samo hⁿ, i dlatego to szybciej nam zbiega im większe n, czyli rząd metody,
pod warunkiem h < T/2pi, gdzie T - jakiś tam średni okres... albo raczej minimalny.
w = df/dt, czyli gdy to jest małe wówczas krok h może być duży, a gdy jest duże wówczas trzeba zmniejszać.

Przypadek silnie eliptycznych orbit: tu w tym punkcie minimalnej odległości 'w' jest potężne - szybki obrót o prawie 180 stopni,
i dlatego algorytmy się na tym sypią.

<quote="975592">Poza tym wzór na błąd bezwzględny nie jest hⁿ - to jest bardzo duże uproszczenie. Jest tylko proporcjonalny do hⁿ, tzn. zmniejszając krok 2 razy zmniejszasz błąd globalny 2ⁿ krotnie (i w drugą stronę). W rzeczywistości tam są stałe, czyli prawdziwy wzór wygląda bardziej tak: błąd = c1 * (c2 * h)ⁿ, gdzie c1 i c2 zależą od rodzaju SR, jego rzędu oraz rodzaju rozwiązywanego problemu i stałych czasowych układu. Oznacza to, że im mniejszy krok tym zwykle opłaca się podnieść rząd metody. Im większy krok, tym rząd powinien być niższy. Natomiast na temat wyboru punktu, gdzie opłaca się przełączyć i jak sterować rzędem/krokiem to niejedną publikację napisano, więc nie jest to takie proste jak piszesz "im wyższy rząd tym szybciej/dokładniej".<quote>

Może być. Wtedy to c2 jest ograniczone tak: c2 < 2pi/T_min = w_max, oraz c2*hc < 1 - wtedy jest to zbieżne i tym szybciej, im większy rząd metody.

Krolik napisał(a):

Wysoki rząd implikuje bardzo mały krok symulacji. Inaczej dokładność będzie gorsza niż dla metod rzędów niższych. Rząd dobiera się na podstawie wielkości kroku, a wielkość kroku na podstawie akceptowalnego błędu obcięcia oraz dopuszczalnych błędów zaokrągleń. Błąd obcięcia maleje wraz ze zmniejszaniem kroku, natomiast błędy zaokrągleń rosną i trzeba znaleźć jakiś kompromis.

I nie ma czegoś takiego jak "błąd double". Błędy zaokrągleń są zależne nie tylko od precyzji użytych liczb, ale też od tego co z nimi robisz i jak to zakodujesz (np. kolejność operacji ma znaczenie). Zwykle im wyższy rząd, tym więcej obliczeń na krok, a zarazem większy błąd zaokrągleń na krok.

Nie wiem skąd te przekonania o mniejszych krokach dla wyższego rzędu.
Wszędzie używano by metod nawet rzędu 100, gdyby tylko takie istniały - były opracowane... no gdyby komputery obliczały na 1000 bitach zamiast 53.

Niskie rzędy są dobre do zgrubnych obliczeń - szacunkowych, ewentualnie z lenistwa, bo one są bardzo proste, więc łatwo to zaprogramować.

Ten szereg ma być trochę:

sin(w(t+h)) = sin(wt + wh) = sin(a + x), i po prostu rozwijamy w punkcie a, czyli otrzymamy:

sina + cosa * x - sina * x²/2 - cosa * x³/6 ...

tak powinno być.

Natomiast te przypadki sztywne można też odnieść do orbit:
powiedzmy że interesuje nas orbita Jowisza, więc wydaje się że można całkować
Układ Słoneczny z krokiem nawet h = 1 tydzień, ponieważ Jowisz ma okres aż 4333 dni.

No, ale Merkury ma krótszy okres: 88 dni, więc tak duży krok może doprowadzić tu wypadku.

Wszędzie używano by metod nawet rzędu 100, gdyby tylko takie istniały

Istnieją. Można zrobić metodę dowolnego rzędu.

no gdyby komputery obliczały na 1000 bitach zamiast 53.

Nie widzę problemu. Można liczyć na dowolnej liczbie bitów.

Powód, dla którego się nie używa metod tak wysokich rzędów jest banalny - im wyższy rząd, tym:

1. większy koszt obliczeniowy pojedynczego kroku
2. większy błąd zaokrągleń dla pojedynczego kroku
3. mniejszy obszar stabilności metody (wymagany mniejszy krok)
4. mniejszy krok przy którym metoda ma przewagę dokładnościową nad metodą rzędu niższego (znowu wymagany mniejszy krok)

Punkty 3. i 4. wynikają wprost z tego, że SR bazują na aproksymacji poszukiwanej funkcji inną funkcją o znanej postaci, którą można całkować analitycznie, np. aproksymacji wielomianowej. Np. aproksymując funkcją stałą (wielomian 0-stopnia) otrzymasz SR Eulera, który jest 1-rzędu. Aproksymując funkcją liniową rozpiętą na f(_n) i f(_n+1) dostajesz SR trapezowy, który jest 2-rzędu. Aproksymacja budowana jest zawsze na podstawie wartości funkcji z kroków poprzednich (SR otwarty) lub na podstawie wartości z kroków poprzednich oraz wartości znajdowanej w bieżącym kroku (SR zamknięty). Zwykle jest niestety tak, że im wyższego rzędu jest aproksymacja, tym więcej danych historycznych jest branych pod uwagę i błąd aproksymacji zaczyna bardzo silnie rosnąć ze wzrostem h.

Poczytaj sobie o efekcie Rungego: http://pl.wikipedia.org/wiki/Efekt_Rungego

Krolik napisał(a):

Istnieją. Można zrobić metodę dowolnego rzędu.

Mi właśnie o to chodzi, zatem podaj przykładowe metody wysokiego rzędu, np. 16 i 32.

I oczywiście żeby tam nie było przesadnie dużo obliczeń, bo takie coś też dyskwalifikuje metodę.

Dla metod rzędu n=4 są zwykle 4 obliczenia wartości f(x,y);
ta standardowa RK4 ma chyba 5, a są też z 3-ma - i to jest chyba minimum;
dla n=2 może być tylko 1, np. LeapFrog, a dla n = 1 też jest 1, bo mniej nie wypada.

Zatem ile ma, musi być wywołań funkcji w metodzie rzędu n=8, i większych: 16, 32?

W przypadku równań Newtona są to obliczenia sił, a dla układu planetarnego
złożonego z N ciał będzie N(N-1)/2 ~ N^2 par sił, które trzeba wyliczać.

Krolik napisał(a):

Nie widzę problemu. Można liczyć na dowolnej liczbie bitów.

Jeśli szybkość obliczeń spadnie milion razy, a pewnie tak będzie dla 1000 bitów precyzji,
no to już w zasadzie po symulacji - nie ma sensu obliczać 5 lat stanu systemu po jednym roku.

Nakład obliczeń wzrasta chyba wykładniczo z precyzją: 2 razy większa precyzja, 2² = 4 razy dłużej obliczamy.
1000/53 = ~20, a 2²⁰ to akurat milion, czyli trafiłem idealnie.

Krolik napisał(a):

Powód, dla którego się nie używa metod tak wysokich rzędów jest banalny - im wyższy rząd, tym:

1. większy koszt obliczeniowy pojedynczego kroku
2. większy błąd zaokrągleń dla pojedynczego kroku
3. mniejszy obszar stabilności metody (wymagany mniejszy krok)
4. mniejszy krok przy którym metoda ma przewagę dokładnościową nad metodą rzędu niższego (znowu wymagany mniejszy krok)

1. większy koszt na krok, ale krok będzie z 1000 razy większy.
   To należy przeliczać na czas - ten rzeczywisty, i tu będzie mniejszy nakład dla większego rzędu.

Mamy krok h oraz obliczamy k kroków na s, czyli szybkość metody:
s = k * h [s/s] - sekund obliczeń na sekundę rzeczywistą - tego symulowanego systemu.

I teraz dla n = 1: krok h = 1s, oraz k = 1 mln kroków /s:
s = 1 milion s/s,

a dla metody rzędu 4 krok będzie minimum h = 1000s, a k - kroków / s tylko z 5 razy mniejsze, czyli wyjdzie:
s = 1000*1mln/5 = 200 mln s/s - 200 razy szybciej.

2. mniej obliczeń = szybciej = mniej błędów zaokrągleń
3. stosujemy adekwatne metody do problemu, a nie przypadkowe
4. niższy rząd ma większy błąd (przypadków anomalnych nie liczy się - p-t 3)

Krolik napisał(a):

Punkty 3. i 4. wynikają wprost z tego, że SR bazują na aproksymacji poszukiwanej funkcji inną funkcją o znanej postaci, którą można całkować analitycznie, np. aproksymacji wielomianowej. Np. aproksymując funkcją stałą (wielomian 0-stopnia) otrzymasz SR Eulera, który jest 1-rzędu. Aproksymując funkcją liniową rozpiętą na f(_n) i f(_n+1) dostajesz SR trapezowy, który jest 2-rzędu. Aproksymacja budowana jest zawsze na podstawie wartości funkcji z kroków poprzednich (SR otwarty) lub na podstawie wartości z kroków poprzednich oraz wartości znajdowanej w bieżącym kroku (SR zamknięty). Zwykle jest niestety tak, że im wyższego rzędu jest aproksymacja, tym więcej danych historycznych jest branych pod uwagę i błąd aproksymacji zaczyna bardzo silnie rosnąć ze wzrostem h.

Poczytaj sobie o efekcie Rungego: http://pl.wikipedia.org/wiki/Efekt_Rungego

To dotyczy tylko interpolacji dla równoodległych węzłów.

Nic się nie psuje w metodach wyższego rzędu - krok nie musi być zmniejszany, lecz ma być zwiększany.

A te przypadki sztywne to po prostu próby obliczania z h >> T, czyli np. symulacja orbity Księżyca z krokiem 1 rok.
Niekompetencja badacza, czy programisty, a nie wada metody.

Mi właśnie o to chodzi, zatem podaj przykładowe metody wysokiego rzędu, np. 16 i 32.

BDF można wygenerować dowolnego rzędu. Problem jest tylko w tym, że BDF rzędu > 2 nie są A-stabilne a BDF rzędu > 6 nie są 0-stabilne i podejrzewam, że problem ten dotyczy też metod RK wyższych rzędów. Czyli w sumie dokładnie to o czym pisałem wcześniej - im wyższy rząd, tym gorsza stabilność i przydatność w praktyce maleje.

http://en.wikipedia.org/wiki/Backward_differentiation_formula

Cytując za wikipedią o metodach RK:

The Gauss–Legendre method with s stages has order 2s, so its stability function is the Padé approximant with m = n = s. It follows that the method is A-stable.[20] This shows that A-stable Runge–Kutta can have arbitrarily high order. In contrast, the order of A-stable linear multistep methods cannot exceed two.

Tu masz RK 10-rzędu: http://sce.uhcl.edu/feagin/courses/rk10.pdf
W artykule jest nawet ładny wykres przedstawiający w jakich sytuacjach metoda RK niższego rzędu jest efektywniejsza (= dokładniejsza przy tym samym nakładzie obliczeń), a w jakich nie. Ten RK 10. rzędu jest raczej pomyślany bardziej o liczbach poczwórnej precyzji.

Dlaczego nie możesz stosować RK 4. rzędu po prostu z mniejszym krokiem? Zauważ, że tylko 10-krotne zmniejszenie kroku powoduje aż 10000-krotne zmniejszenie globalnego błędu (a 100000-krotne lokalnego), a nakład na pojedynczy krok nie jest wcale duży. Dlatego RK rzędu 4 są tak bardzo popularne. W elektronice powszechnie stosuje się SR trapezowy oraz BDF rzędu 2 lub 3 - które mają w końcu niższy rząd niż RK 4, a mimo to w zupełności wystarczają do uzyskania 3-5 miejsc znaczących.

Jeśli szybkość obliczeń spadnie milion razy, a pewnie tak będzie dla 1000 bitów precyzji,
no to już w zasadzie po symulacji - nie ma sensu obliczać 5 lat stanu systemu po jednym roku.

Nakład obliczeń wzrasta chyba wykładniczo z precyzją: 2 razy większa precyzja, 22 = 4 razy dłużej obliczamy.
1000/53 = ~20, a 220 to akurat milion, czyli trafiłem idealnie.

Dawno większej bzdury nie widziałem. Skąd ją wziąłeś?
Gwoli wyjaśnienia: złożoność dodawania i odejmowania jest liniowa względem liczby bitów, złożoność mnożenia jest O(n^1.465) lub w przypadku dłuższych liczb O(n log n log log n). Do tego obecne procesory potrafią w jednym takcie wykonywać kilka operacji na 64-bitowych liczbach, więc obliczenia na liczbach np. 256-bitowych wcale takie powolne być nie muszą. Owszem, będą powolniejsze niż na wbudowanych doublach, ale nie milion razy, a co najwyżej kilka (jak się postarasz) lub kilkadziesiąt (jak olejesz optymalizacje).

Krolik napisał(a):

Dlaczego nie możesz stosować RK 4. rzędu po prostu z mniejszym krokiem? Zauważ, że tylko 10-krotne zmniejszenie kroku powoduje aż 10000-krotne zmniejszenie globalnego błędu (a 100000-krotne lokalnego), a nakład na pojedynczy krok nie jest wcale duży. Dlatego RK rzędu 4 są tak bardzo popularne. W elektronice powszechnie stosuje się SR trapezowy oraz BDF rzędu 2 lub 3 - które mają w końcu niższy rząd niż RK 4, a mimo to w zupełności wystarczają do uzyskania 3-5 miejsc znaczących.

W symulacji Układu Słonecznego potrzeba aż 12 cyfr, przynajmniej, co dla RK4 i długotrwałej symulacji nie przejdzie.
3 cyfry to błąd 1/1000, więc np. dla orbity Ziemi błąd położenia byłby 150 mln / 1000 = 0.15 miliona km.

A gdybyś chciał symulować powiedzmy lot na Jowisza miałbyś z 0.5 miliona km błędu - byłaby to praktycznie losowa wyprawa...
Słońce ma promień 0.7 mln km, a Jowisz jest z 10 razy mniejszy od Słońca, i przy okazji 10 większy od Ziemi, czyli około 0.13 mln km średnicy.

Poza tym tu dopiero dla setek tysięcy lat symulacji widać coś ciekawego:
zmiany inklinacji orbit, obroty peryheliów, kształt orbit - eliptyczność, itp.
Dlatego to musi być perfekt - pełne 16 cyfr, na milion lat, czyli musimy obliczać z 30 cyfr, bo sam błąd zaokrągleń rośnie tu stale.

To nie jest oblizanie równania typu:
y' = ay, gdzie sobie wyliczasz kawałeczek krzywej w zadanym i króciutkim przedziale.