Matematyka Obliczeniowa II

semestr zimowy 2012-13

Skrypt

Plan wykładu

1. Iteracyjne metody rozwiązywania algebraicznych układów równań-
   • liniowych
     - skrypt rozdział: 5: (bez sporej części 5.3.3 - tylko definicja SOR bez analizy) w szczególności - formaty macierzy rzadkich, metody iteracyjne stacjonarne : Jakobi, Gauss-Seidel, Richardsona, metody gradientowe np. najszybszego spadku,
     - rozdział: 6: metody Kryłowa - szczegółowo CG i GMRES z implementacją, analiza zbieżności CG (pełna) i elementy analizy zbieżności GMRES (jak w skrypcie).
     -rozdział 7: Prekonditionery - ogólnie prawostronne, lewostronne, obustronne, przykłady podstawowych najprostszych prekonditionerów np. Jakobi, blokowy Jakobi, SOR, Gaussa-Seidla, prekonditionery bazujące na prostych metodach iteracyjnych, niepełny LU czyli ILU, wzory na PCG.
     -rozdziały 8.2 i 8.3: Idea metody wielosiatkowej i met. strukturalnych.
   • nieliniowych
     - skrypt rozdział 9: - metoda iteracji prostej Banachowskiej i metoda Newtona z analizą zbieżności
     - rozdział 10: wariacje na temat metody Newtona (uproszczona met. Newtona, metoda Newtona z przybliżoną pochodna, niedokładna met. Newtona)
     - rozdział: 12.1: metoda nawrotów
2. Numeryczny problem własny - metody znajdowania przybliżeń wartości i wektorów własnych rozdziały 2,3,4 w szczególności - sprowadzenie macierzy do macierzy podobnej w postaci Hessenberga poprzez odbicia Householdera (trójdiagonalnej jeśli macierz symetryczna), metoda potęgowa, odwrotna potęgowa, metod ilorazów Rayleigha, metoda rónoczesnych iteracji i jej zwiazek z metoda QR, metoda QR "czysta" i z przesunięciami (przesunięcie Rayleigha i Wilkinsona) - idea zbieżności QR dla A=A^T; metoda dziel i rządź (divide and conquer) dla macierzy symetrycznej trójdiagonalnej; rozklad SVD - co to jest i do czego sie stosuje, jak go obliczać
   ta cz. wykladu nie bazowała szczegółowo na skrypcie ale poza met. dziel i rządź - skrypt zawiera wszystko co było na wykładzie
3. Całkowanie wielowymiarowe -- elementy rozdziałów 13, 14, 15 tzn. przekleństwa wymiaru dla wielowymiarowych całek (rozdz. 13) na przykładzie kwadratury wielowymiarowej prostokątów, metody Monte Carlo (MC) -definicja + tw. o zbieżności metody, wady i zalety (rozdz. 14), redukcja wariancji - idea metody warstw losowych i funkcji kontrolnej; idea metod Quasi Monte Carlo (QMC) - (rodz. 15)

Bibliografia

Podręczniki

1. James W. Demmel, Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia 1997.
2. Peter Deuflhard, Andreas Hohmann. Numerical analysis in modern scientific computing, wolu- men 43 serii Texts in Applied Mathematics. Springer-Verlag, New York, wydanie ii, 2003. An introduction. Ogólny podręcznik do analizy numerycznej - w szczególności zawiera opis metod bezpośrednich rozwiązywania układów równań liniowych, metody cg, metody Newtona i kilku metod dla zadania własnego.
3. J.M. Jankowscy, M. Dryja. Przegląd metod i algorytmów numerycznych, tom I i II. Biblioteka inżynierii oprogramowania. Wydawnictwo Naukowo-Techniczne, Warszawa, 1995. Ogólny podręcznik do analizy numeryczne. Zawiera opis bardzo wielu algorytmów nas interesujących - niestety część bez analizy
4. C. T. Kelley. Iterative methods for linear and nonlinear equations, wolumen 16 serii Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1995. Podręcznik zawierający opis zarówno metod iteracyjnych (w tym CG, PCG i GMRES) oraz wielowymiarowej metody Newtona
5. A. Kiełbasiński, H. Schwetlick. Numeryczna algebra liniowa. Wydawnictwa Naukowo-Techniczne, 1992. Podręcznik zawiera opis tylko metod bezpośrednich rozwiązuwania równań liniowych jak również niektóre algorytmy rozwiązywania zadania własnego
6. David Kincaid, Ward Cheney. Analiza numeryczna. Wydawnictwa Naukowo-Techniczne, Warszawa, 2006. Tłum. z ang.: S. Paszkowski. Ogólny podręcznik do analizy numeryczne.
7. J. Stoer, R. Bulirsch. Wstęp do analizy numerycznej. Biblioteka matematyczna. PWN, Warszawa, 1995.
8. Lloyd N. Trefethen, David Bau, III, Numerical linear algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.

