MN01LAB

From Studia Informatyczne


Eksperymenty ze środowiskiem obliczeń numerycznych

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

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

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ć \displaystyle 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 \displaystyle \sin(\pi\cdot i\cdot 0.4), gdzie \displaystyle i=0,\ldots, 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: \displaystyle +,\, -, \, \times, \, \div, wartość funkcji exp(\displaystyle x) = \displaystyle e^x dla dowolnych \displaystyle x rzeczywistych. Naszym kryterium jest, by \displaystyle |e^x - \exp(x)| \leq \epsilon, czyli by błąd bezwzględny aproksymacji nie przekroczył zadanego \displaystyle \epsilon.

Wykonaj eksperymenty w C lub w Octave, pokazujące koszt metody w zależności od \displaystyle x oraz w zależności od \displaystyle \epsilon. 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 \displaystyle \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ą,

\displaystyle  e^x = \sum_{n=0}^\infty \frac{x^n}{n!}

lub o granicy

\displaystyle  e^x = \lim_{n\rightarrow \infty} \left( 1 + \frac{x}{n}\right)^{\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 \displaystyle \exp, jako \displaystyle x^y = \exp(y\cdot \log(x))).

Przyjrzyjmy się więc pierwszemu.

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

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

\displaystyle  |e^x - \sum_{n=0}^N \frac{x^n}{n!} | \leq \frac{|\xi|^{N+1}}{(N+1)!}.

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

\displaystyle  \exp(x,N) = \sum_{n=0}^N \frac{x^n}{n!},
\displaystyle  \left|e^x - \exp(x,N)  \right| \leq \frac{|x|^{N+1}}{(N+1)!},

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

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 \displaystyle 10^{-8}.

Błąd względny aproksymacji wielomianem Taylora.
Enlarge
Błąd względny aproksymacji wielomianem Taylora.

Koszt aproksymacji rośnie wraz z \displaystyle x:

Koszt aproksymacji wielomianem Taylora.
Enlarge
Koszt aproksymacji wielomianem Taylora.

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

  • Błąd bezwzględny przekracza założoną wartość, gdy \displaystyle \epsilon = 2.2\cdot 10^{-15}, a błąd względny od pewnego momentu rośnie z \displaystyle x.
    Błąd względny aproksymacji wielomianem Taylora dla zadanej bardzo małej tolerancji błędu.
    Enlarge
    Błąd względny aproksymacji wielomianem Taylora dla zadanej bardzo małej tolerancji błędu.
  • Nie daje się tak policzyć \displaystyle e^{2006}.
  • Wartość expa() może być ujemna dla \displaystyle 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 \displaystyle \exp(x) dla dużych \displaystyle x.

Wskazówka

Rzecz w tym, że dla dużych \displaystyle 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 \displaystyle x\geq 0.

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

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

Tak więc zadanie redukuje się do wyznaczenia \displaystyle \exp(t) dla małego \displaystyle t oraz do co najwyżej \displaystyle k dodatkowych mnożeń potrzebnych do wyznaczenia całkowitej potęgi \displaystyle e^k (ile mnożeń naprawdę wystarczy?). Pamiętaj, przyjęliśmy, że znamy reprezentację numeryczną liczby \displaystyle 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.
Enlarge
Wersja B jest istotnie tańsza.