Układy równań liniowych - rozkłady typu LU i LL'

W tym rozdziale zapoznamy się z metodami służących do rozwiązywania układów równań liniowych przy pomocy uzyskiwaniu odpowiednich rozkładów macierzy typu $ LU $ i $ LL^T $ oraz obliczaniu macierzy odwrotnej.

Zapoznamy się z odpowiednimi funkcjami octave'a służącymi
znajdowaniu odpowiednich rozkładów, czy rozwiązywaniu układów równań liniowych z pomocą tych rozkładów.

Podstawowe funkcje i operatory octave'a związane z rozkładem LU to:

  • [L,U,P]=lu(A) - funkcja zwracająca czynniki rozkład $ LU $ nieosobliwej macierzy $ A $, czyli $ PA=LU $,
  • [R]=chol(A) - funkcja zwracająca czynnik rozkładu Choleskiego
    macierzy symetrycznej dodatnio określonej $ A=A^T>0 $, czyli $ A=R^TR $, (zazwyczaj w literaturze podaje się ten rozkład w terminach czynnika $ L=R^T $),
  • operator x=A\b - operator rozwiązujący układ równań $ Ax=b $ dla $ A $ macierzy nieosobliwej i $ b $ wektora prawej strony,
  • inv(A) - funkcja zwracająca macierz odwrotną do $ A $ macierzy nieosobliwej,
  • cond(A) - funkcja zwracająca współczynnik uwarunkowania macierzy $ A $, czyli $ \|A\|_2*\|A^{-1}\|_2 $.
  1. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie.

    Przetestuj operator octave'a \ dla układu równań z macierzą $ A=[1,2;3,4] $ i $ x=[1;3] $
    z f=A*x. Czy y=A\f jest rozwiązaniem tego układu?

    Policz normy residualną $ \|Ay-f\|_p $ i normę błędu $ \|x-x\|_p $ dla $ p=1,2 $.

    Powtórz testy dla macierzy osobliwej $ [1,2;3,6] $ i prawie osobliwej $ [1,2;3,6-10^{-6}] $. Co się wtedy dzieje?

  2. Przetestuj solver octave'a, tzn operator backslash \ dla dwóch prostych układów równań liniowych z macierzami nieosobliwymi: $ A_1=[1,1;1,0.98] $ i $ b=[2;1] $, $ A_2=[1,1;1,0.99] $ i wektorem prawej strony układu równań liniowych $ b=[2;1] $ jak w poprzednim zadaniu. Policz
    w normie drugiej różnicę rozwiązań, normę Frobeniusa różnicy $ A_1-A_2 $ oraz uwarunkowania tych macierzy korzystając z funkcji octave'a
    cond().

    Powtórz zadanie dla macierzy: $ A_1=[1,1;1,1-2\epsilon] $ i $ b=[2;1] $, $ A_2=[1,1;1,1-\epsilon] $ dla $ \epsilon=10^{-p} $ z $ p=3,4,<br />
