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

Uwarunkowanie zadania i algorytmy numerycznie poprawne.

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

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

Ćwiczenie

Aby obliczyć $S (a, b) = a^{2} - b^{2}$ można zastosować dwa algorytmy: ${𝐀 𝐋 𝐆}_{1} (a, b) = a * a - b * b$ oraz ${𝐀 𝐋 𝐆}_{2} (a, b) = (a + b) * (a - b)$ . Pokazać, że oba algorytmy są numerycznie poprawne, ale drugi z nich wywołuje mniejszy błąd względny wyniku w przypadku, gdy $r d_{ν} (a) = a$ i $r d_{ν} (b) = b$ .

Rozwiązanie

Rzeczywiście, dla pierwszego algorytmu obliczony w $f l_{ν}$ wynik $\tilde{s}$ spełnia

\tilde{s} = f l_{ν} (f l_{ν} (a \cdot a) - f l_{ν} (b \cdot b)) = (a^{2} (1 + ϵ_{a}) - b^{2} (1 + ϵ_{b})) (1 + ϵ_{-})

,

gdzie $| ϵ_{a} |, | ϵ_{b} |, | ϵ_{-} | \leq ν$ . A więc jesteśmy w sytuacji, gdy --- jeśli tylko $| a | \approx | b |$ --- może nastąpić redukcja cyfr przy odejmowaniu...

Natomiast drugi algorytm w ogóle nie jest na to czuły,

\tilde{s} = f l_{ν} (f l_{ν} (a - b) \times f l_{ν} (a + b)) = ((a - b) (1 + ϵ_{-}) \cdot (a + b) (1 + ϵ_{+})) (1 + ϵ_{\times})

,

gdzie znów $| ϵ_{+} |, | ϵ_{-} |, | ϵ_{\times} | \leq ν$ . Ponieważ ostatecznie

\tilde{s} = (a - b) (a + b) (1 + ϵ_{-}) (1 + ϵ_{+}) (1 + ϵ_{\times}) = = (a^{2} - b^{2}) (1 + E)

,

gdzie $| E | \leq 3 ν$ , algorytm drugi będzie zawsze dawał wynik obarczony małym błędem względnym.

Zwróć uwagę na istotną rolę przyjętego założenia, że $a$ i $b$ są liczbami maszynowymi, reprezentowanymi dokładnie w $f l_{ν}$ . W praktyce obliczeniowej, najczęściej właśnie z takimi danymi będziemy się spotykać...

Ćwiczenie

Pokazać, że naturalny algorytm obliczania cosinusa kąta między dwoma wektorami $a, b \in R^{n}$ ,

\cos (a, b) = \frac{\sum_{j = 1}^{n} a_{j} b_{j}}{\sqrt{(\sum_{j = 1}^{n} a_{j}^{2}) (\sum_{j = 1}^{n} b_{j}^{2})}}

,

jest numerycznie poprawny. Oszacować błąd względny wyniku w $f l_{ν}$ .

Ćwiczenie

Podaj przykład funkcji $f$ , której miejsce zerowe $x^{*}$ ma wspólczynnik uwarunkowania

mały
duży

Rozwiązanie

Ponieważ nasze zadanie to wyznaczenie $x^{*} = f^{- 1} (0)$ , to

