MN04: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
mNie podano opisu zmian
m Zastępowanie tekstu – „,↵</math>” na „</math>,”
 
(Nie pokazano 28 wersji utworzonych przez 3 użytkowników)
Linia 1: Linia 1:
''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki''


==Ćwiczenia: równania nieliniowe skalarne==
<!--
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki
wprowadzi� przykry@mimuw.edu.pl
-->
=Własności zadania obliczeniowego i algorytmu numerycznego=


{{cwiczenie|Metoda Newtona może być zbieżna globalnie||
{{powrot |Metody numeryczne | do strony głównej
przedmiotu <strong>Metody numeryczne</strong>}}


Wykaż, że jeśli <math>\displaystyle f</math> jest rosnąca i wypukła na <math>\displaystyle [a,b]</math> oraz dla <math>\displaystyle x^*\in [a,b]\displaystyle f(x^*) = 0</math>, to metoda Newtona startująca z <math>\displaystyle x_0>x^*</math> jest zbieżna.
==Uwarunkowanie zadania obliczeniowego==
}}


<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 zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując
Powtarzając rachunki z dowodu lokalnej zbieżności, możemy łatwo stwierdzić, że kolejne
zadanie obliczeniowe, są z natury rzeczy wartościami zaburzonymi. Źródłem tych
iteracje dają <math>\displaystyle x_k</math> które także spełniają warunek <math>\displaystyle x_k - x^* > 0</math>. Ponadto
zaburzeń są:
nietrudno sprawdzić, że <math>\displaystyle x_{k+1} < x_k</math>, a więc mamy ciąg malejący i
* błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (np. 0.1 nie jest równe dokładnie <math>1/10</math>)
ograniczony. Jego granicą musi być <math>\displaystyle x^*</math> (dlaczego?).
* błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie <math>f(x) = a</math>, ale <math>a</math> jest rezultatem innej symulacji), a także
</div></div>
* błędy pomiarowe w zadaniach praktycznych (np. 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 <strong>wpływu
zaburzenia danych na wynik</strong> jest kapitalna w świecie obliczeń numerycznych w
ogólności, a w szególności --- inżynierskich.


{{cwiczenie|Fraktale||
Wprowadza się pojęcie <strong>uwarunkowania</strong> zadania, to znaczy jego podatności na
zaburzenia danych. Dla przejrzystości przypuśćmy, że nasze zadanie obliczeniowe
polega na wyznaczeniu <math>f(x)</math> dla danego <math>x</math>.


Ciekawy zbiór o charakterze fraktalnym powstaje w wyniku zastosowania metody
Newtona do rozwiązania równania <math>\displaystyle z^n-1=0</math> w dziedzinie zespolonej. Punkt <math>\displaystyle z_0 =
(x_0,y_0)</math> należy do ''basenu zbiezności metody'', jeśli metoda Newtona jest
zeń zbieżna do (jakiegokolwiek) pierwiastka w/w równania. Kolor odpowiadający
<math>\displaystyle z_0</math> jest określany na podstawie liczby iteracji potrzebnych metodzie do
zbieżności.


Zupełnie miłym (i estetycznie wartościowym) doświadczeniem, jest napisanie
Jak bardzo będzie odległe
programu w Octave, który wyświetla baseny zbieżności metody Newtona jak poniżej.
<math>f(\widetilde{x})</math>, gdy <math>\widetilde{x}\approx x</math>? Rozważa się dwa przypadki:
* <strong>uwarunkowanie względne</strong>: jak względne zaburzenie danych wpływa na błąd względny wyniku: <center><math>\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>\mbox{cond} _{rel}(f,x)</math> spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>.
* <strong>uwarunkowanie bezwzględne</strong>: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku: <center><math>||f(x) - f(\widetilde{x})|| \leq  \mbox{cond} _{abs}(f,x) \cdot ||x - \widetilde{x}||</math></center>
Najmniejszy mnożnik <math>\mbox{cond} _{abs}(f,x)</math>  spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>.
Powiemy, że zadanie <math>f(x)</math> jest
* <strong>dobrze uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \approx 1</math>,
* <strong>źle uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \gg 1</math>,
* <strong>źle postawione</strong> w punkcie <math>x</math>, gdy <math>\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!


[[Image:MNfraktal.png|400px|Basen zbiezności metody Newtona w okolicy początku układu
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
współrzędnych, dla równania <math>\displaystyle z^5-1=0</math>]]
<span  style="font-variant:small-caps;">Przykład: Uwarunkowanie zadania obliczenia sumy</span>
<div class="solution" style="margin-left,margin-right:3em;">


<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">
Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia <math>s(x,y) = x + y</math> ma
Bardzo wygodnie napisać taki program w postaci skryptu Octave. Pamiętaj
<center><math>
jednak, by uniknąć w nim pętli w rodzaju
\mbox{cond} _{abs}(s, (a,b)) = 1, \qquad  \mbox{cond} _{rel}(s, (a,b)) = \frac{|a|+|b|}{|a+b|}
</math></center>


for xpixels <nowiki>=</nowiki> 1:1024<br>
Tak więc, gdy <math>a\approx -b</math>, to <math>\mbox{cond} _{rel}(s, (a,b)) \approx +\infty</math> i zadanie
for ypixels <nowiki>=</nowiki> 1:768<br>
jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może
... wyznacz kolor pixela ....<br>
skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego,
end<br>
najczęściej rzeczywiście tak będzie...
end<br>
</div></div>
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład</span>  
<div class="solution" style="margin-left,margin-right:3em;">
 
Dla różniczkowalnej funkcji skalarnej <math>f : R \rightarrow R</math> mamy
<center><math>
|f(x) - f(\widetilde{x})| \approx |f'(x) | | x - \widetilde{x} |
</math></center>


Taka podwójnie zagnieżdżona pętla może skutecznie zniechęcić Cię do
i w konsekwencji dla zadania obliczenia <math>f(x)</math> dla danego <math>x</math> mamy, przy
Octave'a --- obrazek będzie liczyć się wieki --- zupełnie niepotrzebnie!
założeniu małych zaburzeń,
<center><math>
\mbox{cond} _{abs}( f, x) = |f'(x)|, \qquad  \mbox{cond} _{rel}( f, x) =
\frac{|f'(x)|\cdot|x|}{|f(x)|}</math></center>


</div></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">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
Jest zrozumiałe, że złe uwarunkowanie jest szkodliwe w praktyce numerycznej: zaburzenia danych są nieuniknione, a źle uwarunkowane zadanie tylko je wzmocni na wyjściu. Jednak za jakiś czas zobaczysz wyjątkową sytuację, gdy  [[MN13#Odwrotna metoda potęgowa|złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, ale wręcz <strong>pomaga</strong>]] szybciej rozwiązać zadanie główne!
<div style="font-size:smaller; background-color:#cfd">Kluczem do efektywnej implementacji w Octave jest przeprowadzenie
 
iteracji Newtona na ''wszystkich pikselach jednocześnie''. W tym celu musisz
==Rozkład algorytmu względem informacji==
skorzystać z prowadzenia działań na całych macierzach.</div>
 
</div></div>
<strong>Algorytm</strong> 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


{{cwiczenie|Pierwiastkowanie||
<center><math>{\bf ALG}:\,F\longrightarrow G,
</math></center>


Niech <math>\displaystyle a\geq 0</math>. Aby wyznaczyć <math>\displaystyle \sqrt{a}</math>, można skorzystać z metody Newtona dla
taki że <math>{\bf ALG}(f)</math> jest wynikiem działania algorytmu
równania <math>\displaystyle x^2 = a</math>. Zaprogramuj tę metodę i sprawdź, jak wiele dokładnych cyfr
w arytmetyce idealnej dla danej <math>f</math>.  
wyniku dostajesz w kolejnych krokach. Czy to możliwe, by liczba dokładnych cyfr
znaczących z grubsza podwajała się na każdej iteracji? Wskaż takie <math>\displaystyle a</math> dla
którego to nie będzie prawdą.


<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">
Zauważmy, że wynik <math>{\bf ALG}(f)</math> działania algorytmu nie
Na pewno kusi Cię, by zaprogramować najpierw ogólnie funkcję "metoda
zależy bezpośrednio od <math>f</math>, ale raczej od <strong>informacji</strong>  
Newtona", a następnie przekazać jej jako parametr naszą funkcję. To oczywiście
o <math>f</math> (uzyskanej dzięki poleceniu <math>{\cal IN}</math>). Informacja
będzie działać, i to praktycznie tak samo szybko jak drugi spsoób, w którym po prostu
ta może być <strong>pełna</strong> albo tylko <strong>częściowa</strong>.  
wyliczysz na karteczce wzór na iterację dla ''tej konkretnej funkcji, <math>\displaystyle F(x) =
Informacja jest pełna gdy, np.
x^2-a</math>'' i potem go zaimplementujesz. Jednak
<math>f=(f_1,\ldots,f_n)\in R^n</math> i wczytamy wszystkie
drugie podejście jest tutaj godne polecenia, bo ma na sobie znak charakterystyczny dla numeryki: niech
współrzędne <math>f_i</math>. Informacja może być częściowa, gdy
wszystko, co programujesz będzie tak proste, jak tylko możliwe (ale nie
<math>f</math> jest funkcją. Wtedy wiele danych może posiadać tę
prostsze)!
samą informację, co łatwo zaobserwować na przykładzie
</div></div>
zadania całkowania.


}}
Niech <math>N:F\to \cup_{n=0}^\infty R^n</math> będzie
<strong>operatorem informacji</strong>, tzn.


<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"> 
<center><math>N(f)\,=\,(y_1,y_2,\ldots,y_n)
Nasza iteracja, którą chcemy zaprogramować, to
<center><math>\displaystyle
x_{k+1} = \frac{1}{2}\left( x_k + \frac{a}{x_k}\right).
</math></center>
</math></center>


Dla <math>\displaystyle a\neq 0</math>, spełnione są założenia twierdzenia o zbieżności metody Newtona,
jest informacją o <math>f</math> zebraną przy idealnej realizacji
więc
algorytmu. Zauważmy, że informacja jest pełna, gdy <math>N</math> jest
przekształceniem różnowartościowym, tzn. jeśli
<math>f_1\ne f_2</math> implikuje <math>N(f_1)\ne N(f_2)</math>.
W przeciwnym przypadku mamy do czynienia z informacją
częściową.
 
Każdy algorytm <math>{\bf ALG}</math> może być przedstawiony jako złożenie
operatora informacji i pewnego operatora
<math>\varphi:N(F)\to G</math> zdefiniowanego równością


<center><math>\displaystyle |x_{k+1}-x^*|\leq K |x_{k}-x^*|^2
<center><math>\varphi\left(N(f)\right)\,=\,{\bf ALG}(f).
</math></center>
</math></center>


co właśnie tłumaczy się jako podwajanie liczby dokładnych cyfr
Zauważmy, że w przypadku informacji częściowej zwykle nie
znaczących wyniku.
istnieje algorytm dający dokładne rozwiązanie zadania dla
każdej danej <math>f\in F</math>, ponieważ dla danych o tej samej
informacji mogą istnieć różne rozwiązania.
 
==Problem wyboru algorytmu==


Ale jeśli <math>\displaystyle a=0</math> to metoda zbieżna jest już tylko liniowo (z ilorazem 0.5 --- bo
Wybór algorytmu jest najistotniejszą częścią całego procesu
jaki wtedy jest wzór na iterację?), więc co mniej więcej trzy kroki będziemy
numerycznego rozwiązywania zadania. Kierujemy się przy tym przede
dostawali kolejne zero dokładnego wyniku <math>\displaystyle x^*=0</math>.
wszystkim następującymi kryteriami:
* dokładnością algorytmu,
* złożonością algorytmu,
* własnościami numerycznymi algorytmu.
   
   
</div></div>
Przez dokładność algorytmu rozumiemy różnicę między
rozwiązaniem dokładnym <math>S(f)</math> a rozwiązaniem
<math>{\bf ALG}(f)</math> dawanym przez algorytm w arytmetyce idealnej.
Jeśli <math>{\bf ALG}(f) = S(f)</math>,
<math>\forall f \in F</math>,
to algorytm nazywamy <strong>dokładnym</strong>.
 
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>f</math> składa
się koszt uzyskania infomacji <math>y=N(f)</math> (zwykle jest on
proporcjonalny do liczby wywołań polecenia <math>{\cal IN}</math>), oraz
koszt <strong>kombinatoryczny</strong> przetworzenia tej informacji, aż do
uzyskania wyniku <math>\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>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>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>fl_\nu</math> możemy
otrzymać wynik <math>fl_\nu({\bf ALG}(f))</math> daleko odbiegający od
<math>S(f)</math>. W szczególności, prawie zawsze mamy
 
<center><math>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>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>y=N(f)</math> o danej <math>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 <strong>nieco zaburzonej</strong> <math>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>fl_\nu</math>
będzie <math>(\varphi(y_\nu))_\nu</math> zamiast <math>\varphi(y)</math>. Algorytmy
dające tego rodzaju wyniki uznamy za posiadające najlepsze
własności numeryczne w arytmetyce <math>fl_\nu</math> i nazwiemy numerycznie
poprawnymi.
 
Powiemy, że ciąg rzeczywisty
<math>a_\nu=(a_{\nu,1},\ldots,a_{\nu,n})</math>
(a właściwie rodzina ciągów <math>\{a_\nu\}_\nu</math>) jest
<strong>nieco zaburzonym</strong> ciągiem <math>a=(a_1,\ldots,a_n)</math>, jeśli
istnieje stała <math>K</math> taka, że dla wszystkich dostatecznie
małych <math>\nu</math> zachodzi
 
<center><math>
  |a_{\nu,j} - a_j|\,\le\,K\,\nu\,|a_j|,\qquad 1\le j\le n,
</math></center>
 
albo ogólniej
 
<center><math>
  \|a_\nu - a\|\,\le\,K\,\nu\,\|a\|</math>,</center>
 
gdzie <math>\|\cdot\|</math> jest pewną normą w <math>R^n</math>. W pierwszym
przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim
o zaburzeniu w normie <math>\|\cdot\|</math>.
 
Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają
za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas
 
<center><math>\|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>K_1</math> i <math>K_2</math>
 
<center><math>\|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ść dla zaburzenia w normie, ze stałą <math>K = K_2 K_1 K</math>.
 
{{definicja|Algorytm numerycznie poprawny|Algorytm numerycznie poprawny|
 
Algorytm <math>{\bf ALG}</math> rozwiązywania zadania
nazywamy <strong>numerycznie poprawnym</strong> w zbiorze danych
<math>F_0\subset F</math>, jeśli dla każdej danej <math>f\in F_0</math>
wynik <math>fl_\nu({\bf ALG}(f))</math> działania algorytmu w arytmetyce
<math>fl_\nu</math> można zinterpretować jako nieco zaburzony wynik
algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji
<math>y_\nu=(N(f))_\nu\in N(F)</math> o <math>f</math>, przy czym
poziom zaburzeń nie zależy od <math>f</math>.
 
Formalnie znaczy to, że istnieją stałe <math>K_1</math>, <math>K_2</math>  oraz
<math>\nu_0>0</math> takie, że spełniony jest następujący warunek.
Dla dowolnej <math>\nu\le\nu_0</math> oraz informacji <math>y\in N(F_0)</math>
można dobrać <math>y_\nu\in N(F)</math> oraz
<math>\left(\varphi(y_\nu)\right)_\nu</math> takie, że
 
<center><math>\|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|,
</math></center>
 
<center><math>\|\left(\varphi(y_\nu)\right)_\nu - \varphi(y_\nu)\|\,\le\,
    K_2\,\nu\,\|\varphi(y_\nu)\|</math>,</center>


{{cwiczenie|Odwrotność bez dzielenia||
oraz


Aby wyznaczyć <math>\displaystyle 1/a</math> dla <math>\displaystyle a>0</math> bez dzielenia(!), można zastosować metodę Newtona
<center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\,
do funkcji <math>\displaystyle f(x) = 1/x - a</math>. Pokaż, że na <math>\displaystyle k</math>-tym kroku iteracji,
      fl_\nu\left(\varphi(N(f))\right)\,=\,
<center><math>\displaystyle
      \left(\varphi(y_\nu)\right)_\nu.  
|ax_k - 1| = |ax_{k-1} - 1|^2.
</math></center>
</math></center>


Dla jakich <math>\displaystyle x_0</math> metoda będzie zbieżna do <math>\displaystyle \frac{1}{a}</math>, a dla jakich nie?
Oceń, ile iteracji potrzeba do spełnienia warunku <math>\displaystyle \frac{|x_k-\frac{1}{a}|}{|a|} \leq
\epsilon</math>, gdy <math>\displaystyle \frac{|x_0-\frac{1}{a}|}{|a|} \leq \gamma < 1</math>,
}}
}}


<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">
[[Image:MNcondition7.png|thumb|550px|center|Numerycznie poprawny algorytm daje w arytmetyce <math>fl_\nu</math> wynik <math>ALG(N(x))</math>, który daje
Proste ćwiczenie z analizy. Warunkiem zbieżności jest oczywiście (ze względu na
się zinterpretować jako mało zaburzony wynik <math>f(y)</math> zadania na mało zaburzonych
równość w wyrażeniu na błąd iteracji powyżej) <math>\displaystyle |ax_0 - 1| < 1</math>, to znaczy
danych <math>x</math>.]]
<center><math>\displaystyle 0 < x_0 < \frac{2}{a}.</math></center>
 
Zauważmy,że jeśli <math>f\in R^n</math>,
<math>N(f)=(f_1,\ldots,f_n)</math>, oraz algorytm jest  
dokładny, <math>{\bf ALG}\equiv\varphi\equiv S</math>, to numeryczną
poprawność algorytmu można równoważnie zapisać jako
 
<center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\,
  \left(S(f_\nu)\right)_\nu.  
</math></center>


</div></div>
Numeryczna poprawność algorytmu jest wielce pożądaną cechą.
 
<blockquote  style="background-color: #fefeee; padding:1em;  margin-left,margin-right:2em;  margin-top,margin-bottom: 1em;"> 
Algorytm numerycznie poprawny działa w arytmetyce zmiennoprzecinkowej (niemal) tak, jakby wszystkie obliczenia były wykonywane w arytmetyce dokładnej, a tylko dane i wyniki podlegały reprezentacji w skończonej precyzji.
</blockquote>
 
==Rola uwarunkowania zadania==
 
Niech <math>{\bf ALG}(\cdot)=\varphi(N(\cdot))</math> będzie algorytmem numerycznie
poprawnym dla danych <math>F_0\subset F</math>. Wtedy jego błąd w <math>fl_\nu</math>
można oszacować następująco:
 
<center><math>\begin{align} \|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)\|,
\end{align}</math></center>
 
przy czym <math>\|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>y</math>, to
 
<center><math>\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>fl_\nu</math> prawie tak, jak w arytmetyce idealnej.


{{cwiczenie|||
Z powyższych wzorów wynika, że błąd w <math>fl_\nu</math> algorytmu
Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji
numerycznie poprawnego zależy w dużym stopniu od:
* <math>\displaystyle f(x) = \sin(x)</math>,
* dokładności algorytmu w arytmetyce idealnej,  
* <math>\displaystyle f(x) = \sin^2(x)</math>,
* dokładności <math>\nu</math> arytmetyki <math>fl_\nu</math>,  
* <math>\displaystyle f(x) = (x-1)^5</math> (wzór A),
* wrażliwości algorytmu na małe względne zaburzenia informacji <math>y</math>.  
* <math>\displaystyle f(x) = (x-1)*(x^4 - 4x^3 + 6x^2 - 4x + 1)</math> (wzór B),
* <math>\displaystyle f(x) = (x-2)^13</math> (wzór C),
* <math>\displaystyle f(x) = x^13 - 26*x^12 + 312*x^11 - 2288*x^10 + ... - 8192</math> (wzór D),
   
   
gdy tolerancję błędu zadasz na poziomie <math>\displaystyle 10^{-10}</math>.
Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy
trochę więcej uwagi jedynie trzeciemu.
 
Jeśli <math>\varphi</math> spełnia warunek Lipschitza ze stałą <math>L</math>,
a dokładniej
 
<center><math>\|\varphi(y_\nu)-\varphi(y)\|\,\le\,L\,\|y_\nu-y\|</math>,</center>
 
to


Jak wyjaśnić te wyniki? Czy możesz już być pewna, że masz dobrą implementację?
<center><math>\begin{align} \|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)\|.
\end{align}</math></center>


}}
W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny
algorytmu proporcjonalnie do <math>\nu</math>.
 
Bardziej jednak interesuje nas błąd <strong>względny</strong>. Wybierzmy
"małe" <math>\eta\ge 0</math> i przypuśćmy, że
 
<center><math>\|\varphi(y_\nu)-\varphi(y)\|\;\le\;
    M\,K_1\,\nu\,\max(\eta,\|\varphi(y)\|)</math>,</center>
 
dla pewnej <math>M</math> niezależnej od <math>y</math>, tzn. błąd względny informacji,
<math>\|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>M</math>, albo na błąd bezwzględny ze współczynnikiem <math>M\eta</math>.
(Zauważmy, że gdybyśmy wzięli <math>\eta=0</math>, to dla <math>y</math> takiej, że
<math>\varphi(y)=0</math>, musiałoby być <math>\varphi(y_\nu)=0</math> --- co zwykle, choć
nie zawsze, nie jest prawdą.) Wtedy
 
<center><math>\begin{align} \|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)\|).
\end{align}</math></center>
 
W szczególności, gdy algorytm jest dokładny i korzysta z pełnej
informacji o <math>f</math>, tzn. <math>S\equiv{\bf ALG}\equiv\varphi</math>, to
błąd
 
<center><math>\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>(MK_1+K_2)\nu\ll 1</math>, to błąd względny
algorytmu w <math>fl_\nu</math> jest mały, o ile <math>\|S(f)\|\ge\eta</math>.
Błąd względny jest proporcjonalny do dokładności <math>\nu</math>,
arytmetyki <math>fl_\nu</math>, współczynników proporcjonalności <math>K_i</math>
algorytmu numerycznie poprawnego, oraz do wrażliwości <math>M</math>
zadania <math>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.
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Iloczyn skalarny</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
Załóżmy, że dla danych ciągów rzeczywistych o ustalonej
długości <math>n</math>, <math>a_j</math>, <math>b_j</math>, <math>1\le j\le n</math>, chcemy obliczyć
 
<center><math>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>\tilde a_j</math> i <math>\tilde b_j</math> reprezentacje liczb
<math>a_j</math> i <math>b_j</math> w <math>fl_\nu</math>, <math>\tilde a_j=a_j(1+\alpha_j)</math>,
<math>\tilde b_j=b_j(1+\beta_j)</math>, oraz przez <math>\gamma_j</math> i <math>\delta_j</math>
błędy względne powstałe przy kolejnych mnożeniach i dodawaniach.
Oczywiście <math>|\alpha_j|,|\beta_j|, |\gamma_j|, |\delta_j|\le\nu</math>.
Otrzymujemy
 
<center><math>\begin{align} 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) +\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) +\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),
\end{align}</math></center>
 
