MN03LAB

From Studia Informatyczne


Własności arytmetyki zmiennoprzecinkowej

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Oglądaj wskazówki i rozwiązania
Ukryj wskazówki i rozwiązania

Ćwiczenie: Równe i równiejsze

Wyjaśnij, dlaczego w arytmetyce podwójnej precyzji IEEE 754 mamy

octave:19> 2006/1e309 ans = 0 octave:20> 2.006/1e306 ans = 2.0060e-306 octave:21> (2006/1000)/(1e309/1000) ans = 0

Oczywiście, "teoretycznie" wszystkie trzy liczby powinny być sobie równe (i niezerowe).

Wskazówka

Kluczem jest wyjaśnienie, jaka jest reprezentacja liczby \displaystyle 10^{309} w podwójnej precyzji... A reprezentacja \displaystyle 10^{306}?

Ćwiczenie: Szeregi zbieżne(?)

Podaj przykłady zbieżnych szeregów postaci \displaystyle \sum_{n=1}^{\infty} a_n, których \displaystyle N-te sumy częściowe obliczone w arytmetyce pojedynczej precyzji algorytmem

suma = 0.0;
for n = 1..N
	suma += <math>\displaystyle a_n</math>;

będą

  • ograniczone niezależnie od \displaystyle N, albo
  • numerycznie rozbieżne, to znaczy takie, że dla bardzo dużych \displaystyle N zachodzi suma == Inf.

Wykonaj to samo zadanie, ale podając przykłady szeregów rozbieżnych (w arytmetyce dokładnej).

Rozwiązanie

Rozpatrz na przykład takie szeregi:

  • Szereg zerowy (numerycznie dostajesz oczywiście też 0)
  • \displaystyle 1 + 10^{2006} - 10^{2006} + 0 + \ldots = 1 oraz \displaystyle 1 + 10^{2006} - 1 + 0 + \ldots = 10^{2006} (numerycznie dostaniesz, odpowiednio, NaN (dlaczego?) i Inf.
  • Szereg harmoniczny (numerycznie sumy częściowe w pewnym momencie przestaną rosnąć).
  • Szereg jedynkowy.

Ćwiczenie

Dla kolejnych \displaystyle N, wyznacz \displaystyle N-tą sumę częściową szeregu Taylora dla funkcji wykładniczej, gdy \displaystyle x=-90:

\displaystyle  e^x \approx \sum_{n=0}^N \frac{x^n}{n!},

korzystając z algorytmu podanego w poprzednim zadaniu. Oblicz błąd względny i bezwzględny wyznaczonego przybliżenia, w porównaniu do wartości wyznaczonej z wykorzystaniem funkcji bibliotecznej exp(). Powtórz następnie dla \displaystyle x=10.

Czy --- zgodnie z teorią matematyczną --- sumy te dążą do wartości \displaystyle e^{x}. Objaśnij dokładnie, co się stało.

Rozwiązanie

Zastosowanie naszego algorytmu dla \displaystyle x=-90 daje w wyniku (dla arytmetyki podwójnej precyzji) sumę równą około \displaystyle 10^{30}, gdy oczywiście naprawdę \displaystyle \exp(-90) \approx 0). Problem jest skutkiem zjawiska redukcji cyfr przy odejmowaniu: aby nasza suma była równa zeru, w pewnym momencie musimy odjąć duży składnik od innego podobnego, by w rezultacie dostać coś małego. Spróbuj wychwycić ten moment i wartość dodawanego składnika.

Dla \displaystyle x=10 mamy, dla rosnącego \displaystyle N, całkiem niezły błąd względny ostatecznego

wyniku.

Ćwiczenie

Już wcześniej stwierdziliśmy, że wyznaczanie \displaystyle e \approx (1 + 1/n)^n dla dużego \displaystyle n nie jest dobrym pomysłem. Przeprowadź eksperyment numeryczny potwierdzający to stwierdzenie i objaśnij jego wyniki.

Wskazówka

Kiedy da znać o sobie epsilon maszynowy?

Ćwiczenie

Jak wiadomo, szereg harmoniczny \displaystyle \sum_{n=1}^\infty 1/n jest rozbieżny. Spróbuj przewidzieć, jaki będzie efekt numerycznego wyznaczenia tej sumy w arytmetyce podwójnej precyzji przy użyciu poniższego kodu.

int dlicznik;
	double dsuma, dstarasuma;
	double dskladnik;
	
	dstarasuma = 0.0; dsuma = 1.0; dlicznik = 1;
	while(dstarasuma != dsuma) 
	{
		dskladnik = 1.0/dlicznik;
		dstarasuma = dsuma;
		dsuma += dskladnik;
		dlicznik++;
	}
	printf("Suma = %e. Liczba składników = %d, składnik = %e\n", 
		dsuma, dlicznik-1, dskladnik);
	 

