błędy na floatach i równania różniczkowe

Suma floatów, a raczej błąd, zależy od kolejności sumowania:
(a + b) + c będzie zwykle różne od a + (b+c).

W przypadku sumowania wielu składników błąd rośnie statystycznie zgodnie z sqrt(n), n - liczba składników; ale w najgorszym scenariuszu błąd może być dowolnie duży, nawet ponad 100% poprawnej sumy.

Przykład:
suma 1/i = 1 + 1/2 + 1/3 + ... + 1/N;

Sumujemy w pojedynczej precyzji, czyli float.

float s = 0;
for(i = 1; i <= N; i++) s +=  1.0f/i;

Dla porównania sumujemy tak samo ale na double:

Wyniki:
N = 1 mln; s = 14.357358, błąd 0.25%
N = 10 mln; s = 15.403683, błąd 8%
N = 100 mln; s = 15.403683, błąd 20%

Najmniejszy błąd byłby dla odwrotnej kolejności sumowania: 1/N + ... + 1/2 + 1;
ale suma nadal będzie mniejsza od tej z podwójną precyzją i bez odwracania.

Na double są oczywiście również błędy, ale musimy sumować więcej składników - biliony, czy nawet tryliony.

Niby dużo, bo po co sumować tyle liczb?
Zwykle nie ma takiej potrzeby, no, ale przecież w przypadku równań różniczkowych my właśnie bez przerwy sumujemy:
http://pl.wikipedia.org/wiki/Algorytm_Rungego-Kutty

y_k+1 = y_k + dy; dy - wyliczamy z równań...

W innych metodach podobnie sumujemy, tylko inaczej wyliczamy te dy.

**Pytanie:
Mamy tu taką katastrofalną kumulację błędów, jak w przykładzie, czy też nie?
**

I jeśli nie, to dlaczego?

Weźmy przykładowo numeryczną symulację Układu Słonecznego.
W astronomii często symulują to nawet w miliardach lat - sprawdzają sobie co się stanie z orbitami, np. czy Księżyc nam ucieknie, albo czy orbita Ziemi się powiększy, itp.

Potem rozgłaszają rewelacyjne odkrycia, np. Merkury wyleci z orbity za 6 mln lat, albo

A ostatnio, i coraz częściej, nawet weryfikują takimi symulacjami teorie - modele fizyczne, np. ta słynna precesja Merkurego: zapuszczamy symulację n-body na kilka tysięcy lat, liczymy wszystko wg równań Newtona, i sprawdzamy czy tam faktycznie brakuje 43'' / wiek.

Potem korygujemy te 43'', np. zwiększając troszeczkę dy, być może dokładnie tyle, ile nam ucieka w tym sumowaniu... i już gra!

Można oczywiście sumować lepiej, ale to dość dużo kosztuje.

sama natura ma z tym problem :)
właśnie tak powstaje chaos

Poczytaj jak komputer zapisuje liczby zmiennoprzecinkowe np na wiki, potem zastanów sie co sie dzieje, gdy dodajesz liczby różniące się bardzo wykładnikiem.

MarekR22 napisał(a)

Poczytaj jak komputer zapisuje liczby zmiennoprzecinkowe np na wiki, potem zastanów sie co sie dzieje, gdy dodajesz liczby różniące się bardzo wykładnikiem.

Wiadomo co się dzieje ale pytanie było inne:
jaki jest faktycznie błąd numerycznych metod dla równań różniczkowych?

W szczególności dla problemu Keplera, w którym mamy sumowanie dookoła - średnia zero, a nie takie narastające.
Gdy tu dodamy trzecie ciało, wtedy one wzajemnie się zaburzają, i wtedy właśnie te orbity się obracają.

Błędy obliczeń też są przecież zaburzeniami, zatem również obracają orbitami, co nawet dobrze widać: gdy obliczamy ze zbyt dużym krokiem, wtedy orbita kręci się do tyłu.

Tradycyjnie mówią że metoda rzędu n ma błąd (h/T)^{n, no ale to jest przecież błąd dla obliczeń z nieskończoną precyzją; liczony z Taylora: f(x) = a + bt + ct}2 + ... urywamy to na k-tym wyrazie i mamy błąd rzędu t^k.

Drugi błąd jest z zaokrągleń podczas obliczeń, proporcjonalny chyba do: (T/h)^0.5 * eps.
Ale to jest przecież statystyczny błąd, a nie minimalny, a tylko taki może nam dać jakieś szanse oceny wiarygodności rozwiązania.

W całkowaniu równań różniczkowych błąd numeryczny zwykle nie jest tak istotny jak lokalny błąd obcięcia (wynikający ze skończoności kroku). Pocieszające jest to, że oba te błędy lokalne dają się całkiem ładnie oszacować. Natomiast to, czy błędy będą się kumulować, zależy przede wszystkim od samego problemu - jeśli problem nie jest stabilny, albo jest blisko granicy stabilności (np. układ oscylujący o dużej dobroci), albo używamy nie A-stabilnego schematu różnicowego (np. Runge-Kutty; wszystkie >2 rzędu nie są A-stabilne) to tak, mogą się kumulować i może być to poważny problem. Żeby to sprawdzić, wystarczy wykonać obliczenia z inną precyzją / innym krokiem i sprawdzić, jaką to powoduje różnicę w wyniku.

