Dlaczego w arytmetyki zmiennoprzecinkowej niedokładność liczbowa wynika z dodania małego terminu do różnicy dużych pojęć?

13

Czytałem książkę Computer Simulation of Liquids autorstwa Allena i Tildesleya. Począwszy od strony 71, autorzy omawiają różne algorytmy stosowane do integracji równań ruchu Newtona w symulacjach dynamiki molekularnej (MD). Począwszy od strony 78, autorzy omawiają algorytm Verleta, który jest być może algorytmem integracji kanonicznej w MD. Stanowią one:

Być może najczęściej stosowaną metodą całkowania równań ruchu jest ta, którą początkowo przyjęła Verlet (1967) i przypisała Stormerowi (Gear 1971). Ta metoda jest bezpośrednim rozwiązaniem równania drugiego rzędu . Metoda oparta jest na postions , przyśpieszeniach i pozycjach z poprzedniego kroku. Równanie przesunięcia pozycji brzmi następująco: $m_i \ddot{\textbf{r}}_i = \textbf{f}_i$ $\textbf{r}(t)$ $\textbf{a}(t)$ $\textbf{r}(t - \delta t)$

$\begin{matrix} (3.14) & r (t + δ t) = 2 r (t) - r (t - δ t) + δ t^{2} a (t) . \end{matrix}$ $\tag{3.14}\textbf{r}(t + \delta t) = 2\textbf{r}(t) - \textbf{r}(t - \delta t) + \delta t^2 \textbf{a}(t).$
Należy zwrócić uwagę na kilka punktów dotyczących eqn (3.14). Okaże się, że prędkości w ogóle się nie pojawiają. Zostały one wyeliminowane przez dodanie równań uzyskanych przez rozwinięcie Taylora o : $\textbf{r}(t)$

$r (t + δ t) = r (t) + δ t v (t) + (1 / 2) δ t^{2} a (t) + . . .$ $\textbf{r}(t + \delta t) = \textbf{r}(t) + \delta t \textbf{v}(t) + (1/2) \delta t^2 \textbf{a}(t) + ...$
$\begin{matrix} (3.15) & r (t - δ t) = r (t) - δ t v (t) + (1 / 2) δ t^{2} a (t) - . . . . \end{matrix}$ $\tag{3.15}\textbf{r}(t - \delta t) = \textbf{r}(t) - \delta t \textbf{v}(t) + (1/2) \delta t^2 \textbf{a}(t) - ... .$

Następnie, później (na stronie 80), autorzy stwierdzają:

W przeciwieństwie do algorytmu Verleta ... forma algorytmu może niepotrzebnie wprowadzać niedokładność liczbową. Wynika to z faktu, że w eqn (3.14) mały termin ( ) jest dodawany do różnicy dużych terminów ( ), w celu wygenerowania trajektorii. $\mathcal{O}(\delta t^2)$ $\mathcal{O}(\delta t^0)$

Myślę, że „małym terminem” jest , a „różnica dużych terminów” wynosi . $\delta t^2 \textbf{a}(t)$ $2\textbf{r}(t) - \textbf{r}(t - \delta t)$

Moje pytanie brzmi: dlaczego niedokładność liczbowa wynika z dodania małego terminu do różnicy dużych terminów?

Interesuje mnie raczej podstawowy, konceptualny powód, ponieważ w ogóle nie znam szczegółów arytmetyki zmiennoprzecinkowej. Czy znasz też jakieś odniesienia typu „przegląd” (książki, artykuły lub strony internetowe), które wprowadziłyby mnie w podstawowe idee arytmetyki zmiennoprzecinkowej związane z tym pytaniem? Dziękuję za Twój czas.

floating-point computer-arithmetic Andrzej
źródło

9

Ich obserwacja „forma algorytmu może niepotrzebnie wprowadzać pewną niedokładność liczbową” jest poprawna. Ale ich wyjaśnienie „Powstaje, ponieważ w równaniu (3.14), mały termin ( ) jest dodawany do różnicy dużych terminów ( ), aby wygenerować trajektorię. „” jest fałszywe. $O(δt^2)$ $O(δt^0)$