gdzie w przybliżeniu (tzn. gdy <math>\nu\to 0</math>) mamy <math>|e_1|\leq (n+2)\nu</math>
i <math>|e_j|\leq (n-j+4)\nu</math>, <math>2\le j\le n</math>. Algorytm naturalny jest więc
numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany
w <math>fl_\nu</math> można zinterpretować jako dokładny wynik dla danych
<math>a_{\nu,j}=a_j</math> i <math>b_{\nu,j}=b_j(1+e_j)</math>, przy czym
<math>\|b_\nu-b\|_p\leq (n+2)\nu\|b\|_p</math>.


<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">
Zobaczmy teraz, jak błąd we współrzędnych <math>b_j</math> wpływa
na błąd wyniku. Mamy


Jeśli nie pomylisz się, metoda powinna zbiegać bez problemów do zera funkcji
<center><math>\begin{align} \Big|\sum_{j=1}^n a_jb_j-fl_\nu\Big(\sum_{j=1}^n a_jb_j\Big)\Big|
<math>\displaystyle \sin</math>, dawać komunikat o błędzie dla <math>\displaystyle \sin^2</math> (bo nie ma przeciwnych znaków).
    &= \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|.
\end{align}</math></center>


Zapewne zauważyłaś, że wzory A i B są matematycznie równoważne, podobnie zresztą
Stąd dla <math>\eta\ge 0</math>
jak wzory C i D. Niestety, tylko wzór A i C dają w efekcie dobre przybliżenia
miejsca zerowego, podczas gdy wzory B i D prowadzą do oszacowania na miejsce
zerowe, które ''w ogóle nie zawierają'' prawdziwego miejsca zerowego.