Monografie lub bardzo zaawansowane podręczniki

1. John E. Dennis Jr., Robert B. Schnabel. Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall Series in Computational Mathematics. Prentice-Hall Inc., En- glewood Cliffs, N.J., 1983.
2. Peter Deuflhard. Newton methods for Nonlinear Problems. Affine Invariance and Adaptive Algori- thms. Springer International, 2002.
3. Eugene G. Dyakonov. Optimization in solving elliptic problems. CRC Press, Boca Raton, FL, 1996. Translated from the 1989 Russian original, Translation edited and with a preface by Steve McCormick.
4. Gene H. Golub, Charles F. Van Loan. Matrix computations. Johns Hopkins Studies in the Mathe- matical Sciences. Johns Hopkins University Press, Baltimore, MD, 3rd ed., 1996.
5. J. M. Ortega, W. C. Rheinboldt. Iterative solutions of nonlinear equations in several variables. Academic Press, New York, 1970.
6. Yousef Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied Ma- thematics, Philadelphia, PA, wydanie ii, 2003. On-line
7. Yousef Saad, Numerical Methods for Large Eigenvalue Problems. Society for Industrial and Applied Ma- thematics, Philadelphia, PA, wydanie ii, 2011. On-line
8. A. A. Samarski, J. S. Nikołajew. Metody rozwiązywania równań siatkowych. PWN 1988. Bardzo dobry podręcznik zawiera opis i dokładną analizę metod iteracyjnych z zastosowaniem do rozwiązywania układów równań liniowych powstałych z dyskretyzacji RRcz za pomocą metody różnic skończonych.
9. Barry F. Smith, Petter E. Bjorstad, William D. Gropp. Domain decomposition. Cambridge University Press, Cambridge, 1996. Parallel multilevel methods for elliptic partial differential equations.
10. Andrea Toselli, Olof Widlund. Domain decomposition methods - algorithms and theory, vol. 34 in Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2005.
11. J. F. Traub. Iterative Methods for the Solution of Equations. Englewood Cliffs, New York, 1964.

Projekt komputerowy