Prawdziwym powodem nieznacznego niestabilności numerycznej algorytmu Verlet jest to, że jest tylko nieznacznie stabilny, ponieważ równanie różnica (przede wszystkim w przypadku, gdy zaniedbania w Verlet) ma pasożytnicze rozwiązanie proporcjonalne do , które powoduje, że wprowadzone błędy rosną liniowo w podczas gdy dla w pełni stabilnej metody wieloetapowej stosowanej do rozpraszającego równania różniczkowego wzrost błędu jest ograniczony. $x_{k+1}=2x_k-x_{k-1}$ $a$ $k$ $k$

Edycja: Należy zauważyć, że książka o numerycznej symulacji dynamiki cząsteczkowej i uzyskania wystarczającej dokładności otrzymanych oczekiwania potrzebna jest ogromna liczba etapów zgodnie skali dokładności z tylko . (Często krok czasu jest w pikosekundach, aby podążać za wewnętrzną skalą oscylacji. Ale istotne biologicznie skale czasu są w milisekundach lub większych ( ), chociaż zwykle nie oblicza się tak daleko.) $N$ $O(N^{-1/2})$ $N\sim 10^9$

Aby uzyskać więcej informacji, zobacz http://en.wikipedia.org/wiki/Linear_multistep_method#Stability_and_convergence

Arnold Neumaier
źródło

10

Jeśli szukasz dobrego wprowadzenia, sugerowałbym Davidowi Goldbergowi to, co każdy informatyk powinien wiedzieć o arytmetyki zmiennoprzecinkowej . Może to być trochę zbyt szczegółowe, ale jest dostępne online za darmo.

Jeśli masz dobrą bibliotekę, zasugerowałbym Numeryczne obliczenia Michaela Overtona z arytmetyką zmiennoprzecinkową IEEE lub kilka pierwszych rozdziałów Dokładności i stabilności algorytmów numerycznych Nicka Highama .

Allen i Tildesley odnoszą się konkretnie do anulowania numerycznego . Krótka jest to, że jeśli masz, powiedzmy, tylko trzy cyfry i odjąć 100od 101, można uzyskać 1.00(w trzech cyfr). Liczba wygląda na dokładnie trzy cyfry, ale tak naprawdę tylko pierwsza cyfra jest prawdziwa, a końcowe .00są śmieciami. Dlaczego? Cóż, 100i 101są tylko niedokładnymi reprezentacjami, powiedzmy 100.12345i 101.4321, ale można je przechowywać tylko jako trzycyfrowe liczby.

Pedro
źródło

-1: Gdzie jest anulowanie, które przypisujesz formule Verlet? Zazwyczaj

jest małe, co powoduje, że

, bez skutku anulowania. Spróbuj

!

δ t

$\delta t$

r (\t - δ t) \approx r (t)

$r(\t-\delta t)\approx r(t)$

r (t) = 1

$r(t)=1$

Arnold Neumaier

@ArnoldNeumaier: Tak, przykład Allena i Tildesleya wydaje się nie mieć większego sensu, chciałem jedynie podać odniesienie do problemów pojawiających się, gdy „niewielki termin [..] jest dodawany do różnicy dużych terminów”, co jest tym, co OP zapytał, a nie czy w danym przypadku jest to problem.

Pedro

Ale dodanie krótkiego terminu do dużego terminu to po prostu błąd zaokrąglenia, nic niebezpiecznego. Anulowanie następuje, gdy dwa prawie równe duże terminy są odejmowane, aby uzyskać mały termin. Staje się to problemem tylko wtedy, gdy albo odjęte półprodukty są znacznie większe niż końcowy wynik obliczeń, lub gdy mały wynik pośredni, na który ma wpływ anulowanie, jest dzielony przez inny mały element.

Arnold Neumaier,

@ArnoldNeumaier: Ponieważ, jak sądzę, z mojej odpowiedzi jest dość oczywiste, miałem na myśli problem obliczania różnicy, a nie sumy.

Pedro

1

@ArnoldNeumaier: Punkt zajęty, ale mam nadzieję, że rozumiesz, że uważam to za dość małostkowe jak na „-1”.

Pedro

5

$(3.14)$

r (t) = 101

$\mathbf{r}(t) = 101$

r (t - δ t) = 100

$\mathbf{r}(t - \delta t) = 100$