Wyjaśnieniem tego fenomenu jest zjawisko redukcji cyfr przy odejmowaniu.
<center><math>\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>


[[Image:MNbisekcjawblad.png|400px|Metoda bisekcji ma kłopoty, gdy funkcja zadana jest
gdzie
wzorem D.]]
[[Image:MNbisekcjawresid.png|400px|Residuum też jest duże, gdy <math>\displaystyle f</math> zadana jest
wzorem D.]]


Jeśli chodzi o pewność... No cóż, sprawdziłaś, że działa w przypadkach, gdy
<center><math>K_\eta\,=\,K_\eta(a,b)\,=\,\frac{\sum_{j=1}^n |a_jb_j|}
spodziewałaś się, że będzie działać. Że tam, gdzie spodziewałaś się kłopotów,
            {\max(\eta,|\sum_{j=1}^n a_jb_j|) }.  
lub komunikatu o błędzie --- tak rzeczywiście się stało. Wreszcie,
</math></center>
potwierdziłaś, że zachowanie metody jest zgodne z jej znanymi właściwościami.


Tak więc, można ''przypuszczać'', że implementacja jest poprawna.
Zauważmy, że jeśli iloczyny <math>a_jb_j</math> są wszystkie dodatnie
albo wszystkie ujemne, to <math>K_\eta=1</math>, tzn. zadanie jest dobrze
uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej
<math>n\nu</math>. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile
liczba <math>n</math> składników nie jest horrendalnie duża. W ogólności
jednak <math>K_\eta</math> może być liczbą dowolnie dużą i wtedy nie możemy
być pewni uzyskania dobrego wyniku w <math>fl_\nu</math>.  