{cond}_{a b s} (f^{- 1}, 0) = \frac{1}{f^{'} (x^{*})}

Znaczy to, że im bardziej płaska jest $f$ w otoczeniu pierwiastka $x^{*}$ , tym bardziej nawet małe zaburzenia $f$ mogą spowodować duże przemieszczenie jej miejsca zerowego.

Zauważ, iż dla wielokrotnych miejsc zerowych, ${cond}_{a b s} (f^{- 1}, 0) = \infty$ . Zgadza się to z intuicją, bo może się zdarzyć, że nawet minimalne zaburzenie $f$ spowoduje, iż miejsc zerowych po prostu nie będzie...

MN04LAB: Różnice pomiędzy wersjami

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

Uwarunkowanie zadania i algorytmy numerycznie poprawne.

Menu nawigacyjne

Działania na stronie

Opcje strony

Narzędzia osobiste

Nawigacja

Szukaj

Narzędzia

@@ Linia 1: / Linia 1: @@
+<!--
+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
+-->
+=Uwarunkowanie zadania i algorytmy numerycznie poprawne.=
-==Rozwiązywanie układów równań liniowych==
+{{powrot |Metody numeryczne | do strony głównej
+przedmiotu <strong>Metody numeryczne</strong>}}
+<div class="mw-collapsible mw-made=collapsible mw-collapsed">
+Oglądaj wskazówki i rozwiązania __SHOWALL__<br>
+Ukryj wskazówki i rozwiązania __HIDEALL__
+</div>
 <div style="margin-top:1em; padding-top,padding-bottom:1em;">
@@ Linia 6: / Linia 19: @@
 <div class="exercise">
-W Octave układ równań <math>\displaystyle Ax=b</math> rozwiązujemy korzystając z "operatora
+Aby obliczyć <math>S(a,b)=a^2-b^2</math> można zastosować
-rozwiązywania równania", tzn.
+dwa algorytmy: <math>{\bf ALG}_1(a,b)=a*a-b*b</math> oraz <math>{\bf ALG}_2(a,b)=(a+b)*(a-b)</math>.
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
+Pokazać, że oba algorytmy są numerycznie poprawne, ale drugi
+z nich wywołuje mniejszy błąd względny wyniku w przypadku, gdy
-x <nowiki>=</nowiki> A \ b;
+<math>rd_\nu(a)=a</math> i <math>rd_\nu(b)=b</math>.
-</pre></div>
-W tym celu Octave oczywiście wewenętrznie wykorzystuje funkcję LAPACKa
-<code>DGESV</code>. Ale w Octave jest także funkcja <code>inv</code>, wyznaczająca macierz
-odwrotną, więc niektóre (nie najlepsze, oględnie mówiąc) podręczniki zalecają
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
-x <nowiki>=</nowiki> inv(A)*b;
-</pre></div>
-Przedyskutuj, które podejście jest lepsze i dlaczego. Przeprowadź eksperymenty
-numeryczne weryfikujące Twoją tezę.
 </div></div>
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-Wyznaczenie <math>\displaystyle A^{-1}</math> praktycznie nigdy nie jest konieczne w obliczeniach, a
+Rzeczywiście, dla pierwszego algorytmu obliczony w <math>fl_\nu</math> wynik <math>\tilde{s}</math> spełnia
-zwłaszcza do rozwiązania układu równań. Nie dość, że jest bardziej kosztowne, to
+<center><math>
-dodatkowo nie jest numerycznie poprawne.
+\tilde{s} = fl_\nu(fl_\nu(a\cdot a)-fl_\nu(b\cdot b)) = (a^2(1+\epsilon_a) - b^2(1+\epsilon_b))(1+\epsilon_{-})</math>,</center>
+gdzie <math>|\epsilon_{a}|, |\epsilon_{b}|, |\epsilon_{-}| \leq \nu</math>. A więc jesteśmy w sytuacji, gdy --- jeśli tylko <math>|a|\approx |b|</math> --- może nastąpić redukcja cyfr przy odejmowaniu...
-Przykładowy test mógłby wyglądać w poniższy sposób. Duże macierze losowe są
+Natomiast drugi algorytm w ogóle nie jest na to czuły,
-zazwyczaj dobrze uwarunkowane i najczęściej nieosobliwe.
+<center><math>
+\tilde{s} = fl_\nu(fl_\nu(a - b) \times fl_\nu(a + b)) = ((a-b)(1+\epsilon_{-}) \cdot (a+b)(1+\epsilon_{+}))(1+\epsilon_{\times})</math>,</center>
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
-N <nowiki>=</nowiki> 400; A <nowiki>=</nowiki> rand(N,N);x <nowiki>=</nowiki> rand(N,1); b <nowiki>=</nowiki> A*x;
+gdzie znów <math>|\epsilon_{+}|, |\epsilon_{-}|, |\epsilon_{\times}| \leq \nu</math>.
-tic; X <nowiki>=</nowiki> A \ b; tX <nowiki>=</nowiki> toc;
+Ponieważ ostatecznie
-tic; Y <nowiki>=</nowiki> inv(A)*b; tY <nowiki>=</nowiki> toc;
+<center><math>
-[norm(X-x)/norm(x),  norm(Y-x)/norm(x); tX, tY]
+\tilde{s} = (a-b)(a+b)(1+\epsilon_{-})(1+\epsilon_{+})(1+\epsilon_{\times})=
-</pre></div>
+= (a^2 - b^2)(1+E)</math>,</center>
-<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
+gdzie <math>|E| \leq 3\nu</math>, algorytm drugi będzie <strong>zawsze</strong> dawał wynik obarczony małym błędem względnym.
-ans <nowiki>=</nowiki>
+Zwróć uwagę na istotną rolę przyjętego założenia, że <math>a</math> i <math>b</math> są liczbami maszynowymi, reprezentowanymi dokładnie w <math>fl_\nu</math>. W praktyce obliczeniowej,  najczęściej właśnie z takimi danymi będziemy się spotykać...
-.9861e-13   2.8687e-13
-.3644e-01   4.0609e-01
-</pre></div>
-Wyszło (z grubsza) dwa razy szybciej i dwa razy dokładniej.
 </div></div></div>
@@ Linia 56: / Linia 53: @@
 <div class="exercise">
-Zapisz w Octave algorytm rozkładu LU macierzy (bez wyboru elementu głównego) działający ''in situ''.
+Pokazać, że naturalny algorytm obliczania cosinusa
+kąta między dwoma wektorami <math>a, b\in R^n</math>,
-<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>\cos(a,b)\,=\,\frac{\sum_{j=1}^n a_jb_j}
-<div style="font-size:smaller; background-color:#efe"> Pamiętaj, pętle w Octave wykonują się bardzo wolno! </div>
+      {\sqrt{\Big(\sum_{j=1}^n a_j^2\Big)
-</div></div>
+             \Big(\sum_{j=1}^n b_j^2\Big)}}</math>,</center>
-Wykorzystaj go do napisania funkcji, która rozwiąże układ równań <math>\displaystyle Ax=b</math>.
+jest numerycznie poprawny. Oszacować błąd względny wyniku
+w <math>fl_\nu</math>.
-Przetestuj tę funkcję na kilku macierzach i porównaj czas jego działania z czasem
-wykonania operacji <code>x <nowiki>=</nowiki> A\b</code>.
-Spróbuj zastosować swój algorytm do kilku specjalnych macierzy:
-* Hilberta dużego wymiaru
-* diagonalnej z jednym elementem bardzo małym (a nawet równym zero)
 </div></div>
-<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
+<!--
-Najprostszy wariant to po prostu przepisanie algorytmu z wykładu, który już
-nawet trochę wykorzystywał notację dwukropkową MATLABa:
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
-function [L,U] <nowiki>=</nowiki> lufa(A)
-# Na początku warto dodać trochę sprawdzania danych: czy A jest kwadratowa? a w ogóle,
-# czy jest macierzą? te linie pozostawiamy inwencji Czytelnika
-N <nowiki>=</nowiki> size(A,1);
-for k<nowiki>=</nowiki>1:N-1
-	if (A(k,k) <nowiki>=</nowiki><nowiki>=</nowiki> 0.0)
-		return;
-	end
-	for i<nowiki>=</nowiki>k+1:N #wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
-		A(i,k) <nowiki>=</nowiki> A(i,k)/A(k,k);
-	end
-	for j<nowiki>=</nowiki>k+1:N #modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math>
-		for i<nowiki>=</nowiki>k+1:N
-			A(i,j) -<nowiki>=</nowiki> A(i,k)*A(k,j);
-		end
-	end
-end
-L <nowiki>=</nowiki> tril(A,-1) + eye(N);
-U <nowiki>=</nowiki> triu(A);
-end
-</pre></div>
-Jednak w tej implementacji mamy do czynienia z nawet potójnie zagnieżdżonymi
-pętlami, co dla nawet średnich <math>\displaystyle N</math> może być powodem sporego spowolnienia
-algorytmu. Warto więc uniknąć pętli, odwołując się do instrukcji wektorowych i
-blokowych.
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
-function [L,U] <nowiki>=</nowiki> lufb(A)
-# Na początku warto dodać trochę sprawdzania danych: czy A jest kwadratowa? a w ogóle,
-# czy jest macierzą? te linie pozostawiamy inwencji Czytelnika
-N <nowiki>=</nowiki> size(A,1);
-for k<nowiki>=</nowiki>1:N-1
-	if (A(k,k) <nowiki>=</nowiki><nowiki>=</nowiki> 0.0)
-		return;
-	end
-#wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
-		A(k+1:N,k) <nowiki>=</nowiki> A(k+1:N,k)/A(k,k);
-#modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math>
-		A(k+1:N,k+1:N) -<nowiki>=</nowiki> A(k+1:N,k)*A(k,k+1:N);
-end
-L <nowiki>=</nowiki> tril(A,-1) + eye(N);
-U <nowiki>=</nowiki> triu(A);
-end
-</pre></div>
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
- [Kod testujący lufa(0 i lufb()]
-N <nowiki>=</nowiki> 100;
-A <nowiki>=</nowiki> rand(N) + 1000*eye(N);
-format short e;
-t <nowiki>=</nowiki> ta <nowiki>=</nowiki> tb <nowiki>=</nowiki> Inf;
-K <nowiki>=</nowiki> 5;
-for i <nowiki>=</nowiki> 0:K
-	tic; [L, U] <nowiki>=</nowiki> lu(A); t <nowiki>=</nowiki> min(t,toc);
-	tic; [La, Ua] <nowiki>=</nowiki> lufa(A); ta <nowiki>=</nowiki> min(ta,toc);
-	tic; [Lb, Ub] <nowiki>=</nowiki> lufb(A); tb <nowiki>=</nowiki> min(tb,toc);
-end
-disp("Czasy: lu(), lufa(), lufb()"); [t, ta, tb]
-disp("Błąd:");
-[norm(L*U - A,inf), norm(La*Ua - A,inf), norm(Lb*Ub - A,inf)]
-</pre></div>
-Dla macierzy wymiaru <math>\displaystyle N=100</math> dostaliśmy na Pentium 1.5GHz
-{| border=1
-|+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
-|-
-| Implementacja  ||  lu (oryginalne)  ||   lufa   ||  lufb
-|-
-| Czas (s)  ||  4.0436e-03   ||  1.5877e+01  ||    1.0767e-01
-|-
-| Błąd <math>\displaystyle ||LU-A||_\infty</math>  ||  3.4596e-13  ||   8.0747e-13  ||   8.0747e-13
-|}
-Zwróć uwagę na charakterystyczne fakty:
-* Usunięcie wewnętrznych pętli stukrotnie(!) przyspieszyło <code>lufb</code> w
-stosunku do <code>lufa</code>
-* Gotowa funkcja Octave wciąż jest stukrotnie(!) szybsza od mojej
-najszybszej funkcji.
-To jest typowe, gdy korzystamy ze środowisk typu Octave czy MATLAB. Czasem
-nawet lepiej napisać implementację wykonującą więcej obliczeń, ale korzystającą
-w większm stopniu z wbudowanych funkcji Octave, niż napisać samemu implementację
-liczącą mniej, ale korzystającą z licznych instrukcji <code>for...</code> i podobnych.
-Innym sposobem przyspieszenia działania funkcji Octave jest ich prekompilacja.
-To w znacznym stopniu usuwa problemy związanie z korzystaniem z pętli, ale
-wykracza poza zakres naszego kursu.
-Funkcja wykorzystująca gotowy rozkład LU do rozwiązania układu <math>\displaystyle Ax=b</math> w MATLABie miałaby postać
-<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
-function x <nowiki>=</nowiki> sol(L,U,b)
-x <nowiki>=</nowiki> U \ ( L \ b);
-end
-</pre></div>
-gdyż MATLAB prawidłowo rozpoznaje układy z macierzą trójkątną i stosuje szybszy
-algorytm. Niestety, do tej pory w Octave nie jest to możliwe; nie ma też
-funkcji, która rozwiązywałaby układ z macierzą trójkątną (wykorzystując np.
-LAPACKa). Zrób więc to i... wyślij do twórców Octave'a!
-</div></div></div>
 <div style="margin-top:1em; padding-top,padding-bottom:1em;">
@@ Linia 198: / Linia 70: @@
 <div class="exercise">
-Zrób zadanie poprzednie w C i porównaj z procedurą biblioteczną LAPACKa
+Pokazać, że naturalny algorytm obliczania
-<code>DGESV</code>, najlepiej wspartą dobrze podrasowanymi BLASami.
+<math>\|A x\|_2</math> dla danej macierzy <math>A\inR^{n\times n}</math> i wektora
-</div></div>
+<math>x\inR^n</math> jest numerycznie poprawny. Dokładniej,
-<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Układy równań z wieloma prawymi stronami</span>
-<div class="exercise">
-Podaj sposób taniego wyznaczenia rozwiązania sekwencji <math>\displaystyle k<N</math> układów równań z tą
+<center><math>fl_\nu (\|A x\|_2)\,=\,(A+E) x</math>,</center>
-samą macierzą <math>\displaystyle A\in R^{N\times N}</math> i różnymi prawymi stronami:
-<center><math>\displaystyle
-Ax_i = b_i, \qquad i = 1,\ldots,k.
-</math></center>
-Układy równań z tą samą macierzą, ale ze  zmieniającą się prawą stroną równania
+gdzie <math>\|E\|_2\leq 2(n+2)\sqrt n\nu\|A\|_2</math>. Ponadto, jeśli
-powstają często przy  rozwiązywaniu, np. równań różniczkowych cząstkowych,
+<math>A</math> jest nieosobliwa, to
-gdzie  prawa strona układu odpowiada zmieniającym się warunkom  brzegowym.
-<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(\|A x\|_2)-\|A x\|_2|\,\leq\,2(n+2)\sqrt{n}\,\nu\,
-<div style="font-size:smaller; background-color:#efe"> Wystarczy tylko ''jeden raz'' wyznaczyć rozkład LU macierzy <math>\displaystyle A</math> </div>
+   \left(\|A\|_2\|A^{-1}\|_2\right)\,\|A x\|_2</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">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-No tak, oczywiście,
-wystarczy
-{{algorytm|||
-<pre>
-Znajdź rozkład <math>\displaystyle PA = LU</math>;
-Utwórz macierz <math>\displaystyle B = [b_1,\ldots,b_k]</math>;
-Rozwiąż <math>\displaystyle LY = B</math>;
-Rozwiąż <math>\displaystyle UX = Y</math>;
-</pre>}}
-Rozwiązanie każdego z układów <math>\displaystyle LY = B</math> i <math>\displaystyle UX = Y</math> można przyspieszyć, tworząc
-warianty blokowe znanych algorytmów. Zaimplementowano je w LAPACKu, w tandemie
-funkcji <code>DGETRF</code> (rozkład) i <code>DGETRS</code> (rozwiązanie) oraz w
-<code>DGESV</code>, łączącej funkcjonalność obu.
-</div></div></div>
 <div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Metoda Cholesky'ego</span>
+<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span>
 <div class="exercise">
-Ważnym przykładem macierzy
+Niech <math>{\bf ALG}</math> będzie algorytmem numerycznie
-szczególnej postaci są macierze symetryczne i dodatnio określone.
+poprawnym w zbiorze danych <math>f\in F_0</math>, przy czym dla małych <math>\nu</math>,
-Są to macierze spełniające <math>\displaystyle A=A^T</math> oraz
+<math>fl_\nu({\bf ALG}(f))=\varphi(y_\nu)</math>, gdzie <math>\|y_\nu-y\|\le K\nu\|y\|</math>
+i <math>K</math> nie zależy od <math>\nu</math> i <math>f</math> (<math>y=N(f)</math>). Pokazać, że
-<center><math>\displaystyle x^T A x\,>\,0,\qquad\forall x\ne 0.
+w ogólności <math>{\bf ALG}</math> nie musi być "numerycznie poprawny po
-</math></center>
+współrzędnych", tzn. w ogólności nie istnieje bezwzględna
+stała <math>K_1</math> taka, że dla małych <math>\nu</math> i dla dowolnej
-Dla takich macierzy można nieco zmniejszyć koszt kombinatoryczny
+<math>f\in F_0</math>
-i zużycie pamięci przeprowadzając trochę inny rozkład na macierze trójkątne:
-tak, aby otrzymać
-rozkład
-<center><math>\displaystyle A\,=\,L\cdot D\cdot L^T
-</math></center>
-zamiast <math>\displaystyle PA=LU</math>, przy czym <math>\displaystyle L</math> jest tu jak zwykle macierzą
-trójkątną dolną z jedynkami na przekątnej, a <math>\displaystyle D</math> jest macierzą
-diagonalną z dodatnimi elementami na diagonali. Opracuj taki algorytm. W jego
-implementacji możesz porównywać się z procedurą LAPACKa <code>DPOSV</code>
-Inny wariant tego samego
-rozkładu to tak zwany rozkład Cholesky'ego--Banachiewicza, w którym przy tych
-samych założeniach na <math>\displaystyle A</math>, szukamy rozkładu wykorzystującego tylko jedną macierz
-trójkątną dolną:
-<center><math>\displaystyle A\,=\,\widetilde{L}\cdot \widetilde{L}^T,
+<center><math>|y_{\nu,j}-y_j|\,\le\,K_1\,\nu\,|y_j|, \qquad 1\le j\le n</math>,</center>
-</math></center>
-(oczywiście tym razem nie żądamy, aby <math>\displaystyle \widetilde{L}</math> miała na diagonali jedynki).
+gdzie <math>y=(y_1,\ldots,y_n)</math>.
-Jaka jest relacja między rozkładem <math>\displaystyle LDL^T</math> a <math>\displaystyle \widetilde{L}\widetilde{L}^T</math>?
 </div></div>
-<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
+-->
-Bez zmniejszenia ogólności rozpatrzymy
-tylko pierwszy krok rozkładu <math>\displaystyle LDL^T</math>. W tym celu, zauważmy najpierw, że
-<math>\displaystyle a_{1,1}= e_1^TA e_1>0</math> (gdzie <math>\displaystyle  e_1</math> jest pierwszym
-wersorem), a więc nie musimy przestawiać wierszy,
-bo element na diagonali jest niezerowy. W pierwszym kroku mnożymy
-macierz <math>\displaystyle A</math> z lewej strony przez odpowiednią macierz <math>\displaystyle L_1</math>, a potem
-z prawej przez <math>\displaystyle L_1^T</math>. Kluczem do zrozumienia algorytmu jest uwaga,
-że efektem mnożenia macierzy <math>\displaystyle L_1 A</math> z prawej strony przez <math>\displaystyle L_1^T</math>
-jest wyzerowanie elementów pierwszego wiersza poza <math>\displaystyle a_{1,1}</math>
-i pozostawienie niezmienionych pozostałych elementów. Ponadto
-macierz <math>\displaystyle A^{(1)}=L_1 AL_1^T</math> jest symetryczna i dodatnio określona.
-Rzeczywiście,
-<center><math>\displaystyle \Big(A^{(1)}\Big)^T\,=\,(L_1AL_1^T)^T
-     \,=\,(L_1^T)^TA^TL_1^T\,=\,L_1AL_1^T,
-</math></center>
-oraz dla <math>\displaystyle  x\ne 0</math>
-<center><math>\displaystyle x^T A^{(1)} x\,=\, x^T L_1 A L_1^T  x
-    \,=\,(L_1^T x)^T A(L_1^T x)\,>\,0,
-</math></center>
-bo <math>\displaystyle  x\ne 0</math> implikuje <math>\displaystyle L_1^T x\ne 0</math>. Stąd
-<math>\displaystyle a^{(1)}_{2,2}= e_2^TA^{(1)} e_2>0</math>. Postępując tak
-dalej otrzymujemy
-<center><math>\displaystyle L_{n-1}L_{n-2}\cdots L_2L_1AL_1^TL_2^T\cdots L_{n-2}^TL_{n-1}^T
-    \,=\,D,
-</math></center>
-przy czym macierz <math>\displaystyle D</math> jest diagonalna i dodatnio określona,
-a więc wyrazy na diagonali są dodatnie. Oznaczając
-<center><math>\displaystyle L\,=\,L_1^{-1}L_2^{-1}\cdots L_{n-1}^{-1}
-</math></center>
-dostajemy żądany rozkład.
-Zauważmy, że przy praktycznej realizacji rozkładu <math>\displaystyle A=LDL^T</math>
-wystarczy modyfikować jedynie wyrazy pod i na głównej przekątnej
-macierzy wyjściowej, ponieważ, jak zauważyliśmy, wszystkie kolejne
-macierze <math>\displaystyle A^{(k)}</math> są symetryczne. Pozwala to zmniejszyć --- w stosunku do
-rozkładu LU --- koszt
-obliczeniowy o połowę do <math>\displaystyle 2n^3/3</math> operacji arytmetycznych i zmieścić czynniki
-rozkładu w pamięci przeznaczonej na przechowywanie istotnych elementów <math>\displaystyle A</math>
-(czyli w dolnym trójkącie macierzy).
-Oczywiście, <math>\displaystyle \widetilde{L} = L \sqrt{D}</math>.
-</div></div></div>
 <div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Obliczanie wyznacznika macierzy</span>
+<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span>
 <div class="exercise">
-Bardzo rzadko w praktyce numerycznej zdarza się potrzeba obliczenia wartości
+Podaj przykład funkcji <math>f</math>, której miejsce zerowe <math>x^*</math> ma wspólczynnik
-wyznacznika macierzy <math>\displaystyle A \in R^{N \times N}</math>. Zaproponuj metodę obliczania
+uwarunkowania
-<math>\displaystyle  \mbox{det} (A)</math> oraz wskaż, jakiego rodzaju problemy numeryczne możesz napotkać.
+* mały
+* duży
 </div></div>
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-Wystarczy wykonać rozkład <math>\displaystyle PA=LU</math>. Wtedy
+Ponieważ nasze zadanie to wyznaczenie <math>x^* = f^{-1}(0)</math>, to
+<center><math>
-<center><math>\displaystyle  \mbox{det}  A\,=\,(-1)^p u_{11}u_{22}\cdots u_{nn},
+  \mbox{cond} _{abs} (f^{-1},0) = \frac{1}{f'(x^*)}</math></center>
-</math></center>
-gdzie <math>\displaystyle p</math> jest liczbą przestawień wierszy w eliminacji. Tak więc, koszt
-obliczenia wyznacznika to <math>\displaystyle O(N^3)</math>.
-Kłopoty numeryczne możemy mieć z reprezentacją samej wartości wyznacznika w
+Znaczy to, że im bardziej płaska jest <math>f</math> w otoczeniu pierwiastka <math>x^*</math>, tym
-arytmetyce zmiennoprzecinkowej. Ponieważ
+bardziej nawet małe zaburzenia <math>f</math> mogą spowodować duże przemieszczenie jej
+miejsca zerowego.
-<center><math>\displaystyle  \mbox{det}  (cA) = c^N A,
+Zauważ, iż dla wielokrotnych miejsc zerowych, <math>\mbox{cond} _{abs} (f^{-1},0) = \infty</math>. Zgadza się to z intuicją, bo może się zdarzyć, że nawet minimalne zaburzenie <math>f</math>
-</math></center>
+spowoduje, iż miejsc zerowych po prostu nie będzie...
-to dość łatwo będzie trafić (przy większych niż średnich wartościach <math>\displaystyle N</math>) na wartości
-wyznacznika, które będą albo potwornie duże, albo potwornie małe. Zresztą, zrób
-eksperyment:
-Najpierw z losową macierzą <math>\displaystyle 100\times 100</math> i jej krotnościami...
+[[Image:MNnonlinearcond2.png|thumb|550px|center|Gdy trochę zaburzymy wartości funkcji f, dobrze uwarunkowane miejsce zerowe nie przemieści się zbyt daleko od miejsca zerowego f.]]
-<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
+[[Image:MNnonlinearcond4.png|thumb|550px|center|Gdy trochę zaburzymy wartości funkcji f, źle  uwarunkowane miejsce zerowe może przemieścić się bardzo daleko od miejsca zerowego f.]]
-octave:26> A<nowiki>=</nowiki>rand(100);
-octave:27> det(A)
-ans <nowiki>=</nowiki> -7.5240e+24
-octave:28> det(10*A)
-ans <nowiki>=</nowiki> -7.5240e+124
-octave:29> det(100*A)
-ans <nowiki>=</nowiki> -7.5240e+224
-octave:30> det(1000*A)
-ans <nowiki>=</nowiki> -Inf
-</pre></div>
-oraz z coraz większymi macierzami losowymi:
-<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
-octave:32> det(rand(100))
-ans <nowiki>=</nowiki>  4.3751e+25
-octave:33> det(rand(400))
-ans <nowiki>=</nowiki> -3.0348e+218
-octave:34> det(rand(800))
-ans <nowiki>=</nowiki> -Inf
-</pre></div>
-Można sobie z
-tym poradzić, szacując wpierw rząd wielkości wyznacznika, a następnie
-reprezentując jego wartość w formie cecha--mantysa, np. jako parę liczb <math>\displaystyle m,c</math>
-takich, że <math>\displaystyle  \mbox{det} (A) = m\cdot 10^c</math>, ale
-dopasowując obie tak, by nie nastąpił ani nadmiar, ani niedomiar.
 </div></div></div>