1. Metody iteracyjne dla układów liniowych z metodami dla macierzy niesymetrycznych. Biblioteka w C z kilkoma metodami iteracyjnymi oraz testy tych metod, min. GMRes, CG dla A'Ax=A'Tf, minimalnych residuów, FOM. (A' - to macierz A transponowana) -zajety
2. PCG, GMRes i prekonditionery. Zaimplementować metody PCG i GMRes oraz kolekcję prekonditionerów: blokowy Jakobi, blokowy Gauss-Seidel bez i z zakładkami, Richardson, hybrydowe - kilka kroków jakiejś stacjonarnej metody iteracyjnej np. Richardsona, niepełny rozkład Choleskiego i inne. Testy tych prekonditionerów z PCG i GMRes.
3. Równania nieliniowe - optymalizacja. Implementacja kilku nieliniowych metod rozwiązywania równań nieliniowych lub zagadnień optymalizacyjnych + testy, m.in. nieliniowa metoda CG i metoda Newtona-CG dla zadania minimalizacyjnego.
4. Zagadnienie własne. Implementacja kilku metod rozwiązywania zagadnia własnego - metoda QR z przesunięciami różnego typu, metoda Jakobiego i metoda dziel i rządź.

LAB

Program labu i skrypty octave'a

• Lab 1 - wprowadzenie do octave'a i metody bezpośrednie rozwiązywania układów równań liniowych (rozkład LU)
  mo2basic.m -prosty skrypt z elementarnymi operacjami; przykładami pętli, funkcji etc
• Lab 2 (22 X 2012) - metody bezpośrednie rozwiązywania układów równań liniowych (rozkłady LU i L'L), macierze rzadkie (sparse), iteracyjne metoda Jakobiego
  • Sprawdź jak działają funkcje lu() i chol().Sprawdź przy pomocy chol() czy macierz 3diagonalna z -1 na pod i nad diagonali i zero na diagonali jest dodatnio określona. Porównajmy czas rozwiązania układu równań z macierzą trójdiagonalną z 2 na diagonali i -1 na pod i nad diagonalach przy pomocy / oraz funkcji lu() i chol(). Stwórz tę macierz jako rzadką przy pomocy polecenia sparse ale bez tworzenia wpierw macierzy gęstej. Znów porównaj czas rozwiązania takiego układu. Dla jakiego max N octave rozwiąże taki układ dla macierzy pełnej i rzadkiej?
  • Przetestuj przykłady 5.8 i 5.9, rozwiąż ćwiczenia 5.9, 5.10, 5.11 w rozdziale 5.3.1 ze skryptu.
  • skrypt z kodem octave'a do przykładu 5.8
  • skrypt z kodem octave'a do przykładu 5.9
• Lab 3 (5 XI 2012) - metoda Jakobi cd., metoda Gaussa-Seidla i SOR oraz metoda Richardsona. Metody projekcji - metoda najszybszego spadku i minimalnych residuów.
  • Przetestuj przykłady 5.10, 5.12, rozwiąż ćwiczenia 5.13 (testy met. Jakobi i G-S dla T_N jak w przykladzie 5.10) i 5.15 (testy dla met iter. x_k=Bx_{k-1}+g dla B=[a,1;0,a] dla róznych wartosci a np a=0.51 - narysuj wykres residuum) w rozdziale 5.3 ze skryptu:
    skrypt z kodem octave'a do przykładu 5.10 (testy metod Jacobiego i Gaussa-Seidela dla B+pI - B losowa; p parametr);
    skrypt z kodem octave'a do przykładu 5.12 (testy metod Jacobiego, Gaussa-Seidela i SOR dla T_N - T_N trójdiagonalna [-1,2,-1]; i roznych wartosci parametru \omega );
  • Zaimplementuj m. Richardsona np. modyfikująć funkcje ze skryptu do przykładu 5.8. Przetestuj dla macierzy A' A - A losowa, i T_N dla N=1000 i różnych wartości parametru m Richardsona oraz dla optymalnego parametru i różnych N=100,1000 etc Powtórz testy ale w normie energetycznej.
  • Przetestuj Richardsona dla T_N-2I i różnych wartości parametru z różnymi losowymi wektorami startowymi.
  • Przetestuj przykład 5.13 w rozdziale 5.5.1 ze skryptu:
    skrypt z kodem octave'a do przykładu 5.13 (testy metod Jacobiego, Gaussa-Seidela i najszybszego spadku dla B'B+pI - B losowa; p parametr)
  • Przetestuj metodę najszybszego spadku x_{k+1}=x_k+ p_k r_k dla p_k=(r_k' r_k)/(r_k' Ar_k) dla tych samych macierzach co w testach m. Richardsona. (Metoda najszybszego spadku jest zaimplementowana w przykładzie 5.13)
  • Zaimplementuuj i przetestuj metodę minimalnych residuów x_{k+1}=x_k+ p_k r_k dla p_k=(r_k' A' r_k)/(r_k' A' Ar_k) dla macierzy losowych, A' A dla A losowej, T_N i T_N+ (1/N)*T_f gdzie T_f - trójdiagonalna z zerami na głownej diagonali, jedynkami na naddiagonali i -1 na poddiagonali
• Lab 4 (19 XI 2012)
  Testy metody CG.
  • Przetestuj funkcję pcg(), tzn. przeczytaj help i uruchom metodę dla układu z macierzą [2,1;1,2] z losowym wektorem prawej strony (bez prekonditionera).
  • Uruchom skrypt z przykładu 6.1: (skrypt z kodem octave'a do przykładu 6.1) czyli Kontynuujemy przykład 5.13. Chcąc porównywać cztery metody: Jacobiego, SOR, metodę najszybszego spadku oraz sprzężonych gradientów, będziemy korzystać z macierzy A = B' B + pI, gdzie p > 0 oraz B jest losową macierzą rozrzedzoną. Zwiększanie parametru p nie tylko poprawia diagonalną dominację, ale także poprawia uwarunkowanie A . Jako parametr relaksacji dla SOR wybraliśmy (strzelając w ciemno) \omega= 1.3 .
  • Ćwiczenie 6.7. Sprawdź w przykładzie 6.1 rozdziale 6.1.1 ze skryptu:
    skrypt z kodem octave'a do przykładu 6.1 , czy faktycznie uwarunkowanie macierzy A wpływa na szybkość zbieżności metody CG i najszybszego spadku. Aby zbadać uwarunkowanie macierzy, możesz skorzystać z polecenia cond(A), albo wykorzystać estymator uwarunkowania dostępny w pcg.
  • Ćwiczenie 6.8. Sprawdź, modyfikując kod przykładu 6.1, czy jeśli A nie będzie symetryczna (lub nie będzie dodatnio określona), wpłynie to istotnie na szybkość zbieżności metody CG i najszybszego spadku. Wypróbuj m.in. A = B + pI dla p > 0 (brak symetrii) tak dobranego, by A_{sym} > 0 oraz A = B' B + pI dla p takiego, żeby A miało i dodatnie, i ujemne wartości własne.
  • Uruchom skrypt z przykładu 6.2 ze skryptu:
    skrypt z kodem octave'a do przykładu 6.2 , Chcąc porównywać cztery metody: Jacobiego, SOR, metodę najszybszego spadku oraz sprzężonych gradientów dla macierzy jednowymiarowego laplasjanu T_N . Jako parametr relaksacji dla SOR wybraliśmy wartość optymalną, zgodnie z przykładem 5.11.
  • Zmodyfikuj kod przykładu 6.2 aby przetestować CG dla T_N+pI dla p=0,5,10,20,50 for N=25 and 50 .
  • Przetestuj cg dla macierzy P_N - zdeskretyzowanego laplacianu na kwadracie jednostkowym przy pomocy MRD na jednorodnej siatce: skrypt tworzacy P_N
• Lab 5 (19XI2012) Testy metody GMRES i macierzy sciskajacejtypu niepelny Choleski z progiem (czyli jak ILUT).
  • Załaduj m-plik: gmres() . Przetestuj funkcję gmres(), tj. rozwiąż układ z macierzą [2,1;3,2] z losową prawą stroną
  • GMRES dla dodatnio określonych i symetrycznych. macierzy. Zmodyfikuj skrypt z przykładu 6.1 dodając testy na gmres(). (link do pliku powyżej patrz Lab 4) Porównujemy 5 metod: Jakobi, SOR, najszybszego spadku, cg i gmres dla A = B' B + p I, p > 0 , B losowa. Parametr SOR = \omega= 1.3 .
  • Zmodyfikuj skrypt z przykładu 6.1 dodając testy na gmres(). Weź losową B - 25X25 i testuj Ax=b dla:
    1. A = B + pI p > 0
    2. A = B'B-pI for p>0
    3. A - macierz o wielokrotnych wartosciach wlasnych, np. wymiaru siedem z 3 krotna wartoscia wlasna 2 (pozostale 1 krotne) - dla te jmacierzy tez mozna przetestowac CG.
    4. A - ortogonalna np. znajdź rozkład QR macierzy losowej B i weź A=Q
  • Zmodyfikuj skrypt z przykładu 6.2:
    skrypt z kodem do przykład 6.2. , Porównujemy 5 metod: Jakobi, SOR, najszybszego spadku, cg i gmres dla T_N +pI (T_N - 1D laplacian). Parametr SOR - optymalny por. przykład 1.11. Testuj dla p=0,5,10,20,50 oraz N=25 i 50.
  • Zmodyfikuj skrypt z przykładu 6.2: aby testować GMRES dla T_N + (1/N)T_f + p I; p=0,5,10,20,50 dla N=25 i 50, gdzie T_f - trójdiagonalna z zerami na głównej diagonali, jedynkami na nad-diagonali i -1 na pod-diagonali.
  • Przetestuj GMRES dla macierzy T_N+(1/N)*T_f jak w testach metody MR (minimalnych residuów) patrz lab 3 - ostatnie zadanie.
  • Przetestuj gmres dla macierzy P_N - zdyskretyzowanego laplacianu na kwadracie jednostkowym przy pomocy MRD na jednorodnej siatce: skrypt tworzący P_N
  • Przetestuj macierze sciskające obustronnie typu ilut (niepełny rozkład Choleskiego) dla macierzy losowej symetrycznej dodatnio określonej rzadkiej jak w przykładzie 6.1. A=B'B+pI (p=0,1), tzn. wygeneruj macierz A, czynnik niepełnego rozkładu ILUT (w wersji Choleskiego) L, następnie policz C+L^{-1}A(L')^{-1} (nie liczac oczywiscie m. odwrotnych) i uruchom metody Jakobi, SOR, najszybszego spadku i CG (cg z estymacją uwarunkowania) dla A i C. Przetestuj zbieżności metod dla różnych wartości progu t. Można też policzyć uwarunkowanie A i L^{-1}A(L')^{-1}. Oczywiście skrypt generujący L nie zadziała raczej dla dużych N.
    Macierz niepełnego rozkładu Choleskiego L (LL^T \approx A) można uzyskać poprzez skrypt: icholt.m
• Lab 6 (17 XII 2012) Testy metody wielosiatkowej
  • Metoda wielosiatkowa dla T_N
    • testuj 'smoother' - weź s losowy, policz b=T_Ns (N=100) i policz 3-4 iteracje m. Richardsona x_{k+1}=x_k -a(T_N x_k-b)z parametrem a=0.25 z losowym x0;
    • narysuj wykresy błędów: x0-s i x-s (x - otrzymane przez 2-3 iteracje met. Richardsona).
    • Używając eigs() policz wektory i wartości T_N - narysuj wykresy wektorów własnych dla najmniejszej i największej wartości własnych tj. Q(:,1) i Q(:,N).
    • Narysuj wykresy wartości własnych I- aT_N for a =1,0.5,0.25,0.125 i N=100. (wartosci wlasne albo z wzozu lub policzone numerycznie z poprzedniego podpunktu)
  • ExtensionCreate.m tworzącą macierz E wymiaru R^{N_0} X R^N (N_0>N). utwórz E N=100 i N_0 = 25 i narysuj wykresy kolumn E (pierwszej, 20tej etc)
  • Załaduj m-plik z funkcją: twogrid.m i testuj dla T_Nx=b z N_0=N/2. losowego b - narysuj wykresy norm residuów dla N=24,100,1000. Czy szybkość zbieżności (ilość iteracji) zalezy od N?
  • Testuj fsolve.m na x+y=1;x^2+y^2=1 i równania Allen Cahna -u''+\delta u(1-u^2)=f z u=0 na brzegu. Po dyskretyzacji MRS: h^{-2}T_N U + \delta U.*(1-U^3)=F - delta=1,1e-2,1e2 etc
  • Testuj Newtona dla x+y=1;x^2+y^2=1 i r. Allena Cahna n -u''+\delta u(1-u^2)=f z u=0 na brzegu. Mozna uzyc m-pliku newton.m
  • Testuj Newtona z Jakobianem aproksymowanym różnicami skończonymi dla x+y=1;x^2+y^2=1 i różnych wartości parametru dla różnic skończonych. Mozna uzyc m-pliku fdnewton.m
  • Testuj Newtona z Jakobianem aproksymowanym różnicami skończonymi dla Allena-Cahna.
• Lab 7 testy metody potęgowej, odwrotnej potęgowej, ilorazów Rayleigha i QR.
  • Testuj czy fsolve() znajdzie rozwiązanie (a,x) takie, że [Ax-ax;0.5*x'*x-1]=0 dla A=A^T e.g. A=[2,1;1;2]
  • Załaduj funkcję: powerr() . Testuj dla macierzy A=[2,1;1;2] i A^{-1}. Zmodyfikuj kod aby wyprowadzać na ekran iloraz Rayleigha r_k x_k oraz błedy |r_k-\lambda_1| and \|x_k-q_1\| (dla znanej pary własnej (\lambda_1,q_1)).
  • Powtórz testy dla niesymetrycznej A=CDC^{-1} z D=diag(3,2,1) i losową macierzą C.
  • Zmodyfikuj kod powerr.m pisząc m-plik: inverse.m z odwrotną metodą potęgową Testuj dla A=[2,1;1;2] i macierzy o wartościach własnych: M,2,1 dla M=3,5,10,100. Wydrukuj błedy dla wektora i wartości własnej.
  • Zmodyfikuj kod powerr.m pisząc m-plik: Rayleigh.m z metodą ilorazów Rayleigha Testuj dla A=[2,1;1;2] i macierzy o wartościach własnych: M,2,1 dla M=3,5,10,100. Wydrukuj błedy dla wektora i wartości własnej.
  • Zmodyfikuj kod powerr.m aby metoda zadziałała dla macierzy A posiadającej symetryczne względem zera największe co do modułu wartości własne np dla macierzy symetrycznej o wartościach własnych M,-M,2,1 dla M=3,5,10,100.
  • Zaimplementuj metodę QR z i bez przesunięć (shifts), za przesunięcia weź przesunięcia Rayleigha tzn (A_k)_{N N}. (Z przesunieciami do domu)

Kilka rozwiazan.