δ t^{2} a (t) = 1.49

$\delta t^2 \mathbf{a}(t) = 1.49$

$(3.14)$

r (t + δ t) = 103.49

$\mathbf{r}(t + \delta t) = 103.49$

ale ponieważ możemy użyć tylko trzech cyfr, wynik zostaje obcięty do

r (t + δ t) = 103

$\mathbf{r}(t + \delta t) = 103$

$\mathbf{a}(t)$ $\mathbf{r}(t + 20\delta t) = 331$ $433.90$

Igor F.
źródło

Ale efekt jest taki duży tylko w przypadku 3-cyfrowej arytmetyki dziesiętnej.

Arnold Neumaier

3

Pedro już podaje ważny fakt, a mianowicie odwołanie. Chodzi o to, że każda liczba, z którą obliczasz, ma powiązaną dokładność; na przykład liczba zmiennoprzecinkowa pojedynczej precyzji może reprezentować tylko rzeczy o dokładności do około 8 cyfr. Jeśli masz dwie liczby, które są prawie dokładnie takie same, ale różnią się siódmą cyfrą, różnica będzie ponownie 8-cyfrową liczbą zmiennoprzecinkową z pojedynczą precyzją i wygląda na to, że jest dokładna do 8 cyfr, ale w rzeczywistości tylko pierwsza 1 lub 2 cyfry są dokładne, ponieważ ilości, z których zostały obliczone, nie są dokładne poza pierwszą 1 lub 2 cyframi różnicy.

Teraz cytowana przez ciebie książka pochodzi z 1989 roku. W tamtych czasach obliczenia były najczęściej wykonywane w pojedynczej precyzji, a zaokrąglanie i anulowanie były poważnymi problemami. Obecnie większość obliczeń jest wykonywana przy użyciu podwójnej precyzji z 16 cyframi dokładności, a dzisiaj jest to o wiele mniejszy problem niż kiedyś. Myślę, że warto przeczytać przytoczone akapity z odrobiną soli i wziąć je w kontekście ich czasu.

Wolfgang Bangerth
źródło

anulowanie w arytmetyki podwójnej precyzji może być równie dużym problemem, jak w przypadku pojedynczej precyzji. Przykładem jest eliminacja Gaussa bez obrotu, co często daje bardzo słabe wyniki z powodu anulowania, nawet z podwójną precyzją.

Arnold Neumaier,

-1: Formuła Verleta zazwyczaj zachowuje wszystkie cyfry dokładności, a nie tylko 1 lub 2 z 8 w pojedynczej precyzji.

Arnold Neumaier

@ArnoldNeumaier: Jasne, że możesz mieć podobne problemy z podwójną precyzją. Powiedziałem tylko, że nie spotyka się ich tak często.

Wolfgang Bangerth,

Jeśli stracisz 6 cyfr trzy razy w ciągu obliczeń, stracisz wszystkie cyfry nawet z podwójną precyzją. Algorytmy cierpiące na anulowanie będą zwykle słabe, nawet z podwójną precyzją. Algorytm Verleta jest inny, ponieważ nie ma anulowania, ale niewielki liniowy wzrost błędów. W związku z tym utrata dokładności nie może się zwielokrotnić, przez co nadaje się do znacznie dłuższych czasów integracji. To z pewnością było znane Allenowi i Tildesleyowi.

Arnold Neumaier

Dobrze. Ale mam na myśli to, że jeśli masz algorytm bez anulowania, nadal ponosisz błąd rzędu 1e-8 w pojedynczej precyzji, a jeśli wykonasz kroki czasowe 1e8, możesz mieć problem, nawet jeśli wszystko inne jest dokładne. Kroki czasowe 1e8 to rząd wielkości, jaki możesz mieć dla ODE. Z drugiej strony, z podwójną precyzją, twoja niedokładność na każdym kroku wynosi 1e-16 i wymagałoby kroków czasowych 1e16, aby uzyskać całkowitą utratę dokładności. Jest to szereg kroków, których nie spotkasz w praktyce.

Wolfgang Bangerth,

Dlaczego w arytmetyki zmiennoprzecinkowej niedokładność liczbowa wynika z dodania małego terminu do różnicy dużych pojęć?

Odpowiedzi: