|
|
Linia 1: |
Linia 1: |
|
| |
| <!--
| |
| Konwertowane z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/ pawlik1/latex2mediawiki.php
| |
| -->
| |
|
| |
| =Ćwiczenia: rozwiązywanie układów równań liniowych=
| |
|
| |
| <div style="margin-top:1em; padding-top,padding-bottom:1em;"> | | <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: Mnożyć przez odwrotność to nie zawsze to samo...</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"> | | <div class="exercise"> |
|
| |
|
| W Octave układ równań <math>\displaystyle Ax=b</math> rozwiązujemy korzystając z "operatora
| | Aby obliczyć <math>\displaystyle S(a,b)=a^2-b^2</math> można zastosować |
| rozwiązywania równania"
| | dwa algorytmy: <math>\displaystyle {\bf ALG}_1(a,b)=a*a-b*b</math> oraz <math>\displaystyle {\bf ALG}_2(a,b)=(a+b)*(a-b)</math>. |
| <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 = A \ b;
| | <math>\displaystyle rd_\nu(a)=a</math> i <math>\displaystyle 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 = inv(A)*b;
| |
| </pre></div> | |
|
| |
| Przedyskutuj, które podejście jest lepsze i dlaczego. Przeprowadź eksperymenty
| |
| numeryczne weryfikujące Twoją tezę.
| |
| </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">
| |
| Wyznaczenie <math>\displaystyle A^{-1}</math> praktycznie nigdy nie jest konieczne w obliczeniach, a
| |
| zwłaszcza do rozwiązania układu równań. Nie dość, że jest bardziej kosztowne, to
| |
| dodatkowo nie jest numerycznie poprawne.
| |
|
| |
| Przykładowy test mógłby wyglądać w poniższy sposób. Duże macierze losowe są
| |
| zazwyczaj dobrze uwarunkowane i najczęściej nieosobliwe.
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| N = 400; A = rand(N,N);x = rand(N,1); b = A*x;
| |
| tic; X = A \ b; tX = toc;
| |
| tic; Y = inv(A)*b; tY = toc;
| |
| [norm(X-x)/norm(x), norm(Y-x)/norm(x); tX, tY]
| |
| </pre></div>
| |
|
| |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
| |
|
| |
| ans =
| |
| 1.9861e-13 2.8687e-13
| |
| 2.3644e-01 4.0609e-01
| |
| </pre></div>
| |
|
| |
| Wyszło (z grubsza) dwa razy szybciej i dwa razy dokładniej.
| |
|
| |
| </div></div></div>
| |
|
| |
|
| <div style="margin-top:1em; padding-top,padding-bottom:1em;"> | | <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: Wybór elementu głównego jest naprawdę ważny!</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"> | | <div class="exercise"> |
|
| |
|
| Rozwiąż prosty układ równań
| | Pokazać, że naturalny algorytm obliczania cosinusa |
| | kąta między dwoma wektorami <math>\displaystyle a, b\inR^n</math>, |
|
| |
|
| <center><math>\displaystyle \aligned 10^{-18} x_1 + x_2 &= 1,\\ | | <center><math>\displaystyle \cos(a,b)\,=\,\frac{\sum_{j=1}^n a_jb_j} |
| x_1 + 2x_2 &= 4,
| | {\sqrt{\Big(\sum_{j=1}^n a_j^2\Big) |
| \endaligned</math></center>
| | \Big(\sum_{j=1}^n b_j^2\Big)}}, |
| | </math></center> |
|
| |
|
| czterema sposobami:
| | jest numerycznie poprawny. Oszacować błąd względny wyniku |
| * na kartce papieru (dokładnie!)
| | w <math>\displaystyle fl_\nu</math>. |
| * w komputerze, w arytmetyce podwójnej precyzji, korzystając z rozkładu LU <strong>bez wyboru elementu głównego</strong>
| |
| * jak poprzednio, ale <strong>z wyborem elementu głównego w kolumnie</strong>
| |
| * w LAPACKu lub Octave.
| |
|
| |
| Porównaj uzyskane rozwiązania i wyciągnij wnioski.
| |
| </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">
| |
| Rozwiązanie uzyskane bez wyboru elementu głównego jest bardzo mało precyzyjne <math>\displaystyle x_1</math> jest obarczone bardzo duzym błędem! Prosta analiza poazuje, że winę za to ponosi bardzo wielki współczynnik w macierzy <math>\displaystyle L</math>. Tak nie jest, gdy stosujemy wybór elementu: elementy <math>\displaystyle L</math> są zawsze mniejsze od 1. Dlatego inne rozwiązania numeryczne dają w przypadku naszej macierzy wynik, będący bardzo dobrym przybliżeniem wyniku dokładnego. ''A propos:'' Octave korzysta z <code>DGESV</code> LAPACKa, a <code>DGESV</code> korzysta z wyboru w kolumnie...
| |
| </div></div></div>
| |
|
| |
|
| <div style="margin-top:1em; padding-top,padding-bottom:1em;"> | | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
Linia 83: |
Linia 30: |
| <div class="exercise"> | | <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 |
| | <math>\displaystyle \|A x\|_2</math> dla danej macierzy <math>\displaystyle A\inR^{n\times n}</math> i wektora |
| | <math>\displaystyle x\inR^n</math> jest numerycznie poprawny. Dokładniej, |
|
| |
|
| <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>\displaystyle fl_\nu (\|A x\|_2)\,=\,(A+E) x, |
| <div style="font-size:smaller; background-color:#efe"> Pamiętaj, pętle w Octave wykonują się bardzo wolno! </div>
| | </math></center> |
| </div></div> | |
| | |
| Wykorzystaj go do napisania funkcji, która rozwiąże układ równań <math>\displaystyle Ax=b</math>.
| |
|
| |
|
| Przetestuj tę funkcję na kilku macierzach i porównaj czas jego działania z czasem
| | gdzie <math>\displaystyle \|E\|_2\leq 2(n+2)\sqrt n\nu\|A\|_2</math>. Ponadto, jeśli |
| wykonania operacji <code>x = A\b</code>.
| | <math>\displaystyle A</math> jest nieosobliwa to |
|
| |
|
| Spróbuj zastosować swój algorytm do kilku specjalnych macierzy:
| | <center><math>\displaystyle |fl_\nu(\|A x\|_2)-\|A x\|_2|\,\leq\,2(n+2)\sqrt{n}\,\nu\, |
| * Hilberta dużego wymiaru
| | \left(\|A\|_2\|A^{-1}\|_2\right)\,\|A x\|_2. |
| * 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] = 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 = size(A,1);
| |
| for k=1:N-1
| |
| if (A(k,k) == 0.0)
| |
| return;
| |
| end
| |
| for i=k+1:N #wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
| |
| A(i,k) = A(i,k)/A(k,k);
| |
| end
| |
| for j=k+1:N #modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math>
| |
| for i=k+1:N
| |
| A(i,j) -= A(i,k)*A(k,j);
| |
| end
| |
| end
| |
| end
| |
| L = tril(A,-1) + eye(N);
| |
| U = 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] = 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 = size(A,1);
| |
| for k=1:N-1
| |
| if (A(k,k) == 0.0)
| |
| return;
| |
| end
| |
| #wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
| |
| A(k+1:N,k) = 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) -= A(k+1:N,k)*A(k,k+1:N);
| |
| end
| |
| L = tril(A,-1) + eye(N);
| |
| U = triu(A);
| |
| end
| |
| </pre></div>
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
| [Kod testujący lufa(0 i lufb()]
| |
| N = 100;
| |
| A = rand(N) + 1000*eye(N);
| |
| | |
| format short e;
| |
| | |
| t = ta = tb = Inf;
| |
| | |
| K = 5;
| |
| | |
| for i = 0:K
| |
| tic; [L, U] = lu(A); t = min(t,toc);
| |
| tic; [La, Ua] = lufa(A); ta = min(ta,toc);
| |
| tic; [Lb, Ub] = lufb(A); tb = 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"> </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 = sol(L,U,b)
| |
| x = 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;">
| |
| <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: Czy Twoje programy działają <strong>naprawdę</strong> szybko?</span>
| |
| <div class="exercise">
| |
| | |
| Zrób zadanie poprzednie programując, niczym Zosia-Samosia, wszystko od początku do końca w czystym C. Porównaj czasy działania Twojego programu i programu wywołującego po prostu procedurę biblioteczną LAPACKa
| |
| <code>DGESV</code>, najlepiej wspartą dobrze podrasowanymi BLASami.
| |
| | |
| <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
| |
| <div style="font-size:smaller; background-color:#efe"> To zadanie to wyzwanie! </div>
| |
| </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, namawiamy Cię do zmierzenia się z LAPACKiem i ATLASem....
| |
| | |
| Prościutki program --- po prostu tłumaczącu algorytm z wykładu bezpośrednio
| |
| na C --- u nas spisał
| |
| się, dla macierzy wymiaru <math>\displaystyle N=1\ldots 1024</math>, następująco:
| |
| | |
| [[Image:MNdegsvtiming.png|thumb|450px|center|Prosty solver w C kontra LAPACK i LAPACK z ATLASem.
| |
| Zwróć uwagę na piki zwykłego kodu dla <math>\displaystyle N</math> będących potęgami dwójki]]
| |
| [[Image:MNdegsvtiminglogscale.png|thumb|450px|center|Skala logarytmiczna pozwala lepiej ocenić
| |
| różnice pomiędzy czasami wykonania. Prosty program w C jest około dziesięć razy
| |
| wolniejszy od kodu korzystającego z LAPACKa i około 50 razy wolniejszy od kodu
| |
| korzystającego z LAPACKa i ATLASa.]]
| |
| | |
| </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: 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ą
| |
| 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> | | </math></center> |
|
| |
|
| Układy równań z tą samą macierzą, ale ze zmieniającą się prawą stroną równania
| |
| powstają często przy rozwiązywaniu, np. równań różniczkowych cząstkowych,
| |
| 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">
| |
| <div style="font-size:smaller; background-color:#efe"> Wystarczy tylko <strong>jeden raz</strong> wyznaczyć rozkład LU macierzy <math>\displaystyle A</math> </div>
| |
| </div></div> | | </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 algorytmów podstawienia w przód i w tył, operujące od razu na wszystkich kolumnach macierzy <math>\displaystyle B</math>. 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;"> | | <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"> | | <div class="exercise"> |
|
| |
|
| Ważnym przykładem macierzy
| | Niech <math>\displaystyle {\bf ALG}</math> będzie algorytmem numerycznie |
| szczególnej postaci są macierze symetryczne i dodatnio określone.
| | poprawnym w zbiorze danych <math>\displaystyle f\in F_0</math>, przy czym dla małych <math>\displaystyle \nu</math>, |
| Są to macierze spełniające warunki:
| | <math>\displaystyle fl_\nu({\bf ALG}(f))=\varphi(y_\nu)</math>, gdzie <math>\displaystyle \|y_\nu-y\|\le K\nu\|y\|</math> |
| | | i <math>\displaystyle K</math> nie zależy od <math>\displaystyle \nu</math> i <math>\displaystyle f</math> (<math>\displaystyle y=N(f)</math>). Pokazać, że |
| <center><math>\displaystyle A=A^T \qquad \mbox{oraz} \qquad x^T A x\,>\,0,\quad\forall x\ne 0. | | w ogólności <math>\displaystyle {\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>\displaystyle K_1</math> taka, że dla małych <math>\displaystyle \nu</math> i dla dowolnej |
| Dla takich macierzy można nieco zmniejszyć koszt kombinatoryczny
| | <math>\displaystyle 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>\displaystyle |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>\displaystyle 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></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;"> | | <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"> | | <div class="exercise"> |
|
| |
|
| Bardzo rzadko w praktyce numerycznej zdarza się potrzeba obliczenia wartości
| | Podaj przykład funkcji <math>\displaystyle f</math>, której miejsce zerowe <math>\displaystyle 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></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"> | | <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>\displaystyle x^* = f^{-1}(0)</math>, to |
| | | <center><math>\displaystyle |
| <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>
| |
| | |
| 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
| |
| arytmetyce zmiennoprzecinkowej. Ponieważ
| |
| | |
| <center><math>\displaystyle \mbox{det} (cA) = c^N A,
| |
| </math></center> | | </math></center> |
|
| |
|
| to dość łatwo będzie trafić (przy większych niż średnich wartościach <math>\displaystyle N</math>) na wartości | | Znaczy to, że im bardziej płaska jest <math>\displaystyle f</math> w otoczeniu pierwiastka <math>\displaystyle x^*</math>, tym |
| wyznacznika, które będą albo potwornie duże, albo potwornie małe. Zresztą, zrób
| | bardziej nawet małe zaburzenia <math>\displaystyle f</math> mogą spowodować duże przemieszczenie jej |
| eksperyment:
| | miejsca zerowego. |
| | |
| Najpierw z losową macierzą <math>\displaystyle 100\times 100</math> i jej krotnościami...
| |
| | |
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre> | |
|
| |
| octave:26> A=rand(100);
| |
| octave:27> det(A)
| |
| ans = -7.5240e+24
| |
| octave:28> det(10*A)
| |
| ans = -7.5240e+124
| |
| octave:29> det(100*A)
| |
| ans = -7.5240e+224
| |
| octave:30> det(1000*A)
| |
| ans = -Inf
| |
| </pre></div>
| |
|
| |
| oraz z coraz większymi macierzami losowymi:
| |
|
| |
|
| <div class="output" style="background-color:#e0e8e8; padding:1em"><pre> | | Zauważ, że dla wielokrotnych miejsc zerowych, <math>\displaystyle \mbox{cond} _{abs} (f^{-1},0) = \infty</math>, |
| | | co zgadza się z intuicją, bo może być, że nawet minimalne zaburzenie <math>\displaystyle f</math> |
| octave:32> det(rand(100))
| | spowoduje, iż miejsc zerowych po prostu nie będzie... |
| ans = 4.3751e+25
| |
| octave:33> det(rand(400))
| |
| ans = -3.0348e+218
| |
| octave:34> det(rand(800))
| |
| ans = -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> | | </div></div></div> |