</div></div>
</div></div>


{{cwiczenie|||
<!--
Wskaż ''wszystkie'' wartości <math>\displaystyle x_0</math>, dla jakich metoda Newtona będzie zbieżna
 
do rozwiązania <math>\displaystyle x^*</math> równania
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Pierwiastki trójmianu</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
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>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 omawianych powyżej. Można jednak pokazać
numeryczną poprawność drobnej jego modyfikacji wykorzystującej
wzory Viete'a.
 
{{algorytm|||
<pre>\EATWSDelta = p*p - q;
if  (Delta == 0) 
      OUT(p);
else
if  (Delta > 0)
{
Delta1 = sqrt(d);
if  (p >= 0)
{
x1 = p + Delta1;
x2 = q/z1;
}
else
{
x2 = p - Delta1;
x1 = q/ź2;
}
OUT(x1);  OUT(x2);
}
</pre>}}
 
Mamy bowiem
 
<center><math>\begin{align} 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),
\end{align}</math></center>
 
gdzie <math>|\delta|,|\gamma|\leq 4\nu</math>. Wyróżnik obliczony w <math>fl_\nu</math>
jest więc nieco zaburzonym wyróżnikiem dokładnym dla danych
<math>p</math> i <math>q_\nu=q(1+\delta)</math>. W szczególności
 
<center><math>\mbox{sgn} (fl_\nu(\Delta(p,q)))= \mbox{sgn} (\Delta(p,q_\nu))</math></center>
 
Jeśli <math>p\ge 0</math> to
 
<center><math>\begin{align} 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),
\end{align}</math></center>
 
