MN01LAB

Eksperymenty ze środowiskiem obliczeń numerycznych

W Linuxie czas działania programu można zbadać poleceniem time.

Ćwiczenie

Który z poniższych programów wykona się szybciej?

x = 1.0;
for( i = 0; i < N; i++)
	x = x/3.0; 

czy

x = 1.0; f = 1.0/3.0;
for( i = 0; i < N; i++)
	x = x*f; 

Rozwiązanie

Oczywiście, szybszy będzie program nie wykorzystujący dzielenia. Optymalizujący kompilator (gcc -O3) strawi, a nawet będzie jeszcze bardziej zadowolony z pozornie rozrzutnego kodu

x = 1.0; 
for( i = 0; i < N; i++)
	x = x*(1.0/3.0); 

dlatego, że stałą, przez którą trzeba mnożyć $x$, wyliczy przed wykonaniem programu.

Sprawdź, czy z wyłączoną optymalizacją ten kod okaże się najwolniejszy ze wszystkich... (okazuje się, że nie!)

Ćwiczenie

Napisz program w C, który zapisuje do pliku

• tekstowego
• binarnego

kolejne wartości $\sin (\pi \cdot i\cdot 0.4)$, gdzie $i=0,\dots ,1024$. Następnie porównaj rozmiary plików i możliwości ich odczytania zewnętrznymi narzędziami. Wreszcie, wczytaj liczby z pliku i porównaj je z oryginalnymi wartościami sinusa. Czy możesz wyjaśnić przyczyny różnic?

Powtórz to samo w Octave.

Rozwiązanie

Różnice są skutkiem konwersji liczb zmiennoprzecinkowych do formatu dziesiętnego. Oczywiście, zapis w formacie binarnym daje dokładną kopię zawartości pamięci, nie ma więc żadnych strat.

Ćwiczenie: Implementacja funkcji matematycznych

Pomyśl, jak obliczać, korzystając jedynie z czterech działań podstawowych: $+,-,\times ,÷$, wartość funkcji exp($x$) = $e^x$ dla dowolnych $x$ rzeczywistych. Naszym kryterium jest, by $|e^x-\exp (x)|\le ϵ$, czyli by błąd bezwzględny aproksymacji nie przekroczył zadanego $ϵ$.

Wykonaj eksperymenty w C lub w Octave, pokazujące koszt metody w zależności od $x$ oraz w zależności od $ϵ$. Przeprowadź też sekwencję testów potwierdzających twoje rachunki co do oczekiwanej dokładności (porównując się z funkcją biblioteczną). W C możesz korzystać ze stałej M_E $\approx e=\exp (1)$, zdefiniowanej w pliku nagłówkowym math.h.

Rozwiązanie

W pierwszym odruchu myślimy o szeregu definiującym funkcję wykładniczą,

$e^x=\sum _{n=0}^{\mathrm{\infty }}\frac{x^n}{n!}$

lub o granicy

$e^x=\underset{n\to \mathrm{\infty }}{\lim }{(1+\frac{x}{n})}^{\frac{n}{x}}.$

Drugi pomysł okazuje się niewypałem między innymi dlatego, że musi korzystać z operacji potęgowania z wykładnikiem rzeczywistym (którą skądinąd najprościej implementować korzystając m.in. z funkcji $\exp$, jako $x^y=\exp (y\cdot \log (x))$).

Przyjrzyjmy się więc pierwszemu.

$e^x\approx \sum _{n=0}^N\frac{x^n}{n!},$

a dokładniej, ze wzoru na resztę w postaci Lagrange'a,

$|e^x-\sum _{n=0}^N\frac{x^n}{n!}|\le \frac{|\xi {|}^{N+1}}{(N+1)!}.$

Stąd dostajemy oszacowanie, oznaczając wartość szeregu obciętego na $N$-tym wyrazie,

$\exp (x,N)=\sum _{n=0}^N\frac{x^n}{n!},$ $|e^x-\exp (x,N)|\le \frac{|x{|}^{N+1}}{(N+1)!},$

tak więc $N$ musi być takie, że $\frac{|x{|}^{N+1}}{(N+1)!}\le ϵ$.

Poniżej mamy kod, który to realizuje w Octave (analogiczny program w C napisz samodzielnie):

function [y, N] = expa(x, epsilon)
y = skladnik = 1.0; N = 1;
while (abs(skladnik) > epsilon)
	skladnik *= (x/N);
	y += skladnik;
	N++;
end
end

Możesz go przetestować pod względem osiąganej dokładności i kosztu, porównując się z funkcją biblioteczną:

function [blad, relblad, koszt] = testexp(expX, x)
dokladnie = exp(x);
[wartosc, koszt] = feval(expX,x,1e-8);
blad = abs(wartosc-dokladnie);
relblad = blad/dokladnie;
end

i = 0; zakres = linspace(0,7,100); 
for x = zakres
	fprintf(stderr,'x = %e\n', x);
	i++;
	[blad, relblad, koszt] = testexp('expa',x); 
	koszta(i) = koszt; relblada(i) = relblad;
	blada(i) = blad;
end

xlabel('x');
ylabel('blad');
title(['Epsilon = ', num2str(epsilon)]);
plot(zakres, relblada, ';blad wzgledny;');
plot(zakres, blada, ';blad bezwzgledny;');

Zgodnie z oczekiwaniami, błąd jest poniżej zadanej tolerancji, tutaj $10^{-8}$.

Błąd względny aproksymacji wielomianem Taylora.

Koszt aproksymacji rośnie wraz z $x$:

Koszt aproksymacji wielomianem Taylora.

Niektóre wyniki mogą Cię jednak zaskoczyć:

• Błąd bezwzględny przekracza założoną wartość, gdy $ϵ=2.2\cdot 10^{-15}$, a błąd względny od pewnego momentu rośnie z $x$.

  

  Błąd względny aproksymacji wielomianem Taylora dla zadanej bardzo małej tolerancji błędu.

• Nie daje się tak policzyć $e^{2006}$.
• Wartość expa() może być ujemna dla $x<0$.

Wyjaśnienie tych szokujących faktów (które nie mają nic wspólnego z błędem w implementacji) musisz odłożyć do momentu, gdy bliżej przyjrzymy się temu, jak liczy komputer.

Ćwiczenie: Ciag dalszy

Spróbuj obniżyć koszt wyznaczania $\exp (x)$ dla dużych $x$.

Wskazówka

Rzecz w tym, że dla dużych $x$, trzeba wziąć bardzo dużo wyrazów w szeregu Taylora. Czy można tak wykombinować, by w rezultacie wziąć ich mniej?

Rozwiązanie

Bez zmieniejszenia ogólności, załóżmy, że $x\ge 0$.

Jeśli $x=k+t$, gdzie $t\in [0,1)$, a $k$ jest całkowite, to oczywiście

$e^x=e^{k+t}=e^k{e}^t=e^k\exp (t).$

Tak więc zadanie redukuje się do wyznaczenia $\exp (t)$ dla małego $t$ oraz do co najwyżej $k$ dodatkowych mnożeń potrzebnych do wyznaczenia całkowitej potęgi $e^k$ (ile mnożeń naprawdę wystarczy?). Pamiętaj, przyjęliśmy, że znamy reprezentację numeryczną liczby $e$.

# wersja B
function [y, N] = expb(x, epsilon)
k = floor(x); t = x - k; 
[y, N] = expa(t, epsilon); 
for i = 1:k
 	y *= e;
end
N += (k+2);
end

Wersja B jest istotnie tańsza.