Matematyka Obliczeniowa II

semestr zimowy 2011-12

Plan wykładu

1. Iteracyjne metody rozwiązywania
   • układów równań liniowych - skrypt rozdziały: 1 ($1.3.3. tylko do tw 1.4), 2, 3, 4 (tylko $4.2 i $4.3)
   • algebraicznych układów równań nieliniowych - skrypt rozdziały: 5, 6 (tylko $6.2 i $6.3)
2. Numeryczny problem własny - metody znajdowania przybliżeń wartości i wektorów własnych - 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 QR "czysta" i z przesunięciami (przesunięcie Railegha i Wilkinsona) - idea zbieżności QR dla A=A^T; metoda dziel i rządź (divide and conquer), metoda Hymana
3. Całkowanie wielowymiarowe - niestety nie zdążylismy nic zrobić - w planie było omówienie przekleństwa wymiaru (dla wielowymiarowych całek), metod Monte Carlo (MC) i Quasi Monter Carlo (QMC)...

Skrypt

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. (zajęty) Metody iteracyjne dla układów liniowych z metodami dla macierzy niesymetrycznych. Biblioteka w C z kilkoma metodami iteracyjnymi oraz testy tych metod.
2. PCG GMRes i prekonditionery. Zaimplementować metody PCG i GMRes oraz kolekcję prekonditionerów: blokowy Jakobi, Gauss-Seidel bez iz zakładkami, Richardson, hybrydowe - kilka kroków jakiejś stacjonarnej metody iteracyjnej, niepełny rozkład Choleskiego i inne. Testy tych prekonditionerów z PCG i GMRes.
3. (zajęty) Równania nieliniowe optymalizacja. Implementacja kilku nieliniowych metod rozwiązywania równań nieliniowych lub zagadnień optymalizacyjnych + testy.
4. Zagadnienie własne. Implementacja kilku metod rozwiązywania zagadnia własnego - metoda QR z przesunięciami, metoda Jakobiego i metoda dziel i rządź.

LAB

Program labu i skrypty octave'a

• Lab 1 (10 X 2011) - 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 (24 X 2011) - 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 1.8 i 1.9, rozwiąż ćwiczenia 1.9, 1.10, 1.11 w rozdziale 1.3.1 ze skryptu.
  • skrypt z kodem octave'a do przykładu 1.8
  • skrypt z kodem octave'a do przykładu 1.9
• Lab 3 (14 XI 2011) - metoda Jakobi cd., metoda Gaussa-Seidla i SOR oraz metoda Richardsona. Metody projekcji - metoda najszybszego spadku i minimalnych residuów.
  • Przetestuj przykłady 1.10, 1.12, rozwiąż ćwiczenia 1.13 i 1.15 w rozdziale 1.3.1 ze skryptu:
    skrypt z kodem octave'a do przykładu 1.10 (testy metod Jacobiego i Gaussa-Seidela dla B+pI - B losowa; p parametr);
    skrypt z kodem octave'a do przykładu 1.12 (testy metod Jacobiego, Gaussa-Seidela i SOR dla T_N+pI - T_N trójdiagonalna [-1,2,-1]; p parametr);
  • Zaimplementuj m. Richardsona np. modyfikująć funkcje ze skryptu do przykładu 1.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 1.13 w rozdziale 1.5.1 ze skryptu:
    skrypt z kodem octave'a do przykładu 1.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 1.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 (28 XI 2011)
  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 2.1: Kontynuujemy przykład 1.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 2.7. Sprawdź w przykładzie 2.1 rozdziale 2.1.1 ze skryptu:
    skrypt z kodem octave'a do przykładu 2.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 \lstinline!cond(A)!, albo wykorzystać estymator uwarunkowania dostępny w pcg.
  • Ćwiczenie 2.8. Sprawdź, modyfikując kod przykładu 2.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 < 0 takiego, żeby A miało i dodatnie, i ujemne wartości własne.
  • Uruchom skrypt z przykładu 2.2 ze skryptu:
    skrypt z kodem octave'a do przykładu 2.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 1.11.
  • Zmodyfikuj kod przykładu 2.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 (Dec 12th,2011) Tests of multigrid
  • Metoda wielosiatkowa dla T_N - testuj 'smoother' - weź s losowy, policz b=T_Ns i policz 3-4 iteracje m. Richardsona 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). Narysuj wykresy wartości własnych I- aT_N for a =1,0.5,0.25,0.125 i N=100. 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).
  • ExtensionCreate.m tworzącą macierz E:R^{N_0}-->R^N (N_0>N). utrwó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?
• Lab 6 (2 stycznia 2012) Testy GMRES.
  • 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 2.1 dodając testy na gmres(). (link do pliku powyżej por 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 2.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 - ortogonalna np. znajdź rozkład QR macierzy losowej B i weź A=Q
  • Zmodyfikuj skrypt z przykładu 2.2:
    skrypt z kodem do przykład 2.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 2.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 P_N - zdyskretyzowanego laplacianu na kwadracie jednostkowym przy pomocy MRD na jednorodnej siatce: skrypt tworzący P_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$
  • Implement Newton with FD approximation of Jacobian.
  • Testuj Newtona dla $x+y=1;x^2+y^2=1$ i r. for Allena Cahna n $-u''+\delta u(1-u^2)=f$ with $u=0$ na brzegu.
  • 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.
  Tutaj w pliku tar (zgzipowanym) kilka skryptów: np. zmodyfikowany skrypt przykładu 2.1 z testami GMRES, metoda Newtona klasyczna i z paroksymacją jakobianu z testami rzędu zbieżności etc: lab6.tgz
• Lab 7 (January 16nd, 2012) testy metody potęgowej, odwrotnej potęgowej, ilorazów Rayleigha i QR.
  • Testuj czy fsolve() znajdzie rozwiązanie (x,a)|-->[Ax-ax;0.5*x'*x-1]=0 dla A=A^T e.g. A=[2,1;1;2]
  • Załąduj funkcję: power() . 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 power.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 power.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 power.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}.