gdzie <math>|e_1|\leq 6\nu</math>. Zauważmy, że ostatnia równość
zachodzi dlatego, że dodajemy liczby tego samego znaku. (Inaczej
<math>|e_1|</math> mogłaby być dowolnie duża i tak byłoby w algorytmie
szkolnym.) Dla drugiego pierwiastka mamy


<center><math>\displaystyle \arctan(x) = 0.
<center><math>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>
</math></center>


Wyznacz wartość <math>\displaystyle X_0</math>, z którego startując powinieneś dostać ciąg oscylujący <math>\displaystyle X_0, -X_0,\ldots</math>.
gdzie <math>|e_2|\le 8\nu</math>.
Sprawdź eksperymentalnie, czy tak jest rzeczywiście.
 
Podobny wynik otrzymalibyśmy dla <math>p<0</math>. Algorytm zmodyfikowany
jest więc numerycznie poprawny, gdyż otrzymane w <math>fl_\nu</math> pierwiastki
są nieco zaburzonymi dokładnymi pierwiatkami dla danych
<math>p_\nu=p</math> i <math>q_\nu=q(1+\delta)</math>.
 
Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie
zadania ze względu na zaburzenie danej <math>q</math>, ponieważ pokazaliśmy,  
że zaburzenia <math>p</math> można przenieść na zaburzenia <math>q</math> i wyniku.
Niestety, choć algorytm jest numerycznie poprawny, zaburzenia
<math>q</math> mogą sprawić, że nawet znak wyróżnika <math>\Delta</math> może być
obliczony nieprawidłowo. Na przykład dla <math>p=1</math> i <math>q=1\pm 10^{t+1}</math>
mamy <math>\Delta(p,q)=\mp 10^{t+1}</math>, ale
<math>\Delta(rd_\nu(p),rd_\nu(q))=\Delta(1,1)=0</math>. Ogólnie