Wskazówka

Policz, ile mniej więcej obrotów pętli potrzeba, by suma przestała rosnąć.

Rozwiązanie

Nasza suma przekroczy wartość 20, więc aby 20 + dskladnik == 20, musi być dskładnik \displaystyle \approx 10^{-15} (lub więcej). Tymczasem zakres liczb typu integer wynosi niewiele ponad \displaystyle 10^9, więc po jakimś czasie dlicznik przyjmie wartości ujemne.


Ćwiczenie: Zadanie o wyznaczaniu odwrotności bez dzielenia metodą Newtona

Należy wyznaczyć przybliżenie \displaystyle \frac{1}{a} stosując metodę Newtona do równania \displaystyle \frac{1}{x} - a = 0. Zaproponuj dobre przybliżenie początkowe \displaystyle x_0 wiedząc, że \displaystyle a jest liczbą maszynową typu double. Ile iteracji wystarczy, by osiągnąć \displaystyle \epsilon-zadowalające przybliżenie?

Rozwiązanie

Jak pamiętamy, \displaystyle a = (1+f)2^p, skąd \displaystyle \frac{1}{a} = \frac{1}{1+f}2^{-p}. Wystarczy więc przybliżyć \displaystyle \frac{1}{1+f}. Ponieważ dla \displaystyle 0 \leq f < 1,

\displaystyle \frac{1}{2} \leq \frac{1}{1+f} \leq 1,

to możemy jako pierwsze przybliżenie wybrać \displaystyle x_0 = \frac{3}{4}2^{-p} (jak wybrać jeszcze lepsze \displaystyle x_0?)

Wtedy mamy następujące oszacowania:

\displaystyle |x_0a - 1| = |\frac{3}{4}2^{-p}\cdot (1+f)2^p - 1| = |\frac{3}{4}(1+f) - 1| \leq \frac{1}{2}.

Ponieważ iteracja Newtona spełnia

\displaystyle |x_{k+1}a - 1| = |x_ka-1|^2 = |x_0a-1|^{2^k},

to, jeśli chcemy, by błąd względny przybliżenia \displaystyle x_{k+1} spełniał \displaystyle \frac{|x_{k+1} - \frac{1}{a}|}{\frac{1}{a}} = |x_{k+1}a - 1| \leq \epsilon = 2^{-53}, musimy wykonać co najwyżej \displaystyle 6 iteracji. Ponieważ koszt jednej iteracji to 3 działania arytmetyczne, wyznaczenie odwrotności \displaystyle a byłoby 18 razy droższe od mnożenia i dodawania, które mają podobny koszt.

Tymczasem, we współczesnych realizacjach, dzielenie jest tylko około 4 razy droższe od mnożenia! Możesz sprawdzić, jak to się dzieje np. w procesorach AMD.

Więcej na temat implementacji algorytmów dzielenia oraz wyciągania pierwiastka

kwadratowego możesz dowiedzieć się np. z artykułu Petera Marksteina Software Division and Square Root Using Goldschmidt s Algorithms oraz raportów technicznych Stanford Computer Architecture and Arithmetic Group

Ćwiczenie

Niech \displaystyle 0 < a_1 < a_2 < \cdots < a_n. Czy z punktu widzenia błędów w \displaystyle fl_\nu lepiej jest policzyć sumę tych liczb w kolejności od najmniejszej do największej, czy odwrotnie?

Rozwiązanie

Od najmniejszej do największej.

Ćwiczenie

Dlaczego w algorytmie bisekcji rozwiązywania równania \displaystyle f(x)=0, sprawdzając warunek różnych znaków \displaystyle f w krańcach przedziału \displaystyle [a,b], używamy kryterium

if ( sign(f(x)) != sign(f(xl)) )
...	

a nie, matematycznie równoważnego, wyrażenia

if ( f(x)*f(lewy) < 0 )	
...

Wskazówka

Nie daj się zwieść, nie chodzi tu o zaoszczędzenie jednego mnożenia...

Rozwiązanie

Gdy zbliżamy się do miejsca zerowego \displaystyle f, wartości f(x) i f(lewy) mogą być na tyle małe, że ich iloczyn --- poprzez niedomiar --- będzie miał numeryczną wartość zero i w efekcie program może przeskoczyć do niewłaściwej gałęzi pętli if, a tym samym --- wybrać niewłaściwą lokalizację przedziału zawierającego poszukiwany pierwiastek.