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 $10^{309}$ w podwójnej precyzji... A reprezentacja $10^{306}$?

Ćwiczenie: Szeregi zbieżne(?)

Podaj przykłady zbieżnych szeregów postaci $\sum _{n=1}^{\mathrm{\infty }}{a}_n$, których $N$-te sumy częściowe obliczone w arytmetyce pojedynczej precyzji algorytmem

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

będą

• ograniczone niezależnie od $N$, albo
• numerycznie rozbieżne, to znaczy takie, że dla bardzo dużych $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)
• $1+10^{2006}-10^{2006}+0+\dots =1$ oraz $1+10^{2006}-1+0+\dots =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 $N$, wyznacz $N$-tą sumę częściową szeregu Taylora dla funkcji wykładniczej, gdy $x=-90$:

$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 $x=10$.

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

Rozwiązanie

Zastosowanie naszego algorytmu dla $x=-90$ daje w wyniku (dla arytmetyki podwójnej precyzji) sumę równą około $10^{30}$, gdy oczywiście naprawdę $\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 $x=10$ mamy, dla rosnącego $N$, całkiem niezły błąd względny ostatecznego wyniku.

Ćwiczenie

Już wcześniej stwierdziliśmy, że wyznaczanie $e\approx (1+1/n{)}^n$ dla dużego $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 $\sum _{n=1}^{\mathrm{\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 $\approx 10^{-15}$ (lub więcej). Tymczasem zakres liczb typu integer wynosi niewiele ponad $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 $\frac{1}{a}$ stosując metodę Newtona do równania $\frac{1}{x}-a=0$. Zaproponuj dobre przybliżenie początkowe $x_0$ wiedząc, że $a$ jest liczbą maszynową typu double. Ile iteracji wystarczy, by osiągnąć $ϵ$-zadowalające przybliżenie?

Rozwiązanie

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

$\frac{1}{2}\le \frac{1}{1+f}\le 1$,

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

Wtedy mamy następujące oszacowania:

$|x_{0}a-1|=|\frac{3}{4}{2}^{-p}\cdot (1+f)2^p-1|=|\frac{3}{4}(1+f)-1|\le \frac{1}{2}$

Ponieważ iteracja Newtona spełnia

$|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 $x_{k+1}$ spełniał $\frac{|x_{k+1}-\frac{1}{a}|}{\frac{1}{a}}=|x_{k+1}a-1|\le ϵ=2^{-53}$, musimy wykonać co najwyżej $6$ iteracji. Ponieważ koszt jednej iteracji to 3 działania arytmetyczne, wyznaczenie odwrotności $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 $0<a_1<a_2<\cdots <a_n$. Czy z punktu widzenia błędów w $f{l}_{\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 $f(x)=0$, sprawdzając warunek różnych znaków $f$ w krańcach przedziału $[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 $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.