<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">
<center><math>|fl_\nu(\Delta(p,q))-\Delta(p,q)|\,\leq\,4\nu(p^2+2|q|)</math>,</center>
Aby znaleźć graniczne <math>\displaystyle X_0</math>, czyli takie, dla którego dostaniesz
oscylacje <math>\displaystyle X_0, -X_0,\ldots</math>, musisz rozwiązać równanie nieliniowe:


<center><math>\displaystyle -X_0 = X_0 - \frac{\arctan(X_0)}{\frac{1}{1+X_0^2}}.
a więc tylko dla <math>|\Delta(p,q)|=|p^2-q|>4\nu (p^2+2|q|)</math>
</math></center>
możemy być pewni obliczenia właściwego znaku <math>\Delta</math>. Przy
tym warunku oraz <math>\Delta>0</math> błąd danych przenosi się w
normie euklidesowej na błąd wyniku następująco:


To łatwe.
<center><math>\begin{align} \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}).
\end{align}</math></center>


Stąd widać, że zadanie jest dobrze uwarunkowane dla <math>q/p^2\ll 1</math>
i może być źle uwarunkowane dla <math>q/p^2\approx 1</math>. W ostatnim
przypadku nie możemy być pewni otrzymania dobrego wyniku w <math>fl_\nu</math>.
</div></div>
</div></div>


}}
-->
 
==Literatura==


<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"> 
W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 2.3</b> w
Dla <math>\displaystyle -X_0 < x_0 < X_0</math> mamy zbieżność. Dla <math>\displaystyle |x_0| = |X_0|</math> oscylacje, dla
* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.
większych wartości --- rozbieżność. Jednak w eksperymencie wszystko zależy od
tego, jak dokładnie wyznaczyliśmy <math>\displaystyle X_0</math>. Na przykład, mój solver w Octave
wyznaczył wartość <math>\displaystyle X_0</math>, dla której metoda Newtona dała następujący ciąg:
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
Iteracja:  1,   x <nowiki>=</nowiki> -1.39175
Iteracja:  2,  x <nowiki>=</nowiki> 1.39175
Iteracja:  3,  x <nowiki>=</nowiki> -1.39175
Iteracja:  4,  x <nowiki>=</nowiki> 1.39175
Iteracja:  5,  x <nowiki>=</nowiki> -1.39175


.. itd ..
Warto także przejrzeć rozdział 2 w
* <span style="font-variant:small-caps">P. Deulfhard, A. Hohmann</span>, <cite>Numerical Analysis in Modern Scientific Computing</cite>, Springer, 2003,


Iteracja:  25,  x <nowiki>=</nowiki> -1.39176
omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów.
Iteracja:  26,  x <nowiki>=</nowiki> 1.39178
Nieocenioną monografią na ten temat jest
Iteracja:  27,  x <nowiki>=</nowiki> -1.39183
* <span style="font-variant:small-caps">N. Higham</span>, <cite>Accuracy and Stability of Numerical Algorithms</cite>, SIAM, 2002.
Iteracja:  28,  x <nowiki>=</nowiki> 1.39197
Iteracja:  29,  x <nowiki>=</nowiki> -1.39235
Iteracja:  30,  x <nowiki>=</nowiki> 1.39333
Iteracja:  31,  x <nowiki>=</nowiki> -1.39594
Iteracja:  32,  x <nowiki>=</nowiki> 1.40283
Iteracja:  33,  x <nowiki>=</nowiki> -1.42117
Iteracja:  34,  x <nowiki>=</nowiki> 1.47061
Iteracja:  35,  x <nowiki>=</nowiki> -1.60867
Iteracja:  36,  x <nowiki>=</nowiki> 2.03161
Iteracja:  37,  x <nowiki>=</nowiki> -3.67722
Iteracja:  38,  x <nowiki>=</nowiki> 15.2779
Iteracja: 39,  x <nowiki>=</nowiki> -337.619
Iteracja:  40,  x <nowiki>=</nowiki> 178376
Iteracja:  41,  x <nowiki>=</nowiki> -4.99792e+10
Iteracja:  42,  x <nowiki>=</nowiki> 3.92372e+21
Iteracja:  43,   x <nowiki>=</nowiki> -2.41834e+43
Iteracja:  44,   x <nowiki>=</nowiki> 9.18658e+86
Iteracja:  45,   x <nowiki>=</nowiki> -1.32565e+174
Iteracja:  46,  x <nowiki>=</nowiki> inf
Iteracja:  47,  x <nowiki>=</nowiki> nan
</pre></div>
</div></div>

Aktualna wersja na dzień 21:46, 11 wrz 2023


