|
|
| Linia 1: |
Linia 1: |
|
| |
|
|
| |
| ==Ćwiczenia. Własności arytmetyki zmiennopozycyjnej.==
| |
|
| |
| {{cwiczenie|||
| |
| Podaj przykłady ''zbieżnych'' szeregów postaci <math>\displaystyle \sum_{n=1}^{\infty} a_n</math>, którego <math>\displaystyle N</math>-te sumy częściowe
| |
| obliczone w arytmetyce pojedynczej precyzji algorytmem
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| suma <nowiki>=</nowiki> 0.0;
| |
| for n <nowiki>=</nowiki> 1..N
| |
| suma +<nowiki>=</nowiki> <math>\displaystyle a_n</math>;
| |
| </pre></div>
| |
|
| |
| będą
| |
| * ograniczone niezależnie od <math>\displaystyle N</math>, albo
| |
| * numerycznie rozbieżne, to znaczy takie, że dla bardzo dużych <math>\displaystyle N</math> zachodzi
| |
| <code>suma <nowiki>=</nowiki><nowiki>=</nowiki> Inf</code>.
| |
|
| |
| Wykonaj to samo zadanie, ale podając przykłady szeregów ''rozbieżnych'' (w
| |
| arytmetyce dokładnej).
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
|
| |
| Rozpatrz na przykład takie szeregi zbieżne:
| |
| * Szereg zerowy (numerycznie dostajesz oczywiście też 0)
| |
| * <math>\displaystyle 1 + 10^{2006} - 10^{2006} + 0 + \ldots = 1</math> oraz <math>\displaystyle 1 + 10^{2006} - 1 + 0 + \ldots = 10^{2006}</math>
| |
|
| |
| (numerycznie dostaniesz, odpowiednio, <code>NaN</code> (dlaczego?) i <code>Inf</code>.
| |
| * Szereg harmoniczny (numerycznie sumy częściowe w pewnym momencie przestaną
| |
| rosnąć).
| |
| * Szereg jedynkowy.
| |
|
| |
| </div></div>
| |
|
| |
| {{cwiczenie|||
| |
| Dla kolejnych <math>\displaystyle N</math>, wyznacz <math>\displaystyle N</math>-tą sumę częściową szeregu Taylora dla funkcji wykładniczej, gdy <math>\displaystyle x=-90</math>:
| |
| <center><math>\displaystyle
| |
| e^x \approx \sum_{n=0}^N \frac{x^n}{n!},
| |
| </math></center>
| |
|
| |
| korzystając z algorytmu podanego w poprzednim zadaniu. Porównaj błąd: względny
| |
| i bezwzględny w porównaniu
| |
| do wartości wyznaczonej z wykorzystaniem funkcji bibliotecznej <code>exp()</code>. Powtórz następnie dla <math>\displaystyle x=10</math>.
| |
|
| |
| Czy --- zgodnie z teorią matematyczną --- sumy te dążą do wartości
| |
| <math>\displaystyle e^{x}</math>. Objaśnij dokładnie, co się stało.
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
| Zastosowanie naszego algorytmu dla <math>\displaystyle x=-90</math> daje w wyniku (dla arytmetyki
| |
| podwójnej precyzji) sumę równą około <math>\displaystyle 10^{30}</math>,
| |
| gdy oczywiście naprawdę <math>\displaystyle \exp(-90)
| |
| \approx 0</math>). 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 <math>\displaystyle x=10</math> mamy, dla rosnącego <math>\displaystyle N</math>, całkiem niezły błąd względny ostatecznego
| |
| wyniku.
| |
|
| |
| </div></div>
| |
|
| |
| {{cwiczenie|||
| |
| Już wcześniej stwierdziliśmy, że wyznaczanie <math>\displaystyle e \approx (1 + 1/n)^n</math> dla dużego
| |
| <math>\displaystyle n</math> nie jest dobrym pomysłem. Przeprowadź eksperyment numeryczny potwierdzający
| |
| to stwierdzenie i objaśnij jego wyniki.
| |
| }}
| |
|
| |
| {{cwiczenie|||
| |
| Jak wiadomo, szereg harmoniczny <math>\displaystyle \sum_{n=1}^\infty 1/n</math> 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.
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| int dlicznik;
| |
| double dsuma, dstarasuma;
| |
| double dskladnik;
| |
|
| |
| dstarasuma <nowiki>=</nowiki> 0.0; dsuma <nowiki>=</nowiki> 1.0; dlicznik <nowiki>=</nowiki> 1;
| |
| while(dstarasuma !<nowiki>=</nowiki> dsuma)
| |
| {
| |
| dskladnik <nowiki>=</nowiki> 1.0/dlicznik;
| |
| dstarasuma <nowiki>=</nowiki> dsuma;
| |
| dsuma +<nowiki>=</nowiki> dskladnik;
| |
| dlicznik++;
| |
| }
| |
| printf("Suma <nowiki>=</nowiki> dsuma, dlicznik-1, dskladnik);
| |
|
| |
| return(0);
| |
| }
| |
|
| |
| </pre></div>
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskaz�wka </span><div class="mw-collapsible-content" style="display:none">
| |
| <div style="font-size:smaller; background-color:#efe"> Policz, ile mniej więcej obrotów pętli potrzeba, by suma przestała rosnąć. </div>
| |
| </div></div>
| |
|
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
| Nasza suma przekroczy wartość 20, więc aby <code>20 + dskladnik <nowiki>=</nowiki><nowiki>=</nowiki> 20</code>, musi
| |
| być dskładnik <math>\displaystyle \approx 10^{-15}</math> (lub więcej). Tymczasem zakres liczb typu integer wynosi
| |
| niewiele ponad <math>\displaystyle 10^9</math>, więc po jakimś czasie dlicznik przyjmie wartości ujemne.
| |
|
| |
| Aby doliczyć się "do końca", musielibyśmy licznik potraktować jako liczbę typu
| |
| <code>long long double</code>, którego zakres sięga <math>\displaystyle 10^{18}</math>. Ale i tak mielibyśmy
| |
| kłopot z doczekaniem do końca działalności programu.
| |
|
| |
| Zakładając optymistycznie, że na sekundę zliczamy <math>\displaystyle 10^7</math> składników,
| |
| potrzebowalibyśmy razem około <math>\displaystyle 10^7</math> sekund, czyli ponad trzy miesiące...
| |
|
| |
| </div></div>
| |
|
| |
| {{cwiczenie|||
| |
| Jak szybko i na jakiej wysokości musiałby lecieć samolot npla, aby
| |
| pocisk wystrzeliwany z działka z prędkością 7500 km/h nie trafił w cel, gdy
| |
| potrzebne pierwiastki liczone są wzorem szkolnym?
| |
| }}
| |
|
| |
| {{cwiczenie|Zadanie o wyznaczaniu pierwiastka kwadratowego metodą Newtona||
| |
|
| |
| Dla zadanej liczby <math>\displaystyle a>1</math>, należy wyznaczyć przybliżenie <math>\displaystyle \sqrt{a}</math> stosując metodę Herona. Zaproponuj
| |
| dobre przybliżenie początkowe <math>\displaystyle x_0</math> wiedząc, że <math>\displaystyle a</math> jest liczbą maszynową typu
| |
| <code>double</code>. Ile iteracji wystarczy, by osiągnąć <math>\displaystyle \epsilon</math>-zadowalające przybliżenie?
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
|
| |
| Jak pamiętamy, <math>\displaystyle a = (1+f)2^p = \widetilde{f} 2^{p+1}</math>,
| |
| gdzie <math>\displaystyle \frac{1}{2}\leq \widetilde{f} < 1</math> oraz <math>\displaystyle p\in N</math>. Stąd
| |
|
| |
| <center><math>\displaystyle
| |
| \sqrt{a} = \sqrt{\widetilde{f}}2^{(p+1)/2} = \sqrt{g} 2^k,
| |
| </math></center>
| |
|
| |
| dla pewnego (znanego) <math>\displaystyle k \in N</math>, przy czym <math>\displaystyle \frac{1}{2} \leq \sqrt{g} \leq 1</math>.
| |
| Wobec tego, obliczenie
| |
| <math>\displaystyle \sqrt{a}</math> wymaga po prostu obliczenia <math>\displaystyle \sqrt{g}</math>. Najprostsze przybliżenie <math>\displaystyle \sqrt{g}</math> to
| |
| <math>\displaystyle \sqrt{g}\approx 1</math>. Błąd takiego przybliżenia to
| |
| <math>\displaystyle |\sqrt{g} - 1| \leq 0.5</math>.
| |
|
| |
| Jak można jeszcze bardziej poprawić <math>\displaystyle x_0</math>?
| |
|
| |
| </div></div>
| |
|
| |
| {{cwiczenie|Zadanie o wyznaczaniu odwrotności bez dzielenia metodą Newtona||
| |
|
| |
| Należy wyznaczyć przybliżenie <math>\displaystyle \frac{1}{a}</math> stosując metodę Newtona do równania
| |
| <math>\displaystyle \frac{1}{x} - a = 0</math>. Zaproponuj
| |
| dobre przybliżenie początkowe <math>\displaystyle x_0</math> wiedząc, że <math>\displaystyle a</math> jest liczbą maszynową typu
| |
| <code>double</code>. Ile iteracji wystarczy, by osiągnąć <math>\displaystyle \epsilon</math>-zadowalające przybliżenie?
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
| Jak pamiętamy, <math>\displaystyle a = (1+f)2^p</math>, skąd <math>\displaystyle \frac{1}{a} = \frac{1}{1+f}2^{-p}</math>.
| |
| Wystarczy więc przybliżyć <math>\displaystyle \frac{1}{1+f}</math>. Ponieważ dla <math>\displaystyle 0 \leq f < 1</math>,
| |
|
| |
| <center><math>\displaystyle \frac{1}{2} \leq \frac{1}{1+f} \leq 1,
| |
| </math></center>
| |
|
| |
| to możemy jako pierwsze przybliżenie wybrać <math>\displaystyle x_0 = \frac{3}{4}2^{-p}</math>. (Jak
| |
| wybrać lepsze <math>\displaystyle x_0</math>?)
| |
|
| |
| Wtedy mamy następujące oszacowania:
| |
|
| |
| <center><math>\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}.
| |
| </math></center>
| |
|
| |
| Ponieważ iteracja Newtona spełnia
| |
|
| |
| <center><math>\displaystyle |x_{k+1}a - 1| = |x_ka-1|^2 = |x_0a-1|^{2^k},
| |
| </math></center>
| |
|
| |
| to, jeśli chcemy, by błąd względny przybliżenia <math>\displaystyle x_{k+1}</math> spełniał <math>\displaystyle \frac{|x_{k+1} -
| |
| \frac{1}{a}|}{\frac{1}{a}} = |x_{k+1}a - 1| \leq \epsilon = 2^{-53}</math>, to
| |
| musimy wykonać co najwyżej <math>\displaystyle 6</math> iteracji. Ponieważ koszt jednej iteracji
| |
| to 3 flopy, zatem wyznaczenie odwrotności <math>\displaystyle a</math> 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ć,
| |
| [http://www.cs.utexas.edu/users/moore/publications/divide_paper.pdf 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
| |
| [http://www.informatik.uni-trier.de/Reports/TR-08-2004/rnc6_12_markstein.pdf Software
| |
| Division and Square Root Using Goldschmidt s Algorithms] oraz
| |
| [http://arith.stanford.edu/techrep.html raportów technicznych Stanford Computer Architecture and Arithmetic Group
| |
| ]
| |
| </div></div>
| |
|
| |
| {{cwiczenie|||
| |
| Niech <math>\displaystyle 0<a_1<a_2<\cdots < a_n</math>. Czy z punktu widzenia
| |
| błędów w <math>\displaystyle fl_\nu</math> lepiej jest policzyć sumę tych liczb w kolejności
| |
| od najmniejszej do największej czy odwrotnie?
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
| Od najmniejszej do największej.
| |
| </div></div>
| |
|
| |
| {{cwiczenie|||
| |
| Aby obliczyć <math>\displaystyle S(a,b)=a^2-b^2</math> można zastosować
| |
| dwa algorytmy: <math>\displaystyle {\bf ALG}_1(a,b)=a*a-b*b</math> oraz <math>\displaystyle {\bf ALG}_2(a,b)=(a+b)*(a-b)</math>.
| |
| Pokazać, że oba algorytmy są numerycznie poprawne, ale drugi
| |
| z nich wywołuje mniejszy błąd względny wyniku w przypadku, gdy
| |
| <math>\displaystyle rd_\nu(a)=a</math> i <math>\displaystyle rd_\nu(b)=b</math>.
| |
| }}
| |
|
| |
| {{cwiczenie|||
| |
| Pokazać, że naturalny algorytm obliczania cosinusa
| |
| kąta między dwoma wektorami <math>\displaystyle a, b\inR^n</math>,
| |
|
| |
| <center><math>\displaystyle \cos(a,b)\,=\,\frac{\sum_{j=1}^n a_jb_j}
| |
| {\sqrt{\Big(\sum_{j=1}^n a_j^2\Big)
| |
| \Big(\sum_{j=1}^n b_j^2\Big)}},
| |
| </math></center>
| |
|
| |
| jest numerycznie poprawny. Oszacować błąd względny wyniku
| |
| w <math>\displaystyle fl_\nu</math>.
| |
| }}
| |
|
| |
| {{cwiczenie|||
| |
| Pokazać, że naturalny algorytm obliczania
| |
| <math>\displaystyle \|A x\|_2</math> dla danej macierzy <math>\displaystyle A\inR^{n\times n}</math> i wektora
| |
| <math>\displaystyle x\inR^n</math> jest numerycznie poprawny. Dokładniej,
| |
|
| |
| <center><math>\displaystyle fl_\nu(\|A x\|_2)\,=\,(A+E) x,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \|E\|_2\leq 2(n+2)\sqrt n\nu\|A\|_2</math>. Ponadto, jeśli
| |
| <math>\displaystyle A</math> jest nieosobliwa to
| |
|
| |
| <center><math>\displaystyle |fl_\nu(\|A x\|_2)-\|A x\|_2|\,\leq\,2(n+2)\sqrt{n}\,\nu\,
| |
| \left(\|A\|_2\|A^{-1}\|_2\right)\,\|A x\|_2.
| |
| </math></center>
| |
|
| |
| }}
| |
|
| |
| {{cwiczenie|||
| |
| Niech <math>\displaystyle {\bf ALG}</math> będzie algorytmem numerycznie
| |
| poprawnym w zbiorze danych <math>\displaystyle f\in F_0</math>, przy czym dla małych <math>\displaystyle \nu</math>,
| |
| <math>\displaystyle fl_\nu({\bf ALG}(f))=\varphi(y_\nu)</math>, gdzie <math>\displaystyle \|y_\nu-y\|\le K\nu\|y\|</math>
| |
| i <math>\displaystyle K</math> nie zależy od <math>\displaystyle \nu</math> i <math>\displaystyle f</math> (<math>\displaystyle y=N(f)</math>). Pokazać, że
| |
| w ogólności <math>\displaystyle {\bf ALG}</math> nie musi być "numerycznie poprawny po
| |
| współrzędnych", tzn. w ogólności nie istnieje bezwzględna
| |
| stała <math>\displaystyle K_1</math> taka, że dla małych <math>\displaystyle \nu</math> i dla dowolnej
| |
| <math>\displaystyle f\in F_0</math>
| |
|
| |
| <center><math>\displaystyle |y_{\nu,j}-y_j|\,\le\,K_1\,\nu\,|y_j|, \qquad 1\le j\le n,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle y=(y_1,\ldots,y_n)</math>.
| |
| }}
| |
|
| |
| {{cwiczenie|||
| |
| Podaj przykład funkcji <math>\displaystyle f</math>, której miejsce zerowe <math>\displaystyle x^*</math> ma wspólczynnik
| |
| uwarunkowania
| |
| * mały
| |
| * duży
| |
|
| |
| }}
| |
|
| |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwi�zanie </span><div class="mw-collapsible-content" style="display:none">
| |
| Ponieważ nasze zadanie to wyznaczenie <math>\displaystyle x^* = f^{-1}(0)</math>, to
| |
| <center><math>\displaystyle
| |
| \mbox{cond} _{abs} (f^{-1},0) = \frac{1}{f'(x^*)}.
| |
| </math></center>
| |
|
| |
| Znaczy to, że im bardziej płaska jest <math>\displaystyle f</math> w otoczeniu pierwiastka <math>\displaystyle x^*</math>, tym
| |
| bardziej nawet małe zaburzenia <math>\displaystyle f</math> mogą spowodować duże przemieszczenie jej
| |
| miejsca zerowego.
| |
|
| |
| Zauważ, że dla wielokrotnych miejsc zerowych, <math>\displaystyle \mbox{cond} _{abs} (f^{-1},0) = \infty</math>,
| |
| co zgadza się z intuicją, bo może być, że nawet minimalne zaburzenie <math>\displaystyle f</math>
| |
| spowoduje, iż miejsc zerowych po prostu nie będzie...
| |
| </div></div>
| |