Podstawowe błędy podczas symulacji:
*zły algorytm (algorytm Eulera ma błąd systematyczny, który dla symulacji planet prowadzi do systematycznych błędów, układ nie utrzymuje stałe energii)
*obliczania na zbyt mało precyzyjnych liczbach
*zbyt MAŁY czas kroku całkowania! Gdy krok czasowy jest za mały dodajesz małe poprawki do dużej liczby (np współrzędnej) co prowadzi do dużych błędów numerycznych
*zbyt duży czas kroku czasowego, wtedy przybliżenie jest za mało dokładne

Zwróć uwagę, że krok czasowy może być równocześnie zbyt mały i zbyt duży! Przykład układ planet, w którym w pewnym momencie dystans między dwoma ciałami jest bardzo mały (np przelatujący meteoryt).

Krolik napisał(a)

W całkowaniu równań różniczkowych błąd numeryczny zwykle nie jest tak istotny jak lokalny błąd obcięcia (wynikający ze skończoności kroku). Pocieszające jest to, że oba te błędy lokalne dają się całkiem ładnie oszacować. Natomiast to, czy błędy będą się kumulować, zależy przede wszystkim od samego problemu - jeśli problem nie jest stabilny, albo jest blisko granicy stabilności (np. układ oscylujący o dużej dobroci), albo używamy nie A-stabilnego schematu różnicowego (np. Runge-Kutty; wszystkie >2 rzędu nie są A-stabilne) to tak, mogą się kumulować i może być to poważny problem. Żeby to sprawdzić, wystarczy wykonać obliczenia z inną precyzją / innym krokiem i sprawdzić, jaką to powoduje różnicę w wyniku.

Przecież ta stabilność niczego tu nie gwarantuje.
Energia zachowana, kręt itd.

Z samych zaokrągleń będzie jakiś błąd +/-dE, który niech się zeruje po czasie.

Ale przecież równania ruchu nie są raczej liniowe względem energii, ani krętu, więc błąd położenia nie będzie tak samo symetrycznie oscylował, a tym samym średnia błędu położenia (wyliczanych współrzędnych) się nie wyzeruje.

Niemożliwe?
Można się bez problemu przemieścić na dowolną odległość z zerowym wydatkiem energii:

• rozpędzam się kosztem energii dE
• lecę swobodnie przez dowolny czas t
• hamuję i tym samym odzyskuję dE.

Bilans: -dE + dE = 0, nie ma zmian energii, a przemieszczenie wynosi: d = vt = sqrt(2E/m) t ;

http://pl.wikipedia.org/wiki/IEEE_754

Jeżeli korzystasz z dev'a to przeżuć się na coś innego, np. Code::Blocks

program który liczy średnią liczb dodatnich z przekazanej tablicy liczb całkowitych
Przykładowy problem w Dev'ie i to nie w przypadku pętli N=100 000

MarekR22 napisał(a)

Podstawowe błędy podczas symulacji:
*zły algorytm (algorytm Eulera ma błąd systematyczny, który dla symulacji planet prowadzi do systematycznych błędów, układ nie utrzymuje stałe energii)
*obliczania na zbyt mało precyzyjnych liczbach
*zbyt MAŁY czas kroku całkowania! Gdy krok czasowy jest za mały dodajesz małe poprawki do dużej liczby (np współrzędnej) co prowadzi do dużych błędów numerycznych
*zbyt duży czas kroku czasowego, wtedy przybliżenie jest za mało dokładne

Zwróć uwagę, że krok czasowy może być równocześnie zbyt mały i zbyt duży! Przykład układ planet, w którym w pewnym momencie dystans między dwoma ciałami jest bardzo mały (np przelatujący meteoryt).

Te szczegół z konfliktem długości kroków może być istotny w n-body.
Jest na to jakiś sposób?

Chyba należałoby liczyć każdą planetę z innym krokiem - proporcjonalnym do okresu, czyli Merkury miałby krótszy z 50 razy od Jowisza.
No, ale do wyliczania sił potrzebne są aktualne położenia wszystkich planet.

Może wyliczamy z krótszym krokiem, ale tylko doraźnie:
wyliczamy pozycję wszystkich z mniejszym krokiem, z czego można już wyliczyć co trzeba dla ciała z mniejszym krokiem. Ale nie przesuwamy trwale tych z większymi.

y_k = f(t); to trzymamy

potem potrzebujemy coś pośredniego:

y_k+ = f(t+h'); h' < h;

i odrzucamy to od razu po użyciu

y_k+1 = f(t+h); nowa pozycja liczona z y_k, a nie z tych pomiędzy

Ciekawe jaki byłby efekt takich kombinacji - na szybkość nie powinno mieć to wpływu, a błąd zaokrągleń byłby chyba mniejszy.