|
|
| Linia 1: |
Linia 1: |
|
| |
|
| ==Arytmetyka zmiennopozycyjna==
| |
|
| |
| Metody iteracyjne mają czasem kłopoty, które nie są związane z samą naturą
| |
| problemu matematycznego. Przyrzyjmy się bowiem, jak w dużym zbliżeniu wygląda
| |
| wykres funkcji <math>\displaystyle w(x) = x^4-4x^3+6x^2-4x+1</math>, której wartości zostały obliczone na
| |
| komputerze PC. Nietrudno sprawdzić, że <math>\displaystyle w</math> ma dokładnie jedno miejsce zerowe,
| |
| gdyż <math>\displaystyle w(x)=(x-1)^4</math>. Tymczasem, wykres <math>\displaystyle w</math> (wyznaczony oryginalnym wzorem) zdaje
| |
| się mieć ''mnóstwo'' różnych miejsc zerowych w okolicy <math>\displaystyle x=1</math>. Co gorsza,
| |
| wygląda na to, że <math>\displaystyle w</math> wcale nie jest gładka!
| |
|
| |
| [[Image:MNwielomian4.png|thumb|300px|Wartości funkcji <math>\displaystyle w(x) = x^4-4x^3+6x^2-4x+1</math> obliczone według wzoru. Na marginesie: <math>\displaystyle w(x) =
| |
| (x-1)^4</math>. Prawdziwe wartości zaznaczone na czerwono.]]
| |
|
| |
| Wykonywanie realnych obliczeń na liczbach rzeczywistych w komputerze może być
| |
| źródłem wielu innych zaskoczeń. Na przykład, w komputerze,
| |
|
| |
| <center><math>\displaystyle
| |
| 10\cdot (1.1 - 1) \neq 1
| |
| </math></center>
| |
|
| |
| co możesz łatwo sprawdzić:
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| octave:7> 10 * (1.1 -1) - 1
| |
| ans <nowiki>=</nowiki> 8.8818e-16
| |
| </pre></div>
| |
|
| |
| Przedstawiony wcześniej model obliczeniowy jest modelem idealistycznym, tzn.
| |
| zakłada on, że wszystkie operacje są wykonywane bezbłędnie.
| |
| Dlatego w tym przypadku będziemy mówić o ''arytmetyce idealnej''.
| |
| W praktyce jednak, np. wykonując obliczenia na maszynie cyfrowej,
| |
| operacje arytmetyczne na liczbach rzeczywistych wykonywane są
| |
| z pewnym błędem. Matematycznym modelem arytmetyki maszyny cyfrowej
| |
| jest ''arytmetyka <math>\displaystyle fl_\nu</math>'' (albo arytmetyka
| |
| ''zmiennoprzecinkowa''), którą teraz przedstawimy.
| |
|
| |
| Niech będzie zadana liczba naturalna <math>\displaystyle b</math> (jej znaczenie wyjaśni się w następnym
| |
| rozdziale). Dowolną liczbę rzeczywistą <math>\displaystyle x\ne 0</math> można jednoznacznie przedstawić w postaci
| |
|
| |
| <center><math>\displaystyle x\,=\,s\cdot 2^{c-b}\cdot m, </math></center>
| |
|
| |
| gdzie <math>\displaystyle s\in\{-1,1\}</math> jest znakiem, liczba całkowita
| |
| <math>\displaystyle (c-b)</math> ''cechą'', a liczba rzeczywista <math>\displaystyle m\in [1,2)</math> ''mantysą'' liczby <math>\displaystyle x</math>.
| |
| Zauważmy, że taki
| |
| rozkład jest jednoznaczny i odpowiada przesuwaniu przecinka w rozwinięciu
| |
| binarnym liczby do pierwszej cyfry znaczącej, tj. różnej od zera (stąd nazwa:
| |
| reprezentacja zmiennoprzecinkowa, ang. floating point). Mantysa ma w ogólności
| |
| nieskończenie wiele cyfr binarnych <math>\displaystyle f_j</math> w swoim rozwinięciu dwójkowym,
| |
|
| |
| <center><math>\displaystyle m =
| |
| 1 + f \equiv 1 + \sum_{j=1}^\infty f_j 2^{-j} = (1.f_1f_2f_3\ldots)_2, </math></center>
| |
|
| |
| gdzie
| |
| <math>\displaystyle f_j\in\{0,1\}</math>. Wobec tego najczęściej nie będzie mogła być zapamiętana
| |
| dokładnie w pamięci komputera, gdyż możemy przechować jedynie ''ograniczoną''
| |
| liczbę cyfr cechy i mantysy.
| |
|
| |
| ===Reprezentacja zmiennoprzecinkowa===
| |
|
| |
| W komputerach osobistych mamy do czynienia z reprezentacją liczb rzeczywistych,
| |
| w której do zapisania liczby używa się ściśle określonej liczby bitów <math>\displaystyle t</math> do
| |
| zapisania mantysy i także określonej liczby bitów <math>\displaystyle p</math> do zapisania cechy danej
| |
| liczby niezerowej <math>\displaystyle x</math>:
| |
|
| |
| <center><math>\displaystyle
| |
| s\,c_1c_2\ldots c_p\,f_1f_2\ldots f_t
| |
| </math></center>
| |
|
| |
| (łącznie <math>\displaystyle 1+p+t</math> bitów). Liczby zapisane przy użyciu powyższej sekwencji bitów
| |
| nazywa się ''liczbami maszynowymi''. Są to jedyne dokładnie zapisywalne w
| |
| komputerze liczby rzeczywiste, pozostałe będą musiały zostać wyrażone z
| |
| wykorzystaniem liczb maszynowych.
| |
|
| |
| ''Reprezentacją zmiennoprzecinkową'' niezerowej
| |
| liczby <math>\displaystyle x</math> będziemy nazywać liczbę <math>\displaystyle rd_\nu(x)</math> taką, że
| |
|
| |
| <center><math>\displaystyle
| |
| rd_\nu(x) = (-1)^s \cdot (1+f)\cdot 2 ^{c-b},
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle f</math> jest liczbą dwójkową postaci <math>\displaystyle (0.f_1\ldots f_{t})_2</math>, natomiast <math>\displaystyle c</math>
| |
| jest liczbą naturalną postaci <math>\displaystyle (c_1\ldots c_p)_2</math>. Na znak liczby, <math>\displaystyle s</math>,
| |
| przeznaczony jest jeden bit. Wartości <math>\displaystyle c</math> i <math>\displaystyle f</math> dobiera się tak, żeby <math>\displaystyle rd_\nu(x)</math>
| |
| była tak bliska <math>\displaystyle x</math> jak to możliwe. Stałą całkowitą <math>\displaystyle b</math> dobiera się tak, by
| |
| uzyskać zbalansowany zakres cechy <math>\displaystyle c-b</math> (mniej więcej tyle samo wartości
| |
| ujemnych i dodatnich), a zysk z korzystania z niej jest taki, że nie marnujemy
| |
| dodatkowego bitu na przechowywanie znaku wykładnika potęgi dwójki <math>\displaystyle c-b</math>.
| |
|
| |
| {{przyklad|||
| |
| Rozważmy bardzo prościutki system, w którym zarówno na cechę, jak i mantysę,
| |
| przeznaczone są jedynie po dwa bity, zatem jedna liczba maszynowa zajmuje 5
| |
| bitów. Ponieważ w konsekwencji możliwy zakres wartości <math>\displaystyle c</math>
| |
| to <math>\displaystyle 0,\ldots, 3</math>, rozsądne jest więc przyjęcie korekty <math>\displaystyle b = 1</math>, dzięki czemu
| |
| <math>\displaystyle -1 \leq c-b \leq 2</math>. Z kolei możliwe wartości mantysy to
| |
| <center><math>\displaystyle
| |
| (1.00)_2 = 1,\qquad (1.01)_2 = 1.25,\qquad (1.10)_2 = 1.5,\qquad (1.11)_2 =
| |
| 1.75.
| |
| </math></center>
| |
|
| |
| Wobec tego, jedyne (dodatnie) liczby maszynowe naszej pięciobitowej arytmetyki
| |
| zmiennopozycyjnej to
| |
|
| |
| {| border=1
| |
| |+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
| |
| |-
| |
| | 0.500 || 0.625 || 0.750 || 0.875
| |
| |-
| |
| | 1.000 || 1.250 || 1.500 || 1.750
| |
| |-
| |
| | 2.000 || 2.500 || 3.000 || 3.500
| |
| |-
| |
| | 4.000 || 5.000 || 6.000 || 7.000
| |
|
| |
| |}
| |
|
| |
| [[Image:MNbinarysystem.png|thumb|300px|Liczby maszynowe: reprezentowane dokładnie w
| |
| pięciobitowej arytmetyce o
| |
| precyzji <math>\displaystyle 2^{-2}</math>.]]
| |
|
| |
| }}
| |
|
| |
| Przy takim sposobie reprezentacji, jej błąd względny szacuje się przez
| |
|
| |
| <center><math>\displaystyle
| |
| \left|\frac{rd_\nu(x) - x }{x}\right| \leq \frac{1}{2^{t+1}}.
| |
| </math></center>
| |
|
| |
| Liczbę <math>\displaystyle \nu = \frac{1}{2^{t+1}}</math> nazywa się precyzją arytmetyki. Jak widać, ma
| |
| ma ona wpływ na to, jak wiele cyfr znaczących liczby jest reprezentowanych
| |
| dokładnie. Precyzja arytmetyki zależy wyłącznie od liczby bitów przeznaczonych
| |
| na reprezentację mantysy.
| |
|
| |
| Ostatnią nierówność wygodnie jest zapisać w równoważny
| |
| sposób jako
| |
|
| |
| <center><math>\displaystyle rd_\nu(x)\,=\,x(1+\epsilon), \qquad \mbox{gdzie} \quad |\epsilon|\le\nu.
| |
| </math></center>
| |
|
| |
| ====Standard IEEE 754====
| |
|
| |
| Z nielicznymi egzotycznymi wyjątkami (np. Cray C90), współczesne procesory
| |
| używane w komputerach osobistych lub stacjach roboczych implementują
| |
| [??? IEEE 754 Floating Point Standard], który definiuje dwa zasadnicze
| |
| formaty reprezentacji zmiennoprzecinkowej liczb rzeczywistych:
| |
|
| |
| {| border=1
| |
| |+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
| |
| |-
| |
| | Typ IEEE 754 || Pojedycznej precyzji || Podwójnej precyzji
| |
| |-
| |
| | Nazwa typu w C || float || double
| |
| |-
| |
| | Liczba bitów cechy || 8 || 11
| |
| |-
| |
| | Liczba bitów mantysy || 23 || 52
| |
| |-
| |
| | Liczba bajtów dla typu w C || 4 || 8
| |
| |-
| |
| | Bias (liczba <math>\displaystyle b</math> powyżej) || 127 || 1023
| |
| |-
| |
| | Orientacyjny zakres || <math>\displaystyle 10^{-38}\ldots 10^{+38}</math> || <math>\displaystyle 10^{-308}\ldots 10^{+308}</math>
| |
| |-
| |
| | Orientacyjna precyzja || <math>\displaystyle 6\cdot 10^{-8}</math> || <math>\displaystyle 10^{-16}</math>
| |
| |-
| |
| |
| |
|
| |
| |}
| |
|
| |
| (maksymalna i minimalna wartość cechy <math>\displaystyle c</math> ma specjalne znaczenie). Dodatkowo, w
| |
| procesorach x86 mamy typ podwójnej rozszerzonej precyzji (także
| |
| zdefiniowany w IEEE 754 i odpowiadający dokładnie ówczesnym możliwościom
| |
| procesora Intel 8087). Wszystkie operacje arytmetyczne na procesorach x86
| |
| są faktycznie wykonywane w takiej precyzji (korzystając z 64 bitów dla
| |
| reprezentacji mantysy i 15 bitów dla cechy). Należy pamiętać, że odpowiadający
| |
| mu typ w C <code>long double</code> zajmuje w pamięci 12 bajtów (a nie 80 bitów).
| |
|
| |
| {{uwaga|||
| |
| Producenci niektórych procesorów świadomie rezygnują z implementacji IEEE 754
| |
| dla zwiększenia szybkości działania kosztem niestety dokładności wyniku. Tak
| |
| dawno temu było w procesorach Cray;
| |
| tak obecnie działa np. procesor IBM Cell (stosowany w Playstation); tak też działają niektóre instrukcje wektorowe (z tzw.
| |
| zestawu 3DNow!) w
| |
| procesorach AMD, które np. wynik dzielenia wektorowego zwracają z precyzją tylko
| |
| 14 bitów mantysy.
| |
| }}
| |
|
| |
| W Octave można łatwo podejrzeć reprezentację binarną liczby zmiennoprzecinkowej
| |
| podwójnej precyzji (jest to domyślny typ numeryczny stosowany w MATLABie i
| |
| Octave),
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| octave:9> format bit
| |
| octave:10> x <nowiki>=</nowiki> -2
| |
| x <nowiki>=</nowiki> 1100000000000000000000000000000000000000000000000000000000000000
| |
| octave:11> x <nowiki>=</nowiki> 1/4
| |
| x <nowiki>=</nowiki> 0011111111010000000000000000000000000000000000000000000000000000
| |
| octave:12> x <nowiki>=</nowiki> NaN
| |
| x <nowiki>=</nowiki> 1111111111111000000000000000000000000000000000000000000000000000
| |
| octave:13> x <nowiki>=</nowiki> 0
| |
| x <nowiki>=</nowiki> 0000000000000000000000000000000000000000000000000000000000000000
| |
| octave:14> x <nowiki>=</nowiki> Inf
| |
| x <nowiki>=</nowiki> 0111111111110000000000000000000000000000000000000000000000000000
| |
| octave:15> x <nowiki>=</nowiki> 0.1
| |
| x <nowiki>=</nowiki> 0011111110111001100110011001100110011001100110011001100110011010
| |
| </pre></div>
| |
|
| |
| (w MATLABie możemy zobaczyć tę samą liczbę w zapisie szestnastkowym).
| |
|
| |
| {{przyklad|Nawet liczba 0.1 nie jest reprezentowana dokładnie!||
| |
|
| |
| Rozwinięcie dwójkowe liczby 0.1 jest nieskończone:
| |
|
| |
| <center><math>\displaystyle
| |
| 0.1 = (0.0001 1001 1001 1001 \ldots)_2.
| |
| </math></center>
| |
|
| |
| Ten banalny fakt jest bardzo często przeoczany przez programistów, a w 199x roku
| |
| doprowadził do spektakularnej awarii systemu antyrakietowego Patriot. Okazało
| |
| się, że --- w tajemniczy sposób --- zazwyczaj bezbłędnie trafiające w cel
| |
| rakiety Patriot traciły skuteczność, gdy przez wiele godzin pozostawały w stanie
| |
| gotowości.
| |
|
| |
| Wyjaśnienie zagadki leżało na styku pomiędzy hardware a software rakiety. Jak
| |
| zbadano, w celu pomiaru czasu, zliczano kolejne tyknięcia zegara rakiety, które
| |
| następowały dokładnie co 0.1 sekundy. Następnie, w celu wyznaczenia prawdziwego
| |
| czasu, mnożono liczbę tyknięć zegara przez 0.1 (które właśnie było niedokładnie
| |
| reprezentowane). Gdy cykli zegara było bardzo dużo, błąd bezwzględny wyznaczenia
| |
| czasu stawał się również na tyle poważny, że uniemożliwiał precyzyjne
| |
| wyznaczenie parametrów toru lotu nieprzyjacielskiego obiektu!
| |
|
| |
| Na marginesie zauważmy, że np. liczba <math>\displaystyle 0.125</math> ''jest reprezentowana
| |
| dokładnie'' w arytmetyce zmiennoprzecinkowej (dlaczego?) i nie powodowałaby już
| |
| tego problemu.
| |
|
| |
| }}
| |
|
| |
| {{uwaga|||
| |
| Nie wszystkie maszyny liczące wykorzystują reprezentację dwójkową. Kiedyś
| |
| zdarzały się komputery reprezentujące liczby w postaci
| |
|
| |
| <center><math>\displaystyle x\,=\,s\cdot \beta^c\cdot m,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \beta = 8</math> lub 16, a nawet 3 (sic!). Do dzisiaj, w podręcznych
| |
| kalkulatorach najczęściej spotykaną podstawą reprezentacji liczb jest <math>\displaystyle \beta =
| |
| 10</math>.
| |
|
| |
| Są także takie realizacje arytmetyki zmiennoprzecinkowej, które nie realizują w
| |
| pełni standardu IEEE (np. stare komputery Cray) i np. zamiast zaokrąglenia,
| |
| stosują obcięcie wyniku.
| |
|
| |
| }}
| |
|
| |
| ====Nadmiar i niedomiar====
| |
|
| |
| W maszynie cyfrowej cecha <math>\displaystyle c</math> liczby rzeczywistej
| |
| nie może oczywiście mieć dowolnie dużej wartości bezwzględnej,
| |
| <math>\displaystyle |c|\le c_{\max}</math>, dlatego nie wszystkie liczby rzeczywiste są w ogóle ''reprezentowalne''. Powoduje to powstanie zjawiska ''nadmiaru'' gdy dla liczby
| |
| <math>\displaystyle x\displaystyle c>c_{\max}</math>, oraz zjawiska ''niedomiaru'' gdy <math>\displaystyle c<-c_{\min}</math>. W
| |
| pierwszym przypadku liczba jest tak duża (co do modułu), że
| |
| nie zawiera się w przedziale liczb reprezentowalnych, a w drugim
| |
| jest tak mała, że musi być reprentowana przez zero, przy czym
| |
| błąd względny reprezentacji wynosi wtedy <math>\displaystyle 1</math> a nie <math>\displaystyle \nu</math>.
| |
|
| |
| [[Image:MNbinarysystem1emptyspace.png|thumb|300px|Próżnia wokół zera]]
| |
|
| |
| Arytmetyka IEEE 754 przyjmuje, że liczby dla których następuje overflow są
| |
| reprzentowane przez specjalną wartość <code>Inf</code> (nieskończoność, ze znakiem), która propaguje się w
| |
| obliczeniach.
| |
|
| |
| [[Image:MNbinarysystem2infinity.png|thumb|300px|Wszystkie liczby większe od największej
| |
| zapisywalnej liczby są reprezentowane przez <code>Inf</code>]]
| |
|
| |
| W dalszych rozważaniach zjawiska nadmiaru i niedomiaru będziemy
| |
| dla uproszczenia zaniedbywać, jednak nie zawsze jest to uzasadnione, o czym
| |
| niech świadczy poniższy przykład.
| |
|
| |
| {{przyklad|Wyznaczanie normy euklidesowej wektora||
| |
|
| |
| Jedną z najczęściej wykonywanych operacji na wektorze <math>\displaystyle x = (x_1,\ldots,x_n)^T \in R^n</math> jest obliczenie jego
| |
| normy euklidesowej,
| |
| <center><math>\displaystyle
| |
| ||x||_2 = \sqrt{x_1^2 + \ldots x_n^2}.
| |
| </math></center>
| |
|
| |
| Jak widać, możemy tu łatwo zetknąć się ze zjawiskiem zarówno niedomiaru, jak i
| |
| nadmiaru, gdyż może się na przykład tak złożyć, że mimo iż <math>\displaystyle ||x||_2</math> jest
| |
| reprezentowana, to <math>\displaystyle x_1^2</math> już nie (np. w arytmetyce podwójnej precyzji <math>\displaystyle x_1 =
| |
| 10^{200}</math> i <math>\displaystyle x_2 = 1</math>).
| |
|
| |
| Łatwym wyjściem z tej sytuacji jest wstępna ''normalizacja danych'' tak, by
| |
| wszystkie nie były większe od 1: niech <math>\displaystyle M = \max\{|x_i|: i = 1,\ldots,n\}</math> i
| |
| wtedy
| |
| <center><math>\displaystyle
| |
| ||x||_2 = \sqrt{x_1^2 + \ldots x_n^2} = M\cdot\sqrt{\left(\frac{x_1}{M}\right)^2
| |
| + \ldots + \left(\frac{x_1}{M}\right)^2}.
| |
| </math></center>
| |
|
| |
| i teraz suma pod pierwiastkiem jest zawsze pomiędzy 1 a <math>\displaystyle N</math>. Wadą omówionego
| |
| rozwiązania jest to, że wymaga ono dwukrotnego przejrzenia całego wektora (raz,
| |
| by znaleźć <math>\displaystyle M</math>, drugi raz --- by policzyć sumę. Na szczęście można go
| |
| zmodyfikować tak, by działał w jednym przebiegu. Zupełnie inny algorytm podał
| |
| Moler .
| |
| }}
| |
|
| |
| ====Liczby denormalizowane====
| |
|
| |
| Wymaganie, że mantysa jest postaci <math>\displaystyle 1+f</math>, <math>\displaystyle f\geq 0</math>, powoduje, że wokół zera
| |
| pojawia się coś w rodzaju próżni. Formalnie, liczby mniejsze niż <math>\displaystyle 2^{1-1023}</math>
| |
| powinny być reprezentowane przez 0, lecz zazwyczaj zamiast tego
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| octave:16> format bit
| |
| octave:17> x <nowiki>=</nowiki> 2^(-1022)
| |
| x <nowiki>=</nowiki> 0000000000010000000000000000000000000000000000000000000000000000
| |
| octave:18> x <nowiki>=</nowiki> 2^(-1023)
| |
| x <nowiki>=</nowiki> 0000000000001000000000000000000000000000000000000000000000000000
| |
| octave:19> x <nowiki>=</nowiki> 2^(-1028)
| |
| x <nowiki>=</nowiki> 0000000000000000010000000000000000000000000000000000000000000000
| |
| </pre></div>
| |
|
| |
| W ten sposób można (w podwójnej precyzji) zbliżyć się do zera na odległość około
| |
| <math>\displaystyle 10^{-323}</math>.
| |
|
| |
| [[Image:MNbinarysystem3denormals.png|thumb|300px|Liczby denormalizowane trochę wypełniają
| |
| próżnię wokół zera]]
| |
|
| |
| ====Działania arytmetyczne w arytmetyce <math>\displaystyle fl_\nu</math>====
| |
|
| |
| W arytmetyce <math>\displaystyle fl_\nu</math> implementującej standard IEEE 754, działania arytmetyczne
| |
| na liczbach rzeczywistych (a raczej na ich reprezentacjach) są
| |
| wykonywane dokładnie i tylko wynik jest zaokrąglany. Mamy więc
| |
|
| |
| <center><math>\displaystyle fl_\nu(x\,\Box\,y)\,=\,rd_\nu\left(\,rd_\nu(x)\,\Box\,rd_\nu(y)\,\right),
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \Box\in\{+,-,\times,\div\}</math>, Ogólniej, jeśli <math>\displaystyle {\cal W}_1</math> i
| |
| <math>\displaystyle {\cal W}_2</math> są wyrażeniami o wartościach rzeczywistych, to
| |
| dla dowolnych wartości zmiennych
| |
|
| |
| <center><math>\displaystyle fl_\nu({\cal W}_1\,\Box\,{\cal W}_2)\,=\,
| |
| rd_\nu\left(\,fl_\nu({\cal W}_1)\,\Box\,fl_\nu({\cal W}_2)\right).
| |
| </math></center>
| |
|
| |
| Zwykle dla prostoty będziemy również zakładać podobną
| |
| zależność dla niektórych funkcji standardowych, o ile należą
| |
| one do zbioru operacji elementarnych (chociaż w rzeczywistości
| |
| są one obliczane przez procedury używające czterech podstawowych
| |
| operacji arytmetycznych). I tak będziemy mieć np.
| |
|
| |
| <center><math>\displaystyle \aligned fl_\nu\Big(\sqrt{{\cal W}}\Big) &= \left(\sqrt{fl_\nu({\cal W})}\right)
| |
| (1+\beta_1),\\
| |
| fl_\nu(\cos({\cal W})) &= \left(\cos(fl_\nu({\cal W}))\right)(1+\beta_2),
| |
| \endaligned</math></center>
| |
|
| |
| gdzie <math>\displaystyle |\epsilon_j|\le\nu</math>, oraz <math>\displaystyle \beta_j\le K_j\nu</math> i <math>\displaystyle K_j</math> są
| |
| "niewielkimi" stąłymi.
| |
|
| |
| Przypuśćmy, że w naszym prościutkim pięciobitowym systemie spełniającym wymogi
| |
| standardu IEEE 754 zechcemy wykonać mnożenie
| |
| <center><math>\displaystyle
| |
| 1.3 \cdot 2.4
| |
| </math></center>
| |
|
| |
| Poniżej możemy przekonać się, jak będzie ono przebiegać.
| |
|
| |
| <div class="thumb tright"><div><flash>file=binarysystem4.swf</flash><div.thumbcaption>Mnożenie dwóch liczb rzeczywistych</div></div></div>
| |
|
| |
| \beginnowiki
| |
| [[Image:MNbinarysystem41.png|thumb|300px|Liczba 1.3 nie jest dokładnie reprezentowalna w
| |
| naszym systemie]]
| |
| [[Image:MNbinarysystem42.png|thumb|300px|Jej reprezentacja to najbliższa jej liczba
| |
| maszynowa --- 1.25]]
| |
| [[Image:MNbinarysystem43.png|thumb|300px|Również drugi czynnik, 2.4, nie jest liczbą
| |
| maszynową]]
| |
| [[Image:MNbinarysystem44.png|thumb|300px|A więc jego reprezentacją będzie znów najbliższa mu
| |
| liczba maszynowa.]]
| |
| [[Image:MNbinarysystem45.png|thumb|300px|Mnożenie odbywa się już na reprezentacjach obu
| |
| czynników]]
| |
| [[Image:MNbinarysystem46.png|thumb|300px|Wynik dokładnego mnożenia tych liczb maszynowych to
| |
| 3.125 --- znowu musi być zaokrąglony... ]]
| |
| [[Image:MNbinarysystem47.png|thumb|300px|...do najbliższej liczby maszynowej. Ostatecznie,
| |
| błąd względny wyniku wynosi około <math>\displaystyle 10^{-3}</math> i jest znacznie mniejszy niż
| |
| pesymistyczne oszacowanie (czasem osiągalne, lecz nie tym razem!) <math>\displaystyle 2^{-3}
| |
| \approx 10^{-1}</math>]]
| |
| \endnowiki
| |
|
| |
| Podobnie, jeśli <math>\displaystyle \triangle</math> jest operatorem porównania,
| |
| <math>\displaystyle \triangle\in\{<,\le,=,\ne\}</math>, to wartością wyrażenia
| |
| logicznego <math>\displaystyle {\cal W}_1\triangle {\cal W}_2</math> w <math>\displaystyle fl_\nu</math> jest
| |
| dokładna wartość wyrażenia
| |
| <math>\displaystyle fl_\nu({\cal W}_1)\trianglefl_\nu({\cal W}_2)</math>.
| |
|
| |
| Dosyć dziwnie w porównaniach zachowuje się specjalna liczba <code>NaN</code>
| |
| (ang. not-a-number), dla której zawsze zachodzi, że <code>NaN</code><math>\displaystyle \neq</math><code>NaN</code>.
| |
| Liczba <code>NaN</code> pojawia się jako wynik zabronionych operacji matematycznych,
| |
| np. <math>\displaystyle 0/0, \sqrt{-2},</math> <code>Inf</code> - <code>Inf</code>, itp., i także propaguje się w
| |
| dalszych obliczeniach.
| |
|
| |
| Działania arytmetyczne nie są łączne, co widać na poniższym przykładzie:
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| octave:9> 7.1 - (7+0.1)
| |
| ans <nowiki>=</nowiki> 0
| |
| octave:10> (7.1 - 7) - 0.1
| |
| ans <nowiki>=</nowiki> -3.6082e-16
| |
| </pre></div>
| |
|
| |
| {{cwiczenie|||
| |
| Wyjaśnij, dlaczego w arytmetyce podwójnej precyzji IEEE 754 mamy
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| octave:19> 2006/1e309
| |
| ans <nowiki>=</nowiki> 0
| |
| octave:20> 2.006/1e306
| |
| ans <nowiki>=</nowiki> 2.0060e-306
| |
| </pre></div>
| |
|
| |
| }}
| |
|
| |
| Wbrew pozorom, fakt, że nie mamy dostępu do arytmetyki nieskończonej precyzji
| |
| może mieć daleko idące konsekwencje, o czym przekonaliśmy się na początku
| |
| wykładu.
| |
|
| |
| ====Praktyczne wyznaczanie precyzji arytmetyki====
| |
|
| |
| Aby wyznaczyć precyzję używanej przez nas arytmetyki możemy wykonać prosty test.
| |
| Pomyślmy, jaka jest najmniejsza dodatnia liczba <math>\displaystyle \epsilon_{ \mbox{mach} }</math>, która dodana do jedności da w
| |
| wyniku liczbę ''większą'' od 1.0 (liczbę <math>\displaystyle \epsilon_{ \mbox{mach} }</math> nazywa się czasem epsilonem maszynowym, <tt>macheps</tt>).
| |
| Nietrudno sprawdzić, że liczba ta to <math>\displaystyle 2^{t}</math>, gdzie <math>\displaystyle t</math> to precyzja arytmetyki:
| |
|
| |
| Jasne jest, że w przypadku arytmetyki IEEE-754 jest to liczba równa podwojonej
| |
| precyzji arytmetyki, <math>\displaystyle 2^{-t}</math>, gdzie <math>\displaystyle t</math> jest liczbą cyfr mantysy <math>\displaystyle f</math>. Stąd
| |
| dostajemy prosty algorytm wyznaczania epsilona maszynowego:
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| x <nowiki>=</nowiki> 1.0;
| |
| while ( 1.0 + x > 1.0 )
| |
| {
| |
| x <nowiki>=</nowiki> x / 2.0;
| |
| }
| |
| printf("Macheps <nowiki>=</nowiki> \%g", 2.0*x);
| |
| }
| |
| </pre></div>
| |
|
| |
| Jednak, w rzeczywistości musimy być bardziej ostrożni. Implementując ten
| |
| algorytm w C następująco (pełna implementacja w pliku <code>macheps.c</code>)
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| /* Wyznaczanie epsilona maszynowego, wersja 1 */
| |
| #include <stdio.h>
| |
|
| |
| int main(void)
| |
| {
| |
| int dt;
| |
| double dx;
| |
|
| |
| dt <nowiki>=</nowiki> 0; dx <nowiki>=</nowiki> 1.0;
| |
| while(1.0 + dx > 1.0)
| |
| {
| |
| dx *<nowiki>=</nowiki> 0.5;
| |
| dt++;
| |
| }
| |
| printf("Macheps (double) <nowiki>=</nowiki> \%g. Liczba bitów mantysy <nowiki>=</nowiki> \%d\n", 2*dx, dt);
| |
| return(0);
| |
| }
| |
| </pre></div>
| |
|
| |
| dostajemy wynik niezgodny z oczekiwaniami:
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| Macheps <nowiki>=</nowiki> 1.0842e-19. Liczba bitów mantysy <nowiki>=</nowiki> 64.
| |
| </pre></div>
| |
|
| |
| Wynika to stąd, że w C obliczenia wykonują się zawsze z maksymalną możliwą
| |
| precyzją. W procesorach x86 jest to precyzja arytmetyki ''extended double
| |
| precision'',
| |
| wykorzystującej 80 bitów do reprezentacji liczb. Dlatego działanie
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| 1.0 + dx > 1.0
| |
| </pre></div>
| |
|
| |
| wykona się w arytmetyce nie podwójnej (64-bitowej), ale ''rozszerzonej
| |
| podwójnej'' precyzji. Aby sprawić, by działanie zostało wykonane z wykorzystaniem
| |
| typu <code>double</code>, musimy nasz program trochę zmodyfikować:
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| /* Wyznaczanie epsilona maszynowego, wersja 2 */
| |
| #include <stdio.h>
| |
|
| |
| int main(void)
| |
| {
| |
| int dt;
| |
| double dx, dxp1;
| |
|
| |
| dt <nowiki>=</nowiki> 0; dx <nowiki>=</nowiki> 1.0; dxp1 <nowiki>=</nowiki> 2.0;
| |
| while(dxp1 > 1.0)
| |
| {
| |
| dx *<nowiki>=</nowiki> 0.5;
| |
| dxp1 <nowiki>=</nowiki> 1.0 + dx; /* tym razem wynik działania zostanie zapisany
| |
| do zmiennej typu double */
| |
| dt++;
| |
| }
| |
| printf("Macheps (double) <nowiki>=</nowiki> \%g. Liczba bitów mantysy <nowiki>=</nowiki> \%d\n", 2*dx, dt);
| |
| }
| |
| </pre></div>
| |
|
| |
| Tym razem wynik jest prawidłowy:
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| Macheps <nowiki>=</nowiki> 2.22045e-16. Liczba bitów mantysy <nowiki>=</nowiki> 53
| |
| </pre></div>
| |
|
| |
| {{cwiczenie|||
| |
| Sprawdź, jak zmienią się wyniki, gdy wykorzystasz w swoim programie (zarówno w
| |
| wersji 1, jak w wersji 2) opcje kompilacji:
| |
| * <code>gcc -O3</code>
| |
| * <code>gcc -ffast-math</code>
| |
| * <code>gcc -O3 -ffast-math</code>
| |
|
| |
| Spróbuj objaśnić te wyniki, wspomagając się ewentualnie dokumentacją
| |
| kompilatora.
| |
| }}
| |
|
| |
| LAPACK daje gotową funkcję, <code>DLAMCH</code> (dla liczb podwójnej precyzji) i
| |
| <code>SLAMCH</code> (dla pojedynczej precyzji), pozwalającą stwierdzić
| |
| eksperymentalnie, jakie są parametry używanej arytmetyki, m.in. zakres liczb
| |
| reprezentowalnych, w jakim systemie reprezentowana jest mantysa, oraz oczywiście
| |
| precyzję arytmetyki i liczbę cyfr mantysy. Polecam analizę kodu źródłowego
| |
| <code>LAPACK/dlamch1.f</code> oraz lekturę prac
| |
| * Malcolm M. A. (1972) ''Algorithms to reveal properties of
| |
| floating-point arithmetic.'' Comms. of the ACM, 15, 949-951.
| |
| * Gentleman W. M. and Marovich S. B. (1974) ''More on algorithms
| |
| that reveal properties of floating point arithmetic units.''
| |
| Comms. of the ACM, 17, 276-277.
| |
|
| |
| na których oparto tę funkcję LAPACKa. Poniżej przykład zastosowania tej funkcji
| |
| (kod źródłowy: <code>lamch.c</code>) i wyniki uzyskane na procesorze x86.
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| /* Wyznaczanie epsilona maszynowego i innych charakterystyk arytmetyki podwójnej
| |
| precyzji z wykorzystaniem funkcji DLAMCH z LAPACKa */
| |
| #include <stdio.h>
| |
| #include <math.h>
| |
|
| |
| double dlamch_(char *CMACH); /* funkcja DLAMCH z LAPACKa */
| |
|
| |
| int main(void)
| |
| {
| |
| char CMACH;
| |
|
| |
| CMACH <nowiki>=</nowiki> 'e';
| |
| printf("Epsilon maszynowy: \%g\n", dlamch_(&CMACH));
| |
| CMACH <nowiki>=</nowiki> 'b';
| |
| printf("Podstawa arytmetyki: \%g\n", dlamch_(&CMACH));
| |
| CMACH <nowiki>=</nowiki> 'n';
| |
| printf("Liczba bitów mantysy: \%g\n", dlamch_(&CMACH));
| |
| CMACH <nowiki>=</nowiki> 'u';
| |
| printf("Zakres: \%g ", dlamch_(&CMACH));
| |
| CMACH <nowiki>=</nowiki> 'o';
| |
| printf("... \%g\n", dlamch_(&CMACH));
| |
| CMACH <nowiki>=</nowiki> 'r';
| |
| if(dlamch_(&CMACH) <nowiki>=</nowiki><nowiki>=</nowiki> 1.0)
| |
| printf("Z zaokrąglaniem (dlamch('r') <nowiki>=</nowiki> \%g)\n", dlamch_(&CMACH));
| |
| else
| |
| printf("Bez zaokrąglania (dlamch('r') <nowiki>=</nowiki> \%g)\n", dlamch_(&CMACH));
| |
| return(0);
| |
| }
| |
| </pre></div>
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| Epsilon maszynowy: 2.22045e-16
| |
| Podstawa arytmetyki: 2
| |
| Liczba bitów mantysy: 53
| |
| Zakres: 2.22507e-308 ... 1.79769e+308
| |
| Bez zaokrąglania (dlamch('r') <nowiki>=</nowiki> 0)
| |
| </pre></div>
| |
|
| |
| ===Wpływ błędu zaokrągleń na wyniki obliczeń. Redukcja cyfr i inne
| |
| patologie===
| |
|
| |
| Korzystając z wprowadzonego powyżej modelu arytmetyki zmiennoprzecinkowej możemy
| |
| spróbować uchwycić wpływ błędu reprezentacji i błędu zaokrągleń na wynik
| |
| konkretnego algorytmu.
| |
|
| |
| {{przyklad|||
| |
| Rozważmy banalne zadanie wyznaczenia iloczynu <math>\displaystyle N</math> liczb z tablicy <math>\displaystyle x</math>,
| |
| <center><math>\displaystyle
| |
| s = x_0\cdot \cdots \cdot x_{N-1}.
| |
| </math></center>
| |
|
| |
| W tym celu stosujemy banalny algorytm:
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| s <nowiki>=</nowiki> 1.0;
| |
| for (i<nowiki>=</nowiki>0; i < N; i++)
| |
| s *<nowiki>=</nowiki> x[i];
| |
| </pre></div>
| |
|
| |
| Sprawdźmy, jak będzie on realizowany w arytmetyce zmiennoprzecinkowej. Dla
| |
| uproszczenia założymy, że nie wystąpiło ani zjawisko nadmiaru, ani zjawisko
| |
| niedomiaru (w przeciwnym razie dostaniemy w wyniku, odpowiednio, <math>\displaystyle \pm</math><code>Inf</code>
| |
| lub 0).
| |
|
| |
| Naturalnie, zamiast dokładnych wartości <math>\displaystyle x_0, \ldots x_{N-1}</math>, będziemy mieli w
| |
| komputerze jedynie ich reprezentacje, <math>\displaystyle \widetilde{x}_i = rd_\nu(x_i) = x_i ( 1 +
| |
| \delta_i)</math>, przy czym <math>\displaystyle |\delta_i| \leq \nu</math>.
| |
|
| |
| Oznaczając <math>\displaystyle \widetilde{s}_i</math> wyznaczoną numerycznie wartość iloczynu po <math>\displaystyle i</math>-tym
| |
| kroku pętli, mamy, że
| |
| <center><math>\displaystyle
| |
| \widetilde{s}_{i+1} = fl_\nu(\widetilde{s}_i \times \widetilde{x}_i) = \widetilde{s}_i \cdot
| |
| \widetilde{x}_i \cdot (1 + \epsilon_i),
| |
| </math></center>
| |
|
| |
| gdzie znów <math>\displaystyle |\epsilon_i| \leq \nu</math>. Ostatecznie więc, wyznaczona wartość
| |
| iloczynu, <math>\displaystyle \widetilde{s}</math> spełnia
| |
| <center><math>\displaystyle
| |
| \widetilde{s} = x_0\cdots x_{N-1} \cdot \Pi_{i=0}^{N-1}(1+\epsilon_i)(1+\delta_i).
| |
| </math></center>
| |
|
| |
| Ponieważ <math>\displaystyle \Pi_{i=0}^{N-1}(1+\epsilon_i) =
| |
| (1 + {\cal E})</math>, gdzie, z pominięciem małych wyższego rzędu, <math>\displaystyle |{\cal E}| \leq N\nu</math>, dostajemy
| |
| ostatecznie
| |
| <center><math>\displaystyle
| |
| \widetilde{s} = s \cdot (1+E),
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle |E|\leq 2N\nu</math>. Ponieważ w arytmetyce podwójnej precyzji <math>\displaystyle \nu \approx
| |
| 10^{-16}</math>, to nawet biorąc iloczyn tysiąca (!) liczb, dostajemy wynik, którego
| |
| błąd względny będzie (o ile tylko nie wystąpi nadmiar/niedomiar) bardzo mały,
| |
| rzędu <math>\displaystyle 10^{-13}</math>!
| |
| }}
| |
|
| |
| Powyższe rozumowanie, a także intuicja często wyrażana przez osoby postronne,
| |
| prowadzi do przypuszczenia, że:
| |
|
| |
| <blockquote style="background-color:#fefeee">
| |
| "Duży błąd względny wyniku jest możliwy dopiero po ''kumulacji'' błędów
| |
| zaokrągleń po przeprowadzeniu ''bardzo wielu'' działań arytmetycznych."
| |
| </blockquote>
| |
|
| |
| Jednak to jest to całkowicie '''fałszywy''' pogląd, o czym świadczy kolejny,
| |
| bardzo znamienny
| |
| przykład.
| |
|
| |
| {{przyklad|Redukcja cyfr przy odejmowaniu||
| |
|
| |
| Tym razem nasze zadanie jest znacznie prostsze od poprzedniego. Trzeba wyznaczyć
| |
| po prostu różnicę dwóch liczb:
| |
| <center><math>\displaystyle
| |
| s = a - b.
| |
| </math></center>
| |
|
| |
| Prowadząc analizę jak poprzednio, widzimy, że obliczona numerycznie wartość to
| |
| <center><math>\displaystyle
| |
| \widetilde{s} = fl_\nu(rd_\nu(a) - rd_\nu{b}) = (a(1+\delta_a) - b(1+\delta_b))(1+\epsilon),
| |
| </math></center>
| |
|
| |
| Stąd po prostych oszacowaniach
| |
|
| |
| <center><math>\displaystyle
| |
| \left|\frac{\widetilde{s} - s}{s}\right| \leq 2\frac{|a| + |b|}{|a-b|} \cdot \nu.
| |
| </math></center>
| |
|
| |
| A więc, gdy <math>\displaystyle a\approx b</math>, to <math>\displaystyle \frac{|a| + |b|}{|a-b|} \approx \infty</math> i w
| |
| efekcie możemy utracić nawet ''wszystkie'' znaczące cyfry wyniku! To zjawisko
| |
| właśnie nosi żargonową nazwę "'''utraty cyfr przy odejmowaniu'''", choć
| |
| precyzyjnie powinno się mówić o "zmniejszeniu liczby dokładnych cyfr znaczących
| |
| wyniku przy odejmowaniu dwóch bardzo bliskich sobie liczb".
| |
|
| |
| Przy okazji zauważmy, że prowadząc identyczną analizę dla ''sumy'' dwóch liczb
| |
| <math>\displaystyle a+b</math>, gdzie <math>\displaystyle a</math> i <math>\displaystyle b</math> są
| |
| ''tego samego'' znaku, dostajemy oszacowanie błędu względnego równe <math>\displaystyle 2\nu</math>,
| |
| niezależnie od wartości liczbowych <math>\displaystyle a</math> i <math>\displaystyle b</math>!
| |
|
| |
| }}
| |
|
| |
| Skutki zjawiska redukcji cyfr przy odejmowaniu mogą być dramatyczne i ujawnić
| |
| się w nadspodziewanie elementarnych sytuacjach.
| |
|
| |
| {{przyklad|Wyznaczanie pierwiastków trójmianu kwadratowego||
| |
|
| |
| Niech <math>\displaystyle a,p,q>0</math>. Korzystając ze szkolenego wzoru na pierwiastki równania kwadratowego <math>\displaystyle ax^2 -
| |
| 2px + q = 0</math>
| |
| <center><math>\displaystyle
| |
| x_{1,2} = \frac{1}{a} (p \pm \sqrt{\Delta}),
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \Delta = p^2 - qa > 0</math>, możemy natknąć się na trudności, gdy jeden z
| |
| pierwiastków jest bardzo bliski zera (tzn. gdy <math>\displaystyle p \approx \sqrt{\Delta}</math>).
| |
|
| |
| Taka
| |
| sytuacja np. powstaje gdy rozważać jednostajnie opóźniony ruch pocisku
| |
| wystrzelonego z działka przeciwlotniczego do celu lecącego na małej wysokości. Czas
| |
| trafienia w cel jest --- przy pominięciu oporu powietrza --- rozwiązaniem
| |
| równania kwadratowego, przy czym czs krótki odpowiada bezpośredniemu trafieniu w
| |
| cel, a czas długi --- wystrzeleniu pocisku z wyprzedzeniem wysoko w górę i
| |
| poczekaniu, aż spadając trafi w cel od góry. Oczywiście, żaden artylerzysta nie
| |
| będzie się interesował tym ostatnim przypadkiem: interesujące jest więc ''dokładne'' (bo cel leci szybko) wyznaczenie ''mniejszego'' pierwiastka.
| |
|
| |
| Niestety, skoro <math>\displaystyle p \approx \sqrt{\Delta}</math>, to wyznaczając mniejszy pierwiastek
| |
| <math>\displaystyle x_1</math>
| |
| ryzykujemy utratę cyfr przy odejmowaniu. Ale na szczęście, można to ominąć
| |
| zgodnie z często stosowaną w numeryce regułą:
| |
|
| |
| <blockquote style="background-color:#fefeee">
| |
| Jeśli problem jest trudny --- najlepiej zmienić problem!
| |
| </blockquote>
| |
|
| |
| W naszym wypadku ratunkiem jest matematyczna transformacja problemu tak, by już
| |
| nie było w nim odejmowania bliskich sobie liczb. Rzeczywiście, przecież wciąż
| |
| mamy dobry wzór na ''większy'' z pierwiastków, <math>\displaystyle x_2 = \frac{1}{a} (p +
| |
| \sqrt{\Delta})</math>! Dokładając do tego wzór Viete'a,
| |
|
| |
| <center><math>\displaystyle
| |
| x_1 x_2 = \frac{q}{a},
| |
| </math></center>
| |
|
| |
| dostajemy inny wzór na <math>\displaystyle x_1</math>, nie zawierający feralnego dzielenia. Poniżej
| |
| demonstrujemy program w C, testujący jakość obu podejść.
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
| # include <stdio.h>
| |
| # include <math.h>
| |
|
| |
| /* w(x) <nowiki>=</nowiki> ax^2 - 2px + q <nowiki>=</nowiki> 0 */
| |
| /* delta <nowiki>=</nowiki> 4(p^2 - qa) */
| |
|
| |
| double const a <nowiki>=</nowiki> 2.1, q <nowiki>=</nowiki> 1e-6, p<nowiki>=</nowiki>1.1;
| |
|
| |
| double w(double x) /* wartość wielomianu w punkcie x */
| |
| {
| |
| return(a*x*x - 2.0*p*x + q);
| |
| }
| |
|
| |
| int main(void)
| |
| {
| |
| double x1, x2, x1v, X1, X1v, X2;
| |
| double Delta; /* wartość Delty liczymy w podwójnej precyzji */
| |
| float delta; /* wartość delty liczymy w pojedynczej precyzji */
| |
|
| |
| delta <nowiki>=</nowiki> Delta <nowiki>=</nowiki> sqrt(p*p - q*a);
| |
| printf("Wielomian w(x) <nowiki>=</nowiki> \%e x^2 - \%e x + \%e.\nDelta <nowiki>=</nowiki> \%e\n", a, 2*p, q, delta);
| |
|
| |
| /* pierwiastki liczone wzorem szkolnym, z niedokładną deltą */
| |
| x1 <nowiki>=</nowiki> (p - delta)/a;
| |
| x2 <nowiki>=</nowiki> (p + delta)/a;
| |
|
| |
| /* mniejszy pierwiatek, liczony z mało dokładną deltą, ale lepszym
| |
| wzorem: Viete'a */
| |
| x1v <nowiki>=</nowiki> (q/a)/x2;
| |
|
| |
| /* pierwiastki liczone wzorem szkolnym, z dokładniejszą Deltą */
| |
| X1 <nowiki>=</nowiki> (p - Delta)/a;
| |
| X2 <nowiki>=</nowiki> (p + Delta)/a;
| |
|
| |
| /* mniejszy pierwiatek, liczony z dokładniejszą Deltą, ale lepszym
| |
| wzorem: Viete'a */
| |
| X1v <nowiki>=</nowiki> (q/a)/X2;
| |
|
| |
| printf("\nPierwiastki z mało dokładną deltą:\n");
| |
| printf(" Wzór szkolny: x1 <nowiki>=</nowiki> \%e x2 <nowiki>=</nowiki> \%e\n Wzór Viete'a: x1v <nowiki>=</nowiki> \%e x2 <nowiki>=</nowiki> j.w.\n",
| |
| x1,x2,x1v);
| |
| printf("\nPierwiastki z dokładniejszą Deltą:\n");
| |
| printf(" Wzór szkolny: X1 <nowiki>=</nowiki> \%e X2 <nowiki>=</nowiki> \%e\n Wzór Viete'a: X1v <nowiki>=</nowiki> \%e X2 <nowiki>=</nowiki> j.w.\n",
| |
| X1,X2,X1v);
| |
| printf("\nWzględna zmiana wartości pierwiastka:\n");
| |
| printf(" (x1 - x1v)/x1v <nowiki>=</nowiki> \%e\n", (x1-x1v)/x1v);
| |
| printf(" (x1v -X1v)/X1v <nowiki>=</nowiki> \%e\n", (x1v-X1v)/X1v);
| |
| printf(" (x2 - X2)/X2 <nowiki>=</nowiki> \%e\n", (x2-X2)/X2);
| |
|
| |
| printf("\nWartość wielomianu w wyznaczonych punktach:\n");
| |
| printf(" w(x1) <nowiki>=</nowiki> \%e\n w(x1v) <nowiki>=</nowiki> \%e w(X1v) <nowiki>=</nowiki> \%e\n w(x2) <nowiki>=</nowiki> \%e\n w(X2) <nowiki>=</nowiki> \%e\n ",
| |
| w(x1),w(x1v),w(X1v),w(x2),w(X2));
| |
|
| |
| return(0);
| |
| }
| |
| </pre></div>
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| Wielomian w(x) <nowiki>=</nowiki> 2.100000e+00 x^2 - 2.200000e+00 x + 1.000000e-06.
| |
| Delta <nowiki>=</nowiki> 1.099999e+00
| |
|
| |
| Pierwiastki z mało dokładną deltą:
| |
| Wzór szkolny: x1 <nowiki>=</nowiki> 4.427774e-07 x2 <nowiki>=</nowiki> 1.047619e+00
| |
| Wzór Viete'a: x1v <nowiki>=</nowiki> 4.545456e-07 x2 <nowiki>=</nowiki> j.w.
| |
|
| |
| Pierwiastki z dokładniejszą Deltą:
| |
| Wzór szkolny: X1 <nowiki>=</nowiki> 4.545457e-07 X2 <nowiki>=</nowiki> 1.047619e+00
| |
| Wzór Viete'a: X1v <nowiki>=</nowiki> 4.545457e-07 X2 <nowiki>=</nowiki> j.w.
| |
|
| |
| Względna zmiana wartości pierwiastka:
| |
| (x1 - x1v) / x1v <nowiki>=</nowiki> -2.589022e-02
| |
| (x1v -X1v) / X1v <nowiki>=</nowiki> -1.123337e-08
| |
| (x2 - X2) / X2 <nowiki>=</nowiki> 1.123337e-08
| |
|
| |
| Wartość wielomianu w wyznaczonych punktach:
| |
| w(x1) <nowiki>=</nowiki> 2.589022e-08
| |
| w(x1v) <nowiki>=</nowiki> 1.123337e-14, w(X1v) <nowiki>=</nowiki> -3.194985e-23
| |
| w(x2) <nowiki>=</nowiki> 2.589022e-08
| |
| w(X2) <nowiki>=</nowiki> -1.357688e-17
| |
| </pre></div>
| |
|
| |
| Jak więc widzimy, nawet z niezbyt dokładnie wyznaczoną deltą, mniejszy
| |
| pierwiastek jesteśmy w stanie wyznaczyć bardzo precyzyjnie --- o ile tylko
| |
| unikniemy redukcji cyfr przy odejmowaniu.
| |
|
| |
| }}
| |
|
| |
| {{przyklad|Numeryczna "lupa"||
| |
|
| |
| Ten obrazek już widzieliśmy na początku wykładu.
| |
|
| |
| [[Image:MNwielomian4.png|thumb|300px|Wykres funkcji <math>\displaystyle w(x) = x^4 - 4x^3+6x^2-4x+1 =
| |
| (x-1)^4</math> wyznaczony na dwa sposoby w arytmetyce podwójnej
| |
| precyzji.]]
| |
|
| |
| Wyjaśnieniem tej niepokojącej obserwacji jest znowu zjawisko redukcji cyfr
| |
| przy odejmowaniu: wartości <math>\displaystyle f</math> są bliskie zera, a powstają jako suma dużych
| |
| liczb z przeciwnymi znakami.
| |
| }}
| |
|
| |
| Czasem możemy spotkać się z twierdzeniem, że nie warto zawracać sobie głowy
| |
| metodami numerycznymi, gdyż są dostępne pakiety obliczeń symbolicznych (Maple,
| |
| Mathematica, MuPAD, Maxima), które potrafią "wszystko policzyć z dowolną
| |
| precyzją".
| |
|
| |
| To oczywiście też nie jest do końca prawdą. Precyzja, o której mowa, jest
| |
| jedynie precyzją używanej arytmetyki (rzeczywiście, softwarowo można emulować
| |
| dowolną precyzję), ale dokładność ''wyniku'' nie może być w nich a priori
| |
| zadana. Wiele bowiem może zależeć od właściwości samego zadania obliczeniowego,
| |
| o czym mówimy w następnej części wykładu. Na razie prosty przykład:
| |
|
| |
| {{przyklad|Co oznacza precyzja w pakietach symbolicznych||
| |
|
| |
| Oto fragment sesji pakietu obliczeń symbolicznych MUPAD:
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| >> ((4/3)*3 - 3) - 1
| |
|
| |
| 0
| |
| >> DIGITS :<nowiki>=</nowiki> 10
| |
|
| |
| 10
| |
| >> ((4/3.0)*3 - 3) - 1
| |
|
| |
| -2.168404345e-19
| |
| >> subs(((4/a)*3 - 3) - 1, a <nowiki>=</nowiki> 3.0)
| |
|
| |
| -4.33680869e-19
| |
| >> subs(((4/a)*3 - 3) - 1, a <nowiki>=</nowiki> 3)
| |
|
| |
| 0
| |
| </pre></div>
| |
|
| |
| Najpierw wyznaczona została wartość wyrażenia algebraicznego, korzystając z
| |
| manipulacji symbolicznej --- oczywiście bez trudu system stwierdził, że to
| |
| wyrażenie upraszcza się do zera.
| |
|
| |
| Następnie zażądaliśmy, by <code>DIGITS</code> --- parametr sterujący "liczbą cyfr
| |
| znaczących dla liczb zmiennoprzecinkowych", jak to określa manual MUPADa ---
| |
| przyjął wartość równą 10.
| |
|
| |
| Wymuszając, przez podanie <code>3.0</code> zamiast <code>3</code> stosowanie w
| |
| obliczeniach arytmetyki
| |
| zmiennoprzecinkowej zamiast symbolicznej (pamiętasz, jak to jest w C?) dostajemy
| |
| wynik, który nie ma ani jednej cyfry znaczącej dokładnej. Z drugiej strony,
| |
| widzimy także, iż faktycznie stosowana precyzja obliczeń jest znacznie większa i wynosi
| |
| około 19 cyfr znaczących, chociaż żądaliśmy jedynie 10...
| |
|
| |
| }}
| |
|
| |
| jak wynika z powyższego, w praktyce pakiety symboliczne
| |
| stosują znacznie większą niż żądana precyzję obliczeń, by ustrzec się
| |
| najbardziej typowych patologii.
| |
|
| |
| ===Uwarunkowanie zadania obliczeniowego===
| |
|
| |
| Jak zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując
| |
| zadanie obliczeniowe są z natury rzeczy wartościami zaburzonymi. Źródłem tych
| |
| zaburzeń są:
| |
| * błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (0.1 nie jest
| |
| równe dokładnie <math>\displaystyle 1/10</math>)
| |
| * błędy w parametrach będących rezultatem poprzednich obliczeń (chcemy
| |
| rozwiązać równanie <math>\displaystyle f(x) = a</math>, ale <math>\displaystyle a</math> jest rezultatem innej symulacji), a także
| |
| * błędy pomiarowe w zadaniach praktycznych (chcemy policzyć
| |
| numeryczną prognozę pogody, ale temperaturę wyjściową znamy tylko z dokładnością
| |
| do 0.1 stopnia, co gorsza --- z niektórych stacji w ogóle nie mamy danych)
| |
|
| |
| Okazuje się, że powszechna intuicja, że małe zaburzenia danych powinny dawać
| |
| małe zaburzenia wyniku, nie znajduje potwierdzenia nawet w bardzo prostych
| |
| przypadkach. Z drugiej strony, umiejętność oceny jakościowego ''wpływu
| |
| zaburzenia danych na wynik'' jest kapitalna w świecie obliczeń numerycznych w
| |
| ogólności, a w szególności --- inżynierskich.
| |
|
| |
| Wprowadza się pojęcie ''uwarunkowania'' zadania, to znaczy jego podatności na
| |
| zaburzenia danych. Dla przejrzystości, przypuśćmy, że nasze zadanie obliczeniowe
| |
| polega na wyznaczeniu <math>\displaystyle f(x)</math> dla danego <math>\displaystyle x</math>.
| |
|
| |
| <div class="thumb tright"><div><flash>file=XXX.png</flash><div.thumbcaption>Zadanie obliczeniowe i jego odporność na zaburzenia</div></div></div>
| |
|
| |
| \beginnowiki
| |
| [[Image:MNcondition.png|thumb|300px|Naszym zadaniem jest wyznaczenie, dla <math>\displaystyle x\in X</math>, wartości
| |
| <math>\displaystyle f(x)\in Y</math>.]]
| |
| [[Image:MNcondition2.png|thumb|300px|Jaki będzie rozrzut wyników, gdy ''lekko'' zaburzymy
| |
| dane?]]
| |
| [[Image:MNcondition3.png|thumb|300px|Jeśli równie mały, co zaburzenie, powiemy, że zadanie
| |
| jest dobrze uwarunkowane (jego wynik jest mało podatny na zaburzenia danych).]]
| |
| [[Image:MNcondition4.png|thumb|300px|Może jednak zdarzyć się, że zadanie jest źle
| |
| uwarunkowane, i małe zaburzenie danych skutkuje dużym rozrzutem wyników.]]
| |
| [[Image:MNcondition5.png|thumb|300px|Wtedy nawet bliskie sobie punkty w X, przekształcenie
| |
| <math>\displaystyle f</math> może odwzorowywać w punkty bardzo od siebie odległe. Jest to sytuacja
| |
| skrajnie niekorzystna w zastosowaniach, a zwłaszcza --- w obliczeniach numerycznych.]]
| |
| \endnowiki
| |
|
| |
| Jak bardzo będzie odległe
| |
| <math>\displaystyle f(\widetilde{x})</math>, gdy <math>\displaystyle \widetilde{x}\approx x</math>? Rozważa się dwa przypadki:
| |
| * uwarunkowanie względne: jak względne zaburzenie danych wpływa na błąd
| |
| względny wyniku:
| |
|
| |
| <center><math>\displaystyle
| |
| \frac{||f(x) - f(\widetilde{x})||}{||f(x)||} \leq \mbox{cond} _{rel}(f,x) \cdot \frac{||x - \widetilde{x}||}{||x||}
| |
| </math></center>
| |
|
| |
| Najmniejszy mnożnik <math>\displaystyle \mbox{cond} _{rel}(f,x)</math> spełniający powyższą nierówność
| |
| nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia <math>\displaystyle f(x)</math>
| |
| dla danego <math>\displaystyle x</math>.
| |
| * uwarunkowanie bezwzględne: jak bezwzględne zaburzenie danych wpływa na błąd
| |
| bezwzględny wyniku:
| |
|
| |
| <center><math>\displaystyle
| |
| ||f(x) - f(\widetilde{x})|| \leq \mbox{cond} _{abs}(f,x) \cdot ||x - \widetilde{x}||
| |
| </math></center>
| |
|
| |
| Najmniejszy mnożnik <math>\displaystyle \mbox{cond} _{abs}(f,x)</math> spełniający powyższą nierówność
| |
| nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia <math>\displaystyle f(x)</math>
| |
| dla danego <math>\displaystyle x</math>.
| |
|
| |
| Powiemy, że zadanie jest
| |
| * dobrze uwarunkowane w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle \mbox{cond} (f,x) \approx 1</math>,
| |
| * źle uwarunkowane w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle \mbox{cond} (f,x) >> 1</math>,
| |
| * źle postawione w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle \mbox{cond} (f,x) = +\infty</math>.
| |
|
| |
| Teraz już rozumiemy, że redukcja cyfr przy odejmowaniu jest tylko
| |
| odzwierciedleniem faktu, że zadanie odejmowania dwóch bliskich liczb jest po
| |
| prostu zadaniem źle uwarunkowanym!
| |
|
| |
| {{przyklad|Uwarunkowanie zadania obliczenia sumy||
| |
|
| |
| Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia <math>\displaystyle s(x,y) = x + y</math> ma
| |
| <center><math>\displaystyle
| |
| \mbox{cond} _{abs}(s, (a,b)) = 1, \qquad \mbox{cond} _{rel}(s, (a,b)) = \frac{|a|+|b|}{|a+b|}
| |
| </math></center>
| |
|
| |
| Tak więc, gdy <math>\displaystyle a\approx -b</math>, to <math>\displaystyle \mbox{cond} _{rel}(s, (a,b)) \approx +\infty</math> i zadanie
| |
| jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może
| |
| skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego,
| |
| najczęściej rzeczywiście tak będzie...
| |
| }}
| |
|
| |
| {{przyklad|||
| |
| Dla różniczkowalnej funkcji skalarnej <math>\displaystyle f : R \rightarrow R</math> mamy
| |
| <center><math>\displaystyle
| |
| |f(x) - f(\widetilde{x}| \approx |f'(x) | | x - \widetilde{x} |
| |
| </math></center>
| |
|
| |
| i w konsekwencji dla zadania obliczenia <math>\displaystyle f(x)</math> dla danego <math>\displaystyle x</math> mamy, przy
| |
| założeniu małych zaburzeń,
| |
| <center><math>\displaystyle
| |
| \mbox{cond} _{abs}( f, x) = |f'(x)|, \qquad \mbox{cond} _{rel}( f, x) =
| |
| \frac{|f'(x)|\cdot|x|}{|f(x)|}.
| |
| </math></center>
| |
|
| |
| }}
| |
|
| |
| Możnaby myśleć, że złe uwarunkowanie zawsze jest szkodliwe w praktyce
| |
| numerycznej. Najczęściej właśnie tak jest istotnie. Jednak w praktyce
| |
| numerycznej sporadycznie zdarza się, że [[sec:invit|Uzupe�nij: złe uwarunkowanie pewnego podzadania nie tylko
| |
| nie pogarsza sytuacji, ale wręcz pomaga]] szybciej rozwiązać zadanie główne!
| |
|
| |
| ===Rozkład algorytmu względem informacji===
| |
|
| |
| ''Algorytm'' to dokładnie określona i dozwolona w danym modelu
| |
| obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego
| |
| zadania (w sposób dokładny lub przybliżony).
| |
|
| |
| Z każdym algorytmem związany jest operator
| |
|
| |
| <center><math>\displaystyle {\bf ALG}:\,F\longrightarrowG,
| |
| </math></center>
| |
|
| |
| taki że <math>\displaystyle {\bf ALG}(f)</math> jest wynikiem działania algorytmu
| |
| w arytmetyce idealnej dla danej <math>\displaystyle f</math>.
| |
|
| |
| Zauważmy, że wynik <math>\displaystyle {\bf ALG}(f)</math> działania algorytmu nie
| |
| zależy bezpośrednio od <math>\displaystyle f</math>, ale raczej od ''informacji''
| |
| o <math>\displaystyle f</math> (uzyskanej dzięki poleceniu <math>\displaystyle {\cal IN}</math>). Informacja
| |
| ta może być ''pełna'' albo tylko ''częściowa''.
| |
| Informacja jest pełna gdy, np.
| |
| <math>\displaystyle f=(f_1,\ldots,f_n)\inR^n</math> i wczytamy wszystkie
| |
| współrzędne <math>\displaystyle f_i</math>. Informacja może być częściowa, gdy
| |
| <math>\displaystyle f</math> jest funkcją. Wtedy wiele danych może posiadać tę
| |
| samą informację, co łatwo zaobserwować na przykładzie
| |
| zadania całkowania.
| |
|
| |
| Niech <math>\displaystyle N:F\to\cup_{n=0}^\inftyR^n</math> będzie
| |
| ''operatorem informacji'', tzn.
| |
|
| |
| <center><math>\displaystyle N(f)\,=\,(y_1,y_2,\ldots,y_n)
| |
| </math></center>
| |
|
| |
| jest informacją o <math>\displaystyle f</math> zebraną przy idealnej realizacji
| |
| algorytmu. Zauważmy, że nformacja jest pełna gdy <math>\displaystyle N</math> jest
| |
| przekształceniem różnowartościowym, tzn. jeśli
| |
| <math>\displaystyle f_1\nef_2</math> implikuje <math>\displaystyle N(f_1)\neN(f_2)</math>.
| |
| W przeciwnym przypadku mamy do czynienia z informacją
| |
| częściową.
| |
|
| |
| Każdy algorytm <math>\displaystyle {\bf ALG}</math> może być przedstawiony jako złożenie
| |
| operatora informacji i pewnego operatora
| |
| <math>\displaystyle \varphi:N(F)\toG</math> zdefiniowanego równością
| |
|
| |
| <center><math>\displaystyle \varphi\left(N(f)\right)\,=\,{\bf ALG}(f).
| |
| </math></center>
| |
|
| |
| Zauważmy, że w przypadku informacji częściowej zwykle nie
| |
| istnieje algorytm dający dokładne rozwiązanie zadania dla
| |
| każdej danej <math>\displaystyle f\inF</math>, ponieważ dla danych o tej samej
| |
| informacji mogą istnieć różne rozwiązania.
| |
|
| |
| ===Problem wyboru algorytmu===
| |
|
| |
| Wybór algorytmu jest najistotniejszą częścią całego procesu
| |
| numerycznego rozwiązywania zadania. Kierujemy się przy tym przede
| |
| wszystkim następującymi kryteriami:
| |
| * dokładnością algorytmu,
| |
| * złożonością algorytmu,
| |
| * własnościami numerycznymi algorytmu.
| |
|
| |
| Przez dokładność algorytmu rozumiemy różnicę między
| |
| rozwiązaniem dokładnym <math>\displaystyle S(f)</math>, a rozwiązaniem
| |
| <math>\displaystyle {\bf ALG}(f)</math> dawanym przez algorytm w arytmetyce idealnej.
| |
| Jeśli <math>\displaystyle {\bf ALG}(f)=S(f)</math>, <math>\displaystyle \forallf\inF</math>,
| |
| to algorytm nazywamy ''dokładnym''.
| |
|
| |
| Mówiąc o złożoności, mamy na myśli złożoność pamięciową
| |
| (zwykle jest to liczba stałych i zmiennych używanych przez
| |
| algorytm), jak również złożoność obliczeniową.
| |
| Na złożoność obliczeniową algorytmu dla danej <math>\displaystyle f</math> składa
| |
| się koszt uzyskania infomacji <math>\displaystyle y=N(f)</math> (zwykle jest on
| |
| proporcjonalny do liczby wywołań polecenia <math>\displaystyle {\cal IN}</math>), oraz
| |
| koszt ''kombinatoryczny'' przetworzenia tej informacji, aż do
| |
| uzyskania wyniku <math>\displaystyle \varphi(y)</math>. Koszt kombinatoryczny zwykle
| |
| mierzymy liczbą operacji arytmetycznych wykonywanych przez
| |
| algorytm.
| |
|
| |
| Przez własności numeryczne algorytmu rozumiemy jego
| |
| własności przy realizacji w arytmetyce <math>\displaystyle fl_\nu</math>. Temu
| |
| ważnemu tematowi poświęcimy teraz osobny paragraf.
| |
|
| |
| ===Numeryczna poprawność algorytmu===
| |
|
| |
| Pożądane jest, aby algorytm dawał "dobry" wynik zarówno
| |
| w arytmetyce idealnej, jak i w arytmetyce <math>\displaystyle fl_\nu</math>. Niestety,
| |
| jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm
| |
| jest dokładny to w wyniku jego realizacji w <math>\displaystyle fl_\nu</math> możemy
| |
| otrzymać wynik <math>\displaystyle fl_\nu({\bf ALG}(f))</math> daleko odbiegający od
| |
| <math>\displaystyle S(f)</math>. W szczególności, prawie zawsze mamy
| |
|
| |
| <center><math>\displaystyle S(f)\,\ne\,fl_\nu\left({\bf ALG}(f)\right).
| |
| </math></center>
| |
|
| |
| Zauważmy również, że o ile z reguły znamy dokładne zachowanie
| |
| się algorytmu w arytmetyce idealnej dla danej informacji, to nie
| |
| można tego samego powiedzieć o jego zachowaniu się w arytmetyce
| |
| <math>\displaystyle fl_\nu</math>. W związku z tym powstaje pytanie, jak kontrolować błąd
| |
| algorytmów wynikający z błędów zaokrągleń i jakie algorytmy
| |
| uznamy za te o najwyższej jakości numerycznej.
| |
|
| |
| Istnienie błędów reprezentacji liczb rzeczywistych powoduje,
| |
| że informacja <math>\displaystyle y=N(f)</math> o danej <math>\displaystyle f</math> nie jest w
| |
| ogólności reprezentowana dokładnie. Znaczy to, że zamiast na
| |
| informacji dokładnej, dowolny algorytm będzie operować na
| |
| informacji ''nieco zaburzonej'' <math>\displaystyle y_\nu</math>, tzn. zaburzonej na
| |
| poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm
| |
| będzie w ogólności zaburzony na poziomie błędu reprezentacji.
| |
| W najlepszym więc wypadku wynikiem działania algorytmu w <math>\displaystyle fl_\nu</math>
| |
| będzie <math>\displaystyle (\varphi(y_\nu))_\nu</math> zamiast <math>\displaystyle \varphi(y)</math>. Algorytmy
| |
| dające tego rodzaju wyniki uznamy za posiadające najlepsze
| |
| własności numeryczne w arytmetyce <math>\displaystyle fl_\nu</math> i nazwiemy numerycznie
| |
| poprawnymi.
| |
|
| |
| Dokładniej, powiemy, że ciąg rzeczywisty
| |
| <math>\displaystyle a_\nu=(a_{\nu,1},\ldots,a_{\nu,n})</math>
| |
| (a właściwie rodzina ciągów <math>\displaystyle \{a_\nu\}_\nu</math>) jest
| |
| ''nieco zaburzonym'' ciągiem <math>\displaystyle a=(a_1,\ldots,a_n)</math>, jeśli
| |
| istnieje stała <math>\displaystyle K</math> taka, że dla wszystkich dostatecznie
| |
| małych <math>\displaystyle \nu</math> zachodzi
| |
|
| |
| <center><math>\displaystyle
| |
| |a_{\nu,j} - a_j|\,\le\,K\,\nu\,|a_j|,\qquad 1\le j\le n,
| |
| </math></center>
| |
|
| |
| albo ogólniej
| |
|
| |
| <center><math>\displaystyle
| |
| \|a_\nu - a\|\,\le\,K\,\nu\,\|a\|,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle \|\cdot\|</math> jest pewną normą w <math>\displaystyle R^n</math>. W pierwszym
| |
| przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim
| |
| o zaburzeniu w normie <math>\displaystyle \|\cdot\|</math>.
| |
|
| |
| Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają
| |
| za sobą niewielkie zaburzenia w normie. Rzeczywiście, jeśli
| |
| ([[##powsp|Uzupelnic: powsp ]]) to również
| |
|
| |
| <center><math>\displaystyle \|a_\nu - a\|_\infty \,=\, \max_{1\le j\le n} |a_{\nu,j} - a_j|
| |
| \,\le\,K\,\nu\,\max_{1\le j\le n} |a_j|\,=\,K\,\nu\,\|a\|_\infty,
| |
| </math></center>
| |
|
| |
| i korzystając z faktu, że w przestrzeni skończenie wymiarowej
| |
| wszystkie normy są równoważne otrzymujemy dla pewnych stałych
| |
| <math>\displaystyle K_1</math> i <math>\displaystyle K_2</math>
| |
|
| |
| <center><math>\displaystyle \|a_\nu - a\|\,\le\,K_1\|a_\nu-a\|_\infty\,\le\,
| |
| K_1 K\,\nu\,\|a\|_\infty\,\le\,K_2 K_1 K\,\nu\,\|a\|,
| |
| </math></center>
| |
|
| |
| czyli nierówność ([[##wnorm|Uzupelnic: wnorm ]]) ze stałą <math>\displaystyle K_2 K_1 K</math> zamiast <math>\displaystyle K</math>.
| |
|
| |
| {{definicja|Algorytm <math>\displaystyle {\bf ALG}</math> rozwiązywania zadania
| |
| nazywamy ''numerycznie poprawnym'' w zbiorze danych
| |
| <math>\displaystyle F_0\subsetF</math>, jeśli dla każdej danej <math>\displaystyle f\inF_0</math>
| |
| wynik <math>\displaystyle fl_\nu({\bf ALG}(f))</math> działania algorytmu w arytmetyce
| |
| <math>\displaystyle fl_\nu</math> można zinterpretować jako nieco zaburzony wynik
| |
| algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji
| |
| <math>\displaystyle y_\nu=(N(f))_\nu\inN(F)</math> o <math>\displaystyle f</math>, przy czym
| |
| poziom zaburzeń nie zależy od <math>\displaystyle f</math>.
| |
|
| |
| Formalnie znaczy to, że istnieją stałe <math>\displaystyle K_1</math>, <math>\displaystyle K_2</math>, oraz
| |
| <math>\displaystyle \nu_0>0</math> takie, że spełniony jest następujący warunek.
| |
| Dla dowolnej <math>\displaystyle \nu\le\nu_0</math> oraz informacji <math>\displaystyle y\inN(F_0)</math>
| |
| można dobrać <math>\displaystyle y_\nu\inN(F)</math> oraz
| |
| <math>\displaystyle \left(\varphi(y_\nu)\right)_\nu</math> takie, że
| |
|
| |
| <center><math>\displaystyle \|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|,
| |
| </math></center>
| |
|
| |
| <center><math>\displaystyle \|\left(\varphi(y_\nu)\right)_\nu - \varphi(y_\nu)\|\,\le\,
| |
| K_2\,\nu\,\|\varphi(y_\nu)\|,
| |
| </math></center>
| |
|
| |
| oraz
| |
|
| |
| <center><math>\displaystyle fl_\nu\left({\bf ALG}(f)\right)\,=\,
| |
| fl_\nu\left(\varphi(N(f))\right)\,=\,
| |
| \left(\varphi(y_\nu)\right)_\nu.
| |
| </math></center>
| |
|
| |
| <math>\displaystyle \Box</math> ||
| |
|
| |
| }}
| |
|
| |
| [[Image:MNcondition7.png|thumb|300px|Numerycznie poprawny algorytm daje w arytmetyce <math>\displaystyle fl_\nu</math> wynik <math>\displaystyle ALG(N(x))</math>, który daje
| |
| się zinterpretować jako mało zaburzony wynik <math>\displaystyle f(y)</math> zadania na mało zaburzonych
| |
| danych <math>\displaystyle x</math>.]]
| |
|
| |
| Zauważmy,że jeśli <math>\displaystyle f\inR^n</math>,
| |
| <math>\displaystyle N(f)=(f_1,\ldots,f_n)</math>, oraz algorytm jest
| |
| dokładny, <math>\displaystyle {\bf ALG}\equiv\varphi\equivS</math>, to numeryczną
| |
| poprawność algorytmu można równoważnie zapisać jako
| |
|
| |
| <center><math>\displaystyle fl_\nu\left({\bf ALG}(f)\right)\,=\,
| |
| \left(S(f_\nu)\right)_\nu.
| |
| </math></center>
| |
|
| |
| ===Rola uwarunkowania zadania===
| |
|
| |
| Niech <math>\displaystyle {\bf ALG}(\cdot)=\varphi(N(\cdot))</math> będzie algorytmem numerycznie
| |
| poprawnym dla danych <math>\displaystyle F_0\subsetF</math>. Wtedy jego błąd w <math>\displaystyle fl_\nu</math>
| |
| można oszacować następująco:
| |
|
| |
| <center><math>\displaystyle \aligned \lefteqn{ \|S(f)-fl_\nu({\bf ALG}(f))\| \;=\;
| |
| \|S(f)-\left(\varphi(y_\nu)\right)_\nu\| } \\
| |
| &\le & \|S(f)-\varphi(y)\|\,+\,
| |
| \|\varphi(y)-\varphi(y_\nu)\|\,+\,
| |
| \|\varphi(y_\nu)-\left(\varphi(y_\nu)\right)_\nu\| \\
| |
| &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
| |
| \|\varphi(y)-\varphi(y_\nu)\|\,+\,
| |
| K_2\,\nu\,\|\varphi(y_\nu)\| \\
| |
| &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
| |
| (1 + K_2 \nu) \|\varphi(y)-\varphi(y_\nu)\|\,+\,
| |
| K_2\,\nu\,\|\varphi(y)\|,
| |
| \endaligned</math></center>
| |
|
| |
| przy czym <math>\displaystyle \|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|</math>. Stąd
| |
| w szczególności wynika, że jeśli algorytm jest numerycznie
| |
| poprawny i ciągły ze względu na informację <math>\displaystyle y</math>, to
| |
|
| |
| <center><math>\displaystyle \lim_{\nu\to 0}\,\|S(f)-fl_\nu({\bf ALG}(f))\|\,=\,
| |
| \|S(f)-{\bf ALG}(f)\|.
| |
| </math></center>
| |
|
| |
| To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie
| |
| się zachowywał w <math>\displaystyle fl_\nu</math> prawie tak jak w arytmetyce idealnej.
| |
|
| |
| Z powyższych wzorów wynika, że błąd w <math>\displaystyle fl_\nu</math> algorytmu
| |
| numerycznie poprawnego zależy w dużym stopniu od:
| |
| * dokładności algorytmu w arytmetyce idealnej,
| |
| * dokładności <math>\displaystyle \nu</math> arytmetyki <math>\displaystyle fl_\nu</math>,
| |
| * wrażliwości algorytmu na małe względne zaburzenia
| |
| informacji <math>\displaystyle y</math>.
| |
|
| |
| Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy
| |
| trochę więcej uwagi jedynie trzeciemu.
| |
|
| |
| Jeśli <math>\displaystyle \varphi</math> spełnia warunek Lipschitza ze stałą <math>\displaystyle L</math>,
| |
| a dokładniej
| |
|
| |
| <center><math>\displaystyle \|\varphi(y_\nu)-\varphi(y)\|\,\le\,L\,\|y_\nu-y\|,
| |
| </math></center>
| |
|
| |
| to
| |
|
| |
| <center><math>\displaystyle \aligned \lefteqn{ \|S(f)-fl_\nu({\bf ALG}(f))\|} \\
| |
| &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
| |
| (1+K_2\nu)L\|y_\nu-y\|\,+\,K_2\nu\|\varphi(y)\| \\
| |
| &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
| |
| (1+K_2\nu)LK_1\nu\|y\|\,+\,K_2\nu\|\varphi(y)\|.
| |
| \endaligned</math></center>
| |
|
| |
| W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny
| |
| algorytmu proporcjonalnie do <math>\displaystyle \nu</math>.
| |
|
| |
| Bardziej jednak interesuje nas błąd ''względny''. Wybierzmy
| |
| "małe" <math>\displaystyle \eta\ge 0</math> i przypuśćmy, że
| |
|
| |
| <center><math>\displaystyle \|\varphi(y_\nu)-\varphi(y)\|\;\le\;
| |
| M\,K_1\,\nu\,\max(\eta,\|\varphi(y)\|),
| |
| </math></center>
| |
|
| |
| dla pewnej <math>\displaystyle M</math> niezależnej od <math>\displaystyle y</math>, tzn. błąd względny informscji,
| |
| <math>\displaystyle \|y_\nu-y\|\le K_1\nu\|y\|</math>, przenosi się na błąd względny
| |
| wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia"
| |
| <math>\displaystyle M</math>, albo na błąd bezwzględny ze współczynnikiem <math>\displaystyle M\eta</math>.
| |
| (Zauważmy, że gdybyśmy wzięli <math>\displaystyle \eta=0</math> to dla <math>\displaystyle y</math> takiej, że
| |
| <math>\displaystyle \varphi(y)=0</math> musiałoby być <math>\displaystyle \varphi(y_\nu)=0</math>, co zwykle, choć
| |
| nie zawsze, nie jest prawdą.) Wtedy
| |
|
| |
| <center><math>\displaystyle \aligned \lefteqn{ \|S(f)-fl_\nu({\bf ALG}(f))\| } \\
| |
| & \le & \|S(f)-{\bf ALG}(f)\|+
| |
| (1 + K_2 \nu) M K_1 \nu \max (\eta, \|\varphi(y)\|)+
| |
| K_2 \nu \|\varphi(y)\| \\
| |
| &= \|S(f)-{\bf ALG}(f)\|\,+\,\nu\,
| |
| \Big(\,MK_1(1+K_2\nu)+K_2\Big)\max(\eta,\|\varphi(y)\|).
| |
| \endaligned</math></center>
| |
|
| |
| W szczególności, gdy algorytm jest dokładny i korzysta z pełnej
| |
| informacji o <math>\displaystyle f</math>, tzn. <math>\displaystyle S\equiv{\bf ALG}\equiv\varphi</math>, to
| |
| błąd
| |
|
| |
| <center><math>\displaystyle \frac{\|S(f)-fl_\nu({\bf ALG}(f))\|}
| |
| {\max (\eta, \|S(f)\|)} \;\le\;
| |
| \Big( M K_1 (1+K_2\nu) + K_2\Big)\,\nu
| |
| \,\approx\,(M\,K_1\,+\,K_2)\,\nu.
| |
| </math></center>
| |
|
| |
| Stąd wynika, że jeśli <math>\displaystyle (MK_1+K_2)\nu\ll 1</math> to błąd względny
| |
| algorytmu w <math>\displaystyle fl_\nu</math> jest mały, o ile <math>\displaystyle \|S(f)\|\ge\eta</math>.
| |
| Błąd względny jest proporcjonalny do dokładności <math>\displaystyle \nu</math>,
| |
| arytmetyki <math>\displaystyle fl_\nu</math>, współczynników proporcjonalności <math>\displaystyle K_i</math>
| |
| algorytmu numerycznie poprawnego, oraz do wrażliwości <math>\displaystyle M</math>
| |
| zadania <math>\displaystyle S</math> na małe względne zaburzenia danych.
| |
|
| |
| Zwróćmy uwagę na istotny fakt, że interesują nas właściwie
| |
| tylko te zaburzenia danych (informacji), które powstają przy
| |
| analizie algorytmu numerycznie poprawnego. I tak, jeśli algorytm
| |
| jest numerycznie poprawny z pozornymi zaburzeniami danych w normie,
| |
| to trzeba zbadać wrażliwość zadania ze względu na zaburzenia
| |
| danych w normie. Jeśli zaś mamy pozorne zaburzenia
| |
| "po współrzędnych" (co oczywiście implikuje pozorne zaburzenia
| |
| w normie) to wystarczy zbadać uwarunkowanie zadania ze względu na
| |
| zaburzenia "po współrzędnych", itd.
| |
|
| |
| Zadania, które nie są zbyt wrażliwe na "małe" względne
| |
| zaburzenia danych, tzn. dla których <math>\displaystyle M</math> jest "niewielkie",
| |
| nazywamy ogólnie zadaniami ''dobrze uwarunkowanymi''.
| |
|
| |
| {{przyklad|Iloczyn skalarny||
| |
|
| |
| Załóżmy. że dla danych ciągów rzeczywistych o ustalonej
| |
| długości <math>\displaystyle n</math>, <math>\displaystyle a_j</math>, <math>\displaystyle b_j</math>, <math>\displaystyle 1\le j\le n</math>, chcemy obliczyć
| |
|
| |
| <center><math>\displaystyle S(a,b)\,=\,\sum_{j=1}^n a_j b_j.
| |
| </math></center>
| |
|
| |
| Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem
| |
| i korzystający z pełnej informacji o kolejnych współrzednych.
| |
|
| |
| Oznaczmy przez <math>\displaystyle \tilde a_j</math> i <math>\displaystyle \tilde b_j</math> reprezentacje liczb
| |
| <math>\displaystyle a_j</math> i <math>\displaystyle b_j</math> w <math>\displaystyle fl_\nu</math>, <math>\displaystyle \tilde a_j=a_j(1+\alpha_j)</math>,
| |
| <math>\displaystyle \tilde b_j=b_j(1+\beta_j)</math>, oraz przez <math>\displaystyle \gamma_j</math> i <math>\displaystyle \delta_j</math>
| |
| błędy względne powstałe przy kolejnych mnożeniach i dodawaniach.
| |
| Oczywiście <math>\displaystyle |\alpha_j|,|\beta_j|, |\gamma_j|, |\delta_j|\le\nu</math>.
| |
| Otrzymujemy
| |
|
| |
| <center><math>\displaystyle \aligned \lefteqn{fl_\nu\left(\sum_{j=1}^n a_jb_j\right) \;=\;
| |
| \Big(\,(fl_\nu(\sum_{j=1}^{n-1}a_jb_j)\,+\,\tilde a_n\tilde b_n
| |
| (1+\gamma_n)\,\Big)(1+\delta_n)\,=\,\ldots } \\
| |
| &= \bigg(\cdots\Big(
| |
| \tilde a_1\tilde b_1(1+\gamma_1)+\tilde a_2\tilde b_2
| |
| (1+\gamma_2)\Big)(1+\delta_2) \\
| |
| & & \qquad\qquad\qquad\qquad +\cdots+
| |
| \tilde a_n\tilde b_n(1+\gamma_n)\bigg)(1+\delta_n) \\
| |
| &= \tilde a_1\tilde b_1(1+\gamma_1)(1+\delta_2)
| |
| \cdots(1+\delta_n)\\
| |
| & & \qquad\qquad\qquad\qquad +\cdots+\tilde a_j
| |
| \tilde b_j(1+\gamma_j)(1+\delta_j)\cdots(1+\delta_n) \\
| |
| &= \sum_{j=1}^n a_jb_j(1+e_j),
| |
| \endaligned</math></center>
| |
|
| |
| gdzie w przybliżeniu (tzn. gdy <math>\displaystyle \nu\to 0</math>) mamy <math>\displaystyle |e_1|\leq (n+2)\nu</math>
| |
| i <math>\displaystyle |e_j|\leq (n-j+4)\nu</math>, <math>\displaystyle 2\le j\le n</math>. Algorytm naturalny jest więc
| |
| numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany
| |
| w <math>\displaystyle fl_\nu</math> można zinterpretować jako dokładny wynik dla danych
| |
| <math>\displaystyle a_{\nu,j}=a_j</math> i <math>\displaystyle b_{\nu,j}=b_j(1+e_j)</math>, przy czym
| |
| <math>\displaystyle \|b_\nu-b\|_p\leq (n+2)\nu\|b\|_p</math>.
| |
|
| |
| Zobaczmy teraz, jak błąd we współrzędnych <math>\displaystyle b_j</math> wpływa
| |
| na błąd wyniku. Mamy
| |
|
| |
| <center><math>\displaystyle \aligned \Big|\sum_{j=1}^n a_jb_j-fl_\nu\Big(\sum_{j=1}^n a_jb_j\Big)\Big|
| |
| &= \Big|\sum_{j=1}^n a_jb_j-\sum_{j=1}^n a_jb_j(1+e_j)\Big| \\
| |
| &= \Big|\sum_{j=1}^n e_ja_jb_j\Big|
| |
| \,\le\, \sum_{j=1}^n |e_j||a_jb_j| \\
| |
| &\leq (n+2)\nu\sum_{j=1}^n |a_jb_j|.
| |
| \endaligned</math></center>
| |
|
| |
| Stąd dla <math>\displaystyle \eta\ge 0</math>
| |
|
| |
| <center><math>\displaystyle \frac{|\sum_{j=1}^n a_jb_j-fl_\nu(\sum_{j=1}^n a_jb_j)|}
| |
| {\max(\eta,|\sum_{j=1}^n a_jb_j|)} \,\leq\,
| |
| K_{\eta}\,(n+2)\,\nu,
| |
| </math></center>
| |
|
| |
| gdzie
| |
|
| |
| <center><math>\displaystyle K_\eta\,=\,K_\eta(a,b)\,=\,\frac{\sum_{j=1}^n |a_jb_j|}
| |
| {\max(\eta,|\sum_{j=1}^n a_jb_j|) }.
| |
| </math></center>
| |
|
| |
| Zauważmy, że jeśli iloczyny <math>\displaystyle a_jb_j</math> są wszystkie dodatnie
| |
| albo wszystkie ujemne, to <math>\displaystyle K_\eta=1</math>, tzn. zadanie jest dobrze
| |
| uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej
| |
| <math>\displaystyle n\nu</math>. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile
| |
| liczba <math>\displaystyle n</math> składników nie jest horendalnie duża. W ogólności
| |
| jednak <math>\displaystyle K_\eta</math> może być liczbą dowolnie dużą i wtedy nie możemy
| |
| być pewni uzyskania dobrego wyniku w <math>\displaystyle fl_\nu</math>.
| |
|
| |
| }}
| |
|
| |
| {{przyklad|Pierwiastki trójmianu||
| |
|
| |
| Rozpatrzymy teraz zadanie obliczenia wszystkich pierwiastków
| |
| rzeczywistych równania kwadratowego.
| |
| Będziemy zakładać, że model obliczeniowy dopuszcza obliczanie
| |
| pierwiastków kwadratowych z liczb nieujemnych oraz
| |
| <math>\displaystyle fl_\nu(\sqrt{x})=rd_\nu(\sqrt{rd_\nu(x)})</math>.
| |
|
| |
| Okazuje się, że nie umiemy pokazać numerycznej poprawności
| |
| "szkolnego" algorytmu obliczającego pierwiastki równania
| |
| bezpośrednio ze wzorów ([[##szkolny|Uzupelnic: szkolny ]]). Można jednak pokazać
| |
| numeryczną poprawność drobnej jego modyfikacji wykorzystującej
| |
| wzory Viete'a.
| |
|
| |
| {{algorytm|||
| |
| <pre>
| |
| Delta <nowiki>=</nowiki> p*p - q;
| |
| if (Delta <nowiki>=</nowiki><nowiki>=</nowiki> 0)
| |
| OUT(p);
| |
| else
| |
| if (Delta > 0)
| |
| {
| |
| Delta1 <nowiki>=</nowiki> sqrt(d);
| |
| if (p ><nowiki>=</nowiki> 0)
| |
| {
| |
| x1 <nowiki>=</nowiki> p + Delta1;
| |
| x2 <nowiki>=</nowiki> q/z1;
| |
| }
| |
| else
| |
| {
| |
| x2 <nowiki>=</nowiki> p - Delta1;
| |
| x1 <nowiki>=</nowiki> q/ź2;
| |
| }
| |
| OUT(x1); OUT(x2);
| |
| }
| |
| </pre>}}
| |
|
| |
| Mamy bowiem
| |
|
| |
| <center><math>\displaystyle \aligned fl_\nu(\Delta(p,q)) &= \Big(p^2(1+\alpha)^2(1+\epsilon_1)-q(1+\beta)\Big)
| |
| (1+\epsilon_2) \\
| |
| &= \left( p^2-q\frac{(1+\beta)}{(1+\alpha)^2(1+\epsilon_1)}\right)
| |
| (1+\epsilon_2)(1+\alpha)^2(1+\epsilon_1) \\
| |
| &= \Big(p^2-q(1+\delta)\,\Big)(1+\gamma) \,=\,
| |
| \Delta(p,q(1+\delta))(1+\gamma),
| |
| \endaligned</math></center>
| |
|
| |
| gdzie <math>\displaystyle |\delta|,|\gamma|\leq 4\nu</math>. Wyróżnik obliczony w <math>\displaystyle fl_\nu</math>
| |
| jest więc nieco zaburzonym wyróżnikiem dokładnym dla danych
| |
| <math>\displaystyle p</math> i <math>\displaystyle q_\nu=q(1+\delta)</math>. W szczególności
| |
|
| |
| <center><math>\displaystyle \mbox{sgn} (fl_\nu(\Delta(p,q)))= \mbox{sgn} (\Delta(p,q_\nu)).
| |
| </math></center>
| |
|
| |
| Jeśli <math>\displaystyle p\ge 0</math> to
| |
|
| |
| <center><math>\displaystyle \aligned fl_\nu(x1(p,q)) &= \Big(p(1+\alpha)+
| |
| \sqrt{fl_\nu(\Delta(p,q))}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
| |
| &= \Big(p(1+\alpha)+\sqrt{\Delta(p,q_\nu)
| |
| (1+\gamma)}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
| |
| &= \left(p+\sqrt{\Delta(p,q_\nu)}\frac{\sqrt{1+\gamma}(1+\epsilon_3)}
| |
| {1+\alpha}\right)(1+\epsilon_4)(1+\alpha) \\
| |
| &= \left(p+\sqrt{\Delta(p,q_\nu)}\right)(1+e_1),
| |
| \endaligned</math></center>
| |
|
| |
| gdzie <math>\displaystyle |e_1|\leq 6\nu</math>. Zauważmy, że ostatnia równość
| |
| zachodzi dlatego, że dodajemy liczby tego samego znaku. (Inaczej
| |
| <math>\displaystyle |e_1|</math> mogłaby być dowolnie duża i tak byłoby w algorytmie
| |
| szkolnym.) Dla drugiego pierwiastka mamy
| |
|
| |
| <center><math>\displaystyle fl_\nu(x2(p,q))\,=\,\frac {q(1+\beta)}{fl_\nu(x1(p,q))}(1+\epsilon_5)
| |
| \,=\,\frac{q_\nu}{fl_\nu(x1(p,q))}(1+e_2),
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle |e_2|\le 8\nu</math>.
| |
|
| |
| Podobny wynik otrzymalibyśmy dla <math>\displaystyle p<0</math>. Algorytm zmodyfikowany
| |
| jest więc numerycznie poprawny, gdyż otrzymane w <math>\displaystyle fl_\nu</math> pierwiastki
| |
| są nieco zaburzonymi dokładnymi pierwiatkami dla danych
| |
| <math>\displaystyle p_\nu=p</math> i <math>\displaystyle q_\nu=q(1+\delta)</math>.
| |
|
| |
| Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie
| |
| zadania ze względu na zaburzenie danej <math>\displaystyle q</math>, ponieważ pokazaliśmy,
| |
| że zaburzenia <math>\displaystyle p</math> można przenieść na zaburzenia <math>\displaystyle q</math> i wyniku.
| |
| Niestety, choć algorytm jest numerycznie poprawny, zaburzenia
| |
| <math>\displaystyle q</math> mogą sprawić, że nawet znak wyróżnika <math>\displaystyle \Delta</math> może być
| |
| obliczony nieprawidłowo. Na przykład dla <math>\displaystyle p=1</math> i <math>\displaystyle q=1\pm 10^{t+1}</math>
| |
| mamy <math>\displaystyle \Delta(p,q)=\mp 10^{t+1}</math>, ale
| |
| <math>\displaystyle \Delta(rd_\nu(p),rd_\nu(q))=\Delta(1,1)=0</math>. Ogólnie
| |
|
| |
| <center><math>\displaystyle |fl_\nu(\Delta(p,q))-\Delta(p,q)|\,\leq\,4\nu(p^2+2|q|),
| |
| </math></center>
| |
|
| |
| a więc tylko dla <math>\displaystyle |\Delta(p,q)|=|p^2-q|>4\nu (p^2+2|q|)</math>
| |
| możemy być pewni obliczenia właściwego znaku <math>\displaystyle \Delta</math>. Przy
| |
| tym warunku oraz <math>\displaystyle \Delta>0</math> błąd danych przenosi się w
| |
| normie euklidesowej na błąd wyniku następująco:
| |
|
| |
| <center><math>\displaystyle \aligned \lefteqn{ \Big( (x1(p,q) - x1(p,q_\nu))^2
| |
| +(x2(p,q) - x2(p,q_\nu))^2 \Big)^{1//2} } \\
| |
| &= \frac{\sqrt 2 |\delta q|} {\sqrt{p^2-q}+\sqrt{p^2-q_\nu}}
| |
| \,\leq\, 2\sqrt 2 \nu\frac{|q|}{\sqrt{p^2-q}} \\
| |
| &= 2\sqrt 2 \nu\frac{|q|//p^2}{\sqrt{1-q//p^2}
| |
| \max(\eta//|p|,\sqrt{2(1+(1-q//p^2))}) } \\
| |
| & & \qquad\qquad\qquad\cdot\max(\eta,(x1(p,q)^2+x2(p,q)^2)^{1//2}).
| |
| \endaligned</math></center>
| |
|
| |
| Stąd widać, że zadanie jest dobrze uwarunkowane dla <math>\displaystyle q//p^2\ll 1</math>
| |
| i może być źle uwarunkowane dla <math>\displaystyle q//p^2\approx 1</math>. W ostatnim
| |
| przypadku nie możemy być pewni otrzymania dobrego wyniku w <math>\displaystyle fl_\nu</math>.
| |
| }}
| |