\ldots,16 $.

  3. Funkcja inv().

    Zapoznaj się z pomocą do funkcji: inv(). Przetestuj, obliczając macierz odwrotną do $ A=[1,1;1,-1] $. Czy macierz $ B $ obliczona za pomocą tej funkcji rzeczywiście jest macierzą odwrotną?

    Policz normy pierwszą i Frobeniusa $ \|A*B-I\| $ oraz $ \|B*A-I\| $.

    Zastosuj funkcje do rozwiązywania układu równań liniowych $ Ax=f $ dla znanego $ x $ (liczymy $ f=A*x $). Policz $ y $ jako iloczyn $ B $ z $ f $ i policz błędy residualny $ \|Ay-f\| $ i $ \|x-y\| $ w normie pierwszej i drugiej.

    Powtórz to zadanie dla macierzy $ N\times N $ losowej (funkcja octave'a rand() zwraca macierz losową) dla $ n=10,50,250,1250 $
    ze znanym rozwiązaniem - oszacuj czas przy pomocy funkcji tic i toc. Porównaj czas i błędy w normie pierwszej dla rozwiązania uzyskanego przy pomocy operatora octave'a \.

  4. Funkcje lu() i chol().
    • Zapoznaj się z pomocą do funkcji: lu() i chol().
    • Dla macierzy $ A=[1,1;1,-1] $ znajdź jej rozkłady: $ PA=LU $ i rozkład Choleskiego $ A=L_1^TL_1 $.
    • Sprawdź te rozkłady licząc normy macierzowe pierwszą i Frobeniusa błędów np. $ PA-LU $ i $ A-L_1^TL_1 $.
    • Zastosuj te rozkłady do znalezienia rozwiązania
      układu równań liniowych $ Ax=f $ ze znanym rozwiązaniem np. $ x=[1;1] $ i $ f=[2;0] $.

    • Policz normy pierwszą i drugą wektorowe pomiędzy $ x $, a wynikiem algorytmu $ w $, polegającym na zastosowaniu odpowiedniego rozkładu, oraz takie same normy residualne, tzn. normy $ Aw-f $.
  5. W octavie przetestuj eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy
    $ A=[e, \; 1; 1, \;1] $ z $ e=eps/10 $ (funkcja octave'a eps zwraca epsilon maszynowy)
    i wektorem prawej strony $ f=[1;0]^T $. Trzeba tu zaprogramować samodzielnie eliminację Gaussa bez wyboru elementu głównego, tzn. bez permutacji ani wierszy, ani kolumn dla macierzy $ 2\times 2 $.

  6. (#)
    Testy solvera octave'a dla macierzy Hilberta $ H(N) $

    Polecenie octave'a hilb() tworzy macierz Hilberta.

    • Stwórz macierz $ H(N) $ dla $ N=10,20,30 $,
    • Dla znanego rozwiązania, np. stałego równego jeden na każdej pozycji, czyli $ sol_k=1 $, policz $ H(N)*sol=f $,
    • Rozwiąż w octavie układ z $ H(N) $ i $ f $,
    • Policz normy (1,2 itd) $ \|H(N)x-f\| $ i $ \|x- sol\| $ dla $ x $ rozwiązania wyliczonego przez octave,
    • Policz uwarunkowanie macierzy Hilberta dla powyższych $ N $.
  7. Powtórz zadanie  [link] w arytmetyce pojedynczej precyzji.

    Należy tu sztucznie wymusić, za pomocą funkcji octave'a single(), aby zmienne były zmiennymi pojedynczej precyzji.

  8. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie  [link]).

    • Policz normy: macierzowa normę Frobeniusa i normę macierzową pierwszą różnicy macierzy $ E=H(N)*iH(N)-I $
      dla $ N=10,20,30 $, gdzie $ iH(N) $ macierz odwrotna do macierzy Hilberta obliczona przez invhilb().

    • Wykorzystaj $ iH(N) $ do rozwiązania układu równań z macierzą Hilberta $ H(N) $, tzn.:
      1. stwórz macierze $ H(N) $ oraz $ iH(N) $ dla $ N=10,20,30 $ korzystając z funkcji hilb() i invhilb(),
      2. dla znanego rozwiązania $ sol $, policz $ H(N)*sol=f $,
      3. rozwiąż układ $ H(N)x=f $ mnożąc $ f $ przez $ iH(N) $, tzn. $ x=iH(N)*f $,
      4. policz normy $ \|H(N)x-f\|_p $ i $ \|x- sol\|_p $ dla $ p=1,2,\infty $.
  9. (#)
    Stwórz w octavie macierz trójdiagonalną $ A=(a_{i j})_{i,j=1}^n $ wymiaru $ n\times n $ z $ 2 $ na diagonali i $ -1 $ na pod- i nad diagonalą, tzn.

    $$<br />
a_{i, j}=\left\{\begin{array}{rcl}<br />
2 &  &i=j  \\<br />
-1 &\qquad&  |i-j|=1 \\<br />
0 & &|i-j|>1<br />
\end{array},<br />
\right. \qquad  i,j=1,\ldots,n.<br />
$$

    przy pomocy funkcji octave'a diag(), jak i wprost w pętli. Porównaj czas używając tic i toc dla $ n=10,100,1000 $.

    Policz uwarunkowanie macierzy dla różnych $ N $.

    Policz macierz odwrotną przy pomocy inv() (czy jest trójdiagonalna?).

  10. Sprawdź z wykorzystaniem funkcji octave'a chol(), czy macierz z zadania  [link] jest symetryczna i dodatnio określona.
  11. (#)
    Zaimplementuj rozkład $ LU $ dla macierzy trójdiagonalnej $ N\times N $ bez wyboru elementu głównego, tzn. stwórz własną funkcję octave'a:

    function [x,d,l]=lu3diag(a,b,c,f,N)

    Parametry funkcji:

    • $ a,b,c $ - przekątna, podprzekątna i nadprzekątna macierzy $ A $,
    • $ f $ - wektor prawej strony,
    • $ N $ - wymiar zadania - długość przekątnej $ a $.

    Funkcja powinna zwrócić $ x $ - rozwiązanie $ Ax=f $ oraz
    czynniki rozkładu macierzy $ A=LU $, czyli diagonalę macierzy górnotrójkątnej $ U $ w wektorze $ d $ i poddiagonalę $ L $ - macierzy dolnotrójkątnej, trójdiagonalnej z jedynkami na diagonali, czyli wektor $ l $.

    Naddiagonala $ U $ równa się nadiagonali $ A $, tzn. wektorowi $ c $.

  12. Zastosuj funkcję z zadania  [link] do uzyskania rozkładu $ LU $ macierzy z zadania  [link] dla $ N=4,10,15 $.
    Porównaj z wynikiem uzyskanym przy pomocy funkcji octave'a
    lu().

  13. (#)
    Zastosuj funkcję z zadania  [link] do rozwiązania układu równań z macierzą z zadania  [link] dla $ N=10,100,1000 $, czy $ 10000 $:

    • Rozwiąż układ dla znanego rozwiązania, np. $ x=(1,\ldots,1)^T $,

    • Porównaj czas rozwiązywania własnym algorytmem i algorytmem octave'a, czyli operatorem A\f,

    • Policz błąd rezydualny i
      błąd rzeczywisty w normie drugiej, tzn. $ \|x-\tilde{x}\|_2 $ i $ \|f-A*\tilde{x}\|_2 $, gdzie $ f=A*x $, a $ \tilde{x} $ to rozwiązanie obliczone
      przez octave'a lub obliczone własnym algorytmem.

  14. Zaimplementuj rozkład $ LU $ dla macierzy trójdiagonalnej z częściowym wyborem elementu głównego.

    Następnie testuj ten rozkład jak w zadaniu  [link].

    Zastosuj ten rozkład do rozwiązania układu równań z macierzą z poprzedniego zadania $ A $ i z macierzą $ A-1.5*I $.

  15. Zaprogramuj rozkład Choleskiego typu dla macierzy $ A=A^T>0 $ trójdiagonalnej, tzn.
    policz $ L $ dolnotrójkątną z jedną poddiagonalą (czyli reprezentowaną przez dwa wektory z przekątną i podprzekątną) taką, że $ A=L^TL $.

    Zastosuj do rozwiązania układu równań z zadania  [link].

    Testuj analogicznie do zadania  [link].

    Porównaj macierz $ L $ uzyskaną w ten sposób z macierzą uzyskaną przez
    chol() zastosowaną do tej samej macierzy.