MN03LAB

Własności arytmetyki zmiennoprzecinkowej

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

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

Ć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_{0}a-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_{k}a-1|^{2}=|x_{0}a-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.