Własności zadania obliczeniowego i algorytmu numerycznego

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

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 (np. 0.1 nie jest równe dokładnie 1/10)
  • błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie f(x)=a, ale a jest rezultatem innej symulacji), a także
  • błędy pomiarowe w zadaniach praktycznych (np. 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 f(x) dla danego x.


Jak bardzo będzie odległe f(x~), gdy x~x? Rozważa się dwa przypadki:

  • uwarunkowanie względne: jak względne zaburzenie danych wpływa na błąd względny wyniku:
    ||f(x)f(x~)||||f(x)||condrel(f,x)||xx~||||x||

Najmniejszy mnożnik condrel(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia f(x) dla danego x.

  • uwarunkowanie bezwzględne: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku:
    ||f(x)f(x~)||condabs(f,x)||xx~||

Najmniejszy mnożnik condabs(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia f(x) dla danego x.

Powiemy, że zadanie f(x) jest

  • dobrze uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle postawione w punkcie x, gdy cond(f,x)=+.

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!

Przykład: Uwarunkowanie zadania obliczenia sumy

Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia s(x,y)=x+y ma

condabs(s,(a,b))=1,condrel(s,(a,b))=|a|+|b||a+b|

Tak więc, gdy ab, to condrel(s,(a,b))+ 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...

Przykład

Dla różniczkowalnej funkcji skalarnej f:RR mamy

|f(x)f(x~)||f(x)||xx~|

i w konsekwencji dla zadania obliczenia f(x) dla danego x mamy, przy założeniu małych zaburzeń,

condabs(f,x)=|f(x)|,condrel(f,x)=|f(x)||x||f(x)|

Jest zrozumiałe, że złe uwarunkowanie jest szkodliwe w praktyce numerycznej: zaburzenia danych są nieuniknione, a źle uwarunkowane zadanie tylko je wzmocni na wyjściu. Jednak za jakiś czas zobaczysz wyjątkową sytuację, gdy złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, 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

𝐀𝐋𝐆:FG,

taki że 𝐀𝐋𝐆(f) jest wynikiem działania algorytmu w arytmetyce idealnej dla danej f.

Zauważmy, że wynik 𝐀𝐋𝐆(f) działania algorytmu nie zależy bezpośrednio od f, ale raczej od informacji o f (uzyskanej dzięki poleceniu 𝒩). Informacja ta może być pełna albo tylko częściowa. Informacja jest pełna gdy, np. f=(f1,,fn)Rn i wczytamy wszystkie współrzędne fi. Informacja może być częściowa, gdy f jest funkcją. Wtedy wiele danych może posiadać tę samą informację, co łatwo zaobserwować na przykładzie zadania całkowania.

Niech N:Fn=0Rn będzie operatorem informacji, tzn.

N(f)=(y1,y2,,yn)

jest informacją o f zebraną przy idealnej realizacji algorytmu. Zauważmy, że informacja jest pełna, gdy N jest przekształceniem różnowartościowym, tzn. jeśli f1f2 implikuje N(f1)N(f2). W przeciwnym przypadku mamy do czynienia z informacją częściową.

Każdy algorytm 𝐀𝐋𝐆 może być przedstawiony jako złożenie operatora informacji i pewnego operatora φ:N(F)G zdefiniowanego równością

φ(N(f))=𝐀𝐋𝐆(f).

Zauważmy, że w przypadku informacji częściowej zwykle nie istnieje algorytm dający dokładne rozwiązanie zadania dla każdej danej fF, 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 S(f) a rozwiązaniem 𝐀𝐋𝐆(f) dawanym przez algorytm w arytmetyce idealnej. Jeśli 𝐀𝐋𝐆(f)=S(f), fF, 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 f składa się koszt uzyskania infomacji y=N(f) (zwykle jest on proporcjonalny do liczby wywołań polecenia 𝒩), oraz koszt kombinatoryczny przetworzenia tej informacji, aż do uzyskania wyniku φ(y). 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 flν. 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 flν. Niestety, jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm jest dokładny, to w wyniku jego realizacji w flν możemy otrzymać wynik flν(𝐀𝐋𝐆(f)) daleko odbiegający od S(f). W szczególności, prawie zawsze mamy

S(f)flν(𝐀𝐋𝐆(f)).

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 flν. 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 y=N(f) o danej f 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 yν, 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 flν będzie (φ(yν))ν zamiast φ(y). Algorytmy dające tego rodzaju wyniki uznamy za posiadające najlepsze własności numeryczne w arytmetyce flν i nazwiemy numerycznie poprawnymi.

Powiemy, że ciąg rzeczywisty aν=(aν,1,,aν,n) (a właściwie rodzina ciągów {aν}ν) jest nieco zaburzonym ciągiem a=(a1,,an), jeśli istnieje stała K taka, że dla wszystkich dostatecznie małych ν zachodzi

|aν,jaj|Kν|aj|,1jn,

albo ogólniej

aνaKνa,

gdzie jest pewną normą w Rn. W pierwszym przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim o zaburzeniu w normie .

Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas

aνa=max1jn|aν,jaj|Kνmax1jn|aj|=Kνa,

i korzystając z faktu, że w przestrzeni skończenie wymiarowej wszystkie normy są równoważne, otrzymujemy dla pewnych stałych K1 i K2

aνaK1aνaK1KνaK2K1Kνa,

czyli nierówność dla zaburzenia w normie, ze stałą K=K2K1K.

Definicja Algorytm numerycznie poprawny

Algorytm 𝐀𝐋𝐆 rozwiązywania zadania nazywamy numerycznie poprawnym w zbiorze danych F0F, jeśli dla każdej danej fF0 wynik flν(𝐀𝐋𝐆(f)) działania algorytmu w arytmetyce flν można zinterpretować jako nieco zaburzony wynik algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji yν=(N(f))νN(F) o f, przy czym poziom zaburzeń nie zależy od f.

Formalnie znaczy to, że istnieją stałe K1, K2 oraz ν0>0 takie, że spełniony jest następujący warunek. Dla dowolnej νν0 oraz informacji yN(F0) można dobrać yνN(F) oraz (φ(yν))ν takie, że

yνyK1νy,
(φ(yν))νφ(yν)K2νφ(yν),

oraz

flν(𝐀𝐋𝐆(f))=flν(φ(N(f)))=(φ(yν))ν.
Numerycznie poprawny algorytm daje w arytmetyce flν wynik ALG(N(x)), który daje się zinterpretować jako mało zaburzony wynik f(y) zadania na mało zaburzonych danych x.

Zauważmy,że jeśli fRn, N(f)=(f1,,fn), oraz algorytm jest dokładny, 𝐀𝐋𝐆φS, to numeryczną poprawność algorytmu można równoważnie zapisać jako

flν(𝐀𝐋𝐆(f))=(S(fν))ν.

Numeryczna poprawność algorytmu jest wielce pożądaną cechą.

Algorytm numerycznie poprawny działa w arytmetyce zmiennoprzecinkowej (niemal) tak, jakby wszystkie obliczenia były wykonywane w arytmetyce dokładnej, a tylko dane i wyniki podlegały reprezentacji w skończonej precyzji.

Rola uwarunkowania zadania

Niech 𝐀𝐋𝐆()=φ(N()) będzie algorytmem numerycznie poprawnym dla danych F0F. Wtedy jego błąd w flν można oszacować następująco:

S(f)flν(𝐀𝐋𝐆(f))=S(f)(φ(yν))νS(f)φ(y)+φ(y)φ(yν)+φ(yν)(φ(yν))νS(f)𝐀𝐋𝐆(f)+φ(y)φ(yν)+K2νφ(yν)S(f)𝐀𝐋𝐆(f)+(1+K2ν)φ(y)φ(yν)+K2νφ(y),

przy czym yνyK1νy. Stąd w szczególności wynika, że jeśli algorytm jest numerycznie poprawny i ciągły ze względu na informację y, to

limν0S(f)flν(𝐀𝐋𝐆(f))=S(f)𝐀𝐋𝐆(f).

To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie się zachowywał w flν prawie tak, jak w arytmetyce idealnej.

Z powyższych wzorów wynika, że błąd w flν algorytmu numerycznie poprawnego zależy w dużym stopniu od:

  • dokładności algorytmu w arytmetyce idealnej,
  • dokładności ν arytmetyki flν,
  • wrażliwości algorytmu na małe względne zaburzenia informacji y.

Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy trochę więcej uwagi jedynie trzeciemu.

Jeśli φ spełnia warunek Lipschitza ze stałą L, a dokładniej

φ(yν)φ(y)Lyνy,

to

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)Lyνy+K2νφ(y)S(f)𝐀𝐋𝐆(f)+(1+K2ν)LK1νy+K2νφ(y).

W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny algorytmu proporcjonalnie do ν.

Bardziej jednak interesuje nas błąd względny. Wybierzmy "małe" η0 i przypuśćmy, że

φ(yν)φ(y)MK1νmax(η,φ(y)),

dla pewnej M niezależnej od y, tzn. błąd względny informacji, yνyK1νy, przenosi się na błąd względny wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia" M, albo na błąd bezwzględny ze współczynnikiem Mη. (Zauważmy, że gdybyśmy wzięli η=0, to dla y takiej, że φ(y)=0, musiałoby być φ(yν)=0 --- co zwykle, choć nie zawsze, nie jest prawdą.) Wtedy

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)MK1νmax(η,φ(y))+K2νφ(y)=S(f)𝐀𝐋𝐆(f)+ν(MK1(1+K2ν)+K2)max(η,φ(y)).

W szczególności, gdy algorytm jest dokładny i korzysta z pełnej informacji o f, tzn. S𝐀𝐋𝐆φ, to błąd

S(f)flν(𝐀𝐋𝐆(f))max(η,S(f))(MK1(1+K2ν)+K2)ν(MK1+K2)ν.

Stąd wynika, że jeśli (MK1+K2)ν1, to błąd względny algorytmu w flν jest mały, o ile S(f)η. Błąd względny jest proporcjonalny do dokładności ν, arytmetyki flν, współczynników proporcjonalności Ki algorytmu numerycznie poprawnego, oraz do wrażliwości M zadania S 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.

Przykład: Iloczyn skalarny

Załóżmy, że dla danych ciągów rzeczywistych o ustalonej długości n, aj, bj, 1jn, chcemy obliczyć

S(a,b)=j=1najbj

Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem i korzystający z pełnej informacji o kolejnych współrzednych.

Oznaczmy przez a~j i b~j reprezentacje liczb aj i bj w flν, a~j=aj(1+αj), b~j=bj(1+βj), oraz przez γj i δj błędy względne powstałe przy kolejnych mnożeniach i dodawaniach. Oczywiście |αj|,|βj|,|γj|,|δj|ν. Otrzymujemy

flν(j=1najbj)=((flν(j=1n1ajbj)+a~nb~n(1+γn))(1+δn)==((a~1b~1(1+γ1)+a~2b~2(1+γ2))(1+δ2)++a~nb~n(1+γn))(1+δn)=a~1b~1(1+γ1)(1+δ2)(1+δn)++a~jb~j(1+γj)(1+δj)(1+δn)=j=1najbj(1+ej),

gdzie w przybliżeniu (tzn. gdy ν0) mamy |e1|(n+2)ν i |ej|(nj+4)ν, 2jn. Algorytm naturalny jest więc numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany w flν można zinterpretować jako dokładny wynik dla danych aν,j=aj i bν,j=bj(1+ej), przy czym bνbp(n+2)νbp.

Zobaczmy teraz, jak błąd we współrzędnych bj wpływa na błąd wyniku. Mamy

|j=1najbjflν(j=1najbj)|=|j=1najbjj=1najbj(1+ej)|=|j=1nejajbj|j=1n|ej||ajbj|(n+2)νj=1n|ajbj|.

Stąd dla η0

|j=1najbjflν(j=1najbj)|max(η,|j=1najbj|)Kη(n+2)ν,

gdzie

Kη=Kη(a,b)=j=1n|ajbj|max(η,|j=1najbj|).

Zauważmy, że jeśli iloczyny ajbj są wszystkie dodatnie albo wszystkie ujemne, to Kη=1, tzn. zadanie jest dobrze uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej nν. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile liczba n składników nie jest horrendalnie duża. W ogólności jednak Kη może być liczbą dowolnie dużą i wtedy nie możemy być pewni uzyskania dobrego wyniku w flν.


Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 2.3 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Warto także przejrzeć rozdział 2 w

  • P. Deulfhard, A. Hohmann, Numerical Analysis in Modern Scientific Computing, Springer, 2003,

omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów. Nieocenioną monografią na ten temat jest

  • N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002.