WSTĘP DO OBLICZEŃ NAUKOWYCH

Wykład fakultatywny
Semestr letni 2000/2001

Nowości

Kolokwium na żywo z Octave'a odbędzie się w piątek, 25.05.01. Proszę koniecznie się nie spóźnić! Kolokwium będzie dotyczyć utrwalonej wiedzy o wykorzystaniu Octave, na pewno też trzeba będzie rozwiązać jakieś równanie nieliniowe.

Przykładowe rozwiązanie zagadnienia własnego metodą Newtona w Octavie. Uwaga: są tam nowe pomysły!

Terminy

Wykład: wtorek, 10-12 (2100)

Ćwiczenia/Lab: piątek, 8-10 (3150/Lab2)

Egzamin pisemny: wtorek, 05.06.01, 9-12, s.2070. Wyniki nastpnego dnia rano. Prócz wyników zostanie wywieszona lista dla chętnych na egzamin ustny. Proszę się na nią zawczasu wpisywać, aby uniknąć niepotrzebnych kolejek.

Egzamin ustny: sroda,06.06.01, 9-12, s.1020. Osoby, którym z ważnych powodów nie odpowiada termin, proszę o wcześniejszy kontakt, w celu ustalenia dogodnego terminu.

Spis treści

Opis wykładu

Wykład szczególnie polecany studentom, którzy w ramach pracy magisterskiej bądź zawodowej będą prowadzić rozmaite symulacje numeryczne. Na wykładzie zostaną omówione zaawansowane metody numeryczne rozwiązywania wybranych zagadnień matematyki stosowanej. Zasadniczy ciężar wykładu spocznie na teorii i praktyce współcześnie stosowanych metod iteracyjnych dla równań nieliniowych i liniowych. Zostaną omówione z dowodami różne warianty wielowymiarowej metody Newtona (m.in. klasyczna, z różnicową aproksymacją jakobianu, Broydena) dla układów równań nieliniowych, oraz sposoby globalizacji zbieżności. Ponadto, zostaną także przedstawione ważne metody iteracyjne (m.in. CG, GMRES) dla wielkich układów równań liniowych z macierzami rozrzedzonymi. Obliczenia naukowe wymagają zazwyczaj dużej mocy obliczeniowej, dlatego na wykładzie zajmiemy się także konsekwencjami pewnych detali architektury komputera (np. pamięć hierarchiczna, superskalarność i wieloprocesorowość) dla efektywności algorytmów numerycznych. Istotnym punktem zajęć laboratoryjnych będzie nauczenie słuchaczy skutecznego wykorzystania oprogramowania naukowego.

Wymagania wstępne

Dla dobrego zrozumienia wykładanego materiału wymagana będzie wiedza z następujących przedmiotów podstawowych:

Ponadto, dużo zależeć będzie od zaangażowania i chęci poznania słuchaczy.

Zasady zaliczenia

Do egzaminu w pierwszym terminie mogą przystąpić osoby, które zaliczyły ćwiczenia. Drugi termin jest dla tych (niewyrejestrowanych) studentów, którzy nie zaliczyli ćwiczeń, lub nie zdali egzaminu w I terminie, lub chcą poprawić ocenę uzyskaną w pierwszym terminie. Podczas egzaminu będzie można korzystać z własnych notatek zebranych na jednej kartce formatu A4. Egzamin będzie miał formę pisemną, z ewentualną poprawką na ustnym, gdy są uzasadnione nadzieje na lepszą ocenę.

Na egzaminie pojawią się zadania typu robionych na ćwiczeniach, a także pytania natury bardziej teoretycznej (typu pytań ``na egzamin ustny''). Przykładowo, egzamin mógłby składać się z następujących pytań:

1.: Zbadaj rząd zbieżności metody X.
2.: Udowodnij własność Y metody X lub podaj przykład, gdy tej własności nie ma.
3.: Porównaj metodę X z metodą X^*.
4.: Co to jest pamięć podręczna cache? Jakie są konsekwencje stosowania cache'u (pozytywne i negatywne) dla implementacji algorytmów numerycznych?
5.: Wyjaśnij, na czym polega backtracking.

Zadania egzaminacyjne z r.akad. 1999/2000

Zadania egzaminacyjne z r.akad. 2000/2001

Część ćwiczeń odbędzie się w Laboratorium Komputerowym. Zgodnie z powszechnie panującymi regułami, obecność na ćwiczeniach jest obowiązkowa. Na ćwiczeniach, ocenie podlegać będą prace domowe i praca laboratoryjna. W roku 2000/2001 będzie to kolokwium na żywo z Octave'a. W roku 1999/2000 był to projekt zaliczeniowy.

Ocena ostateczna będzie wypadkową z egzaminu i z ćwiczeń.

Literatura

1.: Kelley, C.T., Iterative methods for linear and nonlinear equations, SIAM 1995
2.: Golub, G; Ortega, J.M. Scientific computing. An introduction with parallel computing, Academic Press, 1993
3.: Ortega, J. M.; Rheinboldt, W. C. Iterative solutions of nonlinear equations in several variables, Academic Press, 1970
4.: MathWorks, Inc., MATLAB User's Guide
5.: Fosdick et al, Introduction to High Performance Scientific Computing
6.: Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, and Henk Van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, 1994

Linki

Netlib: mnóstwo oprogramowania numerycznego, dokumentacja i artykuły, także rozmaite bazy danych, m.in. Top 500.
Tygodnik numeryczny NA Digest - warto zaprenumerować!
Darmowe klony Matlaba: Octave (łudząco podobny do Matlaba) oraz Scilab (mniej kompatybilny z Matlabem, ale bardziej rozbudowany niż Octave).
Podręcznik Matlaba dla początkujących (po angielsku).
Lokalnie dostępna z Laboratorium Studenckiego dokumentacja Octave: /usr/local/share/doc/octave-2.0.16/octave_toc.html
Lokalnie dostępna z Laboratorium Studenckiego dokumentacja Vigie: /usr/local/share/doc/vigie/VigieUserGuide.ps.gz (gzipped postscript, mozna ogladac po ściągnięciu przez gunzip VigieUserGuide.ps.gz; ghostview VigieUserGuide.ps) Pliki demonstracyjne Vigie znajdują się w Labie Studenckim w katalogu /usr/local/share/vigie.

Konspekt wykładu 1999/2000

W1, 15.02.00
Rejestracja. Czym zajmują się OBLICZENIA NAUKOWE? Wymagania/możliwości algorytmiczne i sprzętowe; trochę historii.

Lab
Wprowadzenie do użytkowania systemu Unix.

W2, 22.02.00
Wpływ hardware'u na algorytmy numeryczne. Stacje robocze RISC. Specyfika procesorów RISC. Superskalarność, pipelining i chaining. Procesory wektorowe. Superkomputery: SMP, PVP, MPP. Uwaga o PVM i MPI. Zasoby komputerowe dostępne w Polsce i na świecie (lista 500 najsilniejszych komputerów). Prawo Amdahl'a w wersji sekwencyjnej i równoległej.

Ćw
Wpływ hardware'u na algorytmy numeryczne - ważne uzupełnienia do wykładu. Pamięć hierarchiczna. Lokalność odwołań do danych w czasie i przestrzeni. Pamięć podręczna: rodzaje, specyfika, zalety i niebezpieczeństwa. Algorytmy efektywnie współpracujące z pamięcią hierarchiczną na przykładzie mnożenia dwóch macierzy. KONKURS!

W3, 29.02.00
Środowisko obliczeń numerycznych Matlab na przykładzie pakietu Octave: przyjazny użytkownikowi interfejs do profesjonalnych bibliotek numerycznych LAPACK i LINPACK (algebra liniowa), ODEPACK (równania różniczkowe zwyczajne, także równania rózniczkowo-algebraiczne), MINPACK (równania nieliniowe), QUADPACK (kwadratury) oraz FFTPACK (szybkie przekształcenia Fouriera). Tu można obejrzeć przykładowy zapis sesji Octave.

Lab
Dalsze informacje o Octave: więcej o algebrze liniowej, proste operacje wejścia-wyjścia.

W4, 07.03.00
Wybrane pakiety numeryczne i ich wykorzystanie w praktyce numerycznej. Wizualizacja danych numerycznych.

Ćw
Przypomnienie jednowymiarowej metody Newtona. Praca domowa.

W5, 14.03.00
Przykłady zagadnień matematycznych prowadzących do układu równań nieliniowych. Równania nieliniowe skalarne: przypomnienie. Definicja rzędu zbieżności. Twierdzenie o kwadratowej zbieżności metody Newtona dla f' Lipschitzowskiej. Twierdzenie o rzędzie zbieżności metody Newtona i jego zastosowanie: przypadek zbieżności tylko liniowej oraz przypadek zbieżności sześciennej.

Lab
Wizualizacja funkcji 1- i 2-wymiarowych w Octave. Praca z plikami w Octave. Format danych dla Vigie. Skrypt Octave'a eksportujący pole wektorowe do Vigie.
Praca domowa.

W6, 21.03.00
Podstawowe fakty dotyczące funkcji wielu zmiennych. Twierdzenia o wartości średniej w różnych wariantach. Założenia standardowe. Twierdzenie o lokalnym oszacowaniu pochodnej, jej odwrotności i funkcji.

Ćw
Różne zadania wokół metody Newtona dla równań skalarnych.

W7, 28.03.00
Twierdzenie o zbieżności wielowymiarowej metody Newtona. Kryteria stopu metody Newtona. Implementacja wielowymiarowej metody Newtona i dyskusja kosztu. Metoda cięciw, jej implementacja i koszt.

Lab
Sprawdzenie pracy domowej. Wizualizacja normy pola wektorowego w Vigie. Test jednowymiarowej metody Newtona. Modyfikacje jednowymiarowej metody Newtona: m.in. historia zbieżności.
Praca domowa

W8, 04.04.00
Twierdzenie o lokalnej stabilności metody Newtona. Twierdzenie o zbieżności metody cięciw. Twierdzenie o zbieżności metody Newtona z przybliżoną odwrotnością pochodnej (bez dowodu). Aproksymacja F'(x) przez ilorazy różnicowe: motywacja (koszt) i idea wyboru długości kroku w warunkach niedokładnego obliczania wartości F(x).

Ćw
Zbieżność metody Steffensena. Liniowa zamiana zmiennych w metodzie Newtona: analiza, implementacja i motywy. Dowód twierdzenia o otwartości zbioru macierzy nieosobliwych.
Praca domowa

W9, 11.04.00
Lemat o aproksymacji F'(x) ilorazami różnicowymi. Dyskusja wyboru długości kroku h w obecności szumu. Operator różnicowej aproksymacji pochodnej. Koszt. Twierdzenie o zbieżności metody Newtona z różnicową aproksymacją pochodnej (przypadek dokładnego F).

Podstawowe idee metod globalizacji zbieżności metody Newtona: przykład motywujący. Prototyp algorytmu tropienia (ang. backtracking).

Lab
Krótki sprawdzian znajomości Octave. Dalsze testy metody Newtona. Zera prawie-wielokrotne.
Praca domowa.

W10, 18.04.00
Algorytm backtracking z dostateczną redukcją residuum. Metody wielomianowe: dwu- i trzypunktowe.

Heurystyka metody Broydena i wyprowadzenie wzoru. Twierdzenie o zbieżności metody Broydena (bez dowodu).

Ćw
Metoda Broydena dla równania skalarnego to metoda siecznych. Własność minimalnej zmiany pochodnej w metodzie Broydena. Implementacja i koszt. Specjalne metody rozwiązywania układu z uaktualnioną macierzą pochodnej (formuła Sherman-Morrison-Woodbury).

W11, 09.05.00
Macierze rozrzedzone: przykłady i specyfika (pamięć, koszt mnożenia, itp) w porównaniu z macierzami gęstymi. Formaty przechowywania macierzy rzadkich (aij oraz R-B). Algorytmy iteracyjne dla macierzy rozrzedzonych, oparte na metodach przestrzeni Kryłowa. Uwaga o algorytmach bezpośrednich wykorzystujących reordering niewiadomych.

Metoda gradientów sprzężonych (CG) dla macierzy symetrycznej, dodatnio określonej, jako algorytm minimalizacji. Wielomian residualny. Metoda CG jako metoda bezpośrednia: twierdzenie o zbieżności w N krokach (w arytmetyce dokładnej). Wady CG jako metody bezpośredniej.

Ćw
Macierze symetryczne dodatnio określone (SPD) i ich pierwiastki - różne własności (m.in. norma wielomianu od macierzy SPD). Metody iteracyjne oparte na podziale macierzy; przykład: metoda Richardsona. Twierdzenie Vargi o zbieżności.

W12, 16.05.00
Metoda CG jako metoda iteracyjna i twierdzenie o oszacowaniu błędu. Implementacja (częściowe wyprowadzenie) metody CG i jej koszt.

Metoda GMRES dla macierzy niesymetrycznych i/lub nieokreślonych. GMRES jako algorytm minimalizacji. Sprowadzenie zadania minimalizacyjnego do zadania najmniejszych kwadratów przy użyciu procedury ortogonalizacyjnej Arnoldiego.

Lab
Porównanie metody stycznych i siecznych dla równania $\cos x - x = 0$ . Naprawdę duży układ równań nieliniowych: Implementacja w Octave metody Newtona dla dyskretyzacji w N punktach równania określającego H-funkcję Chandrasekhar'a (równania transferu radiacyjnego).

W13, 23.05.00
Metoda GMRES jako metoda bezpośrednia: twierdzenie o zbieżności w N krokach (w arytmetyce dokładnej). Twierdzenie o zbieżności metody GMRES jako metody iteracyjnej (bez dowodu). Uwagi dotyczące metody GMRES (reortogonalizacja, restart, uwarunkowanie).

Poprawianie uwarunkowania macierzy (ang. preconditioning). Kryteria, jakie powinien spełniać dobry preconditioner. Przykłady prostych preconditionerów (Jacobiego, Richardsona, ILU(0)). Preconditionery oparte na specyficznych własnościach macierzy - idea, na przykładzie metody dekompozycji obszaru.

Ćw
Dobór optymalnego parametru metody Richardsona. Metoda Jacobiego. Liczba kroków CG potrzebna do zredukowania błędu. Metody CGNR i CGNE. Równoważność spektralna i wykorzystanie w preconditioningu.
Praca domowa.

Nie zrobione:

Wielkie układy liniowe: algorytmy blokowe

Uwagi o programowaniu równoległym

Prace domowe 1999/2000

Prototyp metody Newtona
Materiał do samodzielnego opracowania dla grupy środowej i osób nieobecnych na zajęciach czwartkowych (w razie potrzeby, zapraszam na konsultacje).
1.
Zaprogramować w Matlabie funkcję newton, obliczającą rozwiązanie równania f(x)=0 dla zadanej funkcji f (przyjmujemy, że zarówno f jak jej pochodna są zadane jako funkcje w Matlabie, czyli, że są to tzw. M-funkcje).
Przykładowe rozwiązanie zadania: definiujemy trzy M-funkcje: funkcję newton - określającą właściwą metodę Newtona, moja_f - definiującą funkcję, której zera poszukujemy: u nas będzie to f(x) = x² - 1 oraz moja_df - jej pochodną, f'(x) = 2x. Należy pamiętać, aby nazwa pliku odpowiadała nazwie definiowanej w nim funkcji! Przy tak określonych funkcjach, mamy:
```
octave:15> [x,i] = newton('moja_f','moja_df',2,10,1e-15)
x = 1

i = 6

octave:16>
```
2.
Podać metodę Newtona przybliżonego znajdowania liczby 1/a, gdzie a>0, nie korzystającą z operacji dzielenia. (Jeśli na przykład nasz procesor nie ma implementowanej operacji dzielenia liczb podwyższonej precyzji (a ma mnożenie i dodawanie), to można używać tej metody do software'owej emulacji takiego dzielenia: jako dobre(!) przybliżenie startowe wybieralibyśmy odwrotność obliczoną w arytmetyce standardowej precyzji.) Znaleźć wyrażenie na błąd i warunek na przybliżenie początkowe, gwarantujący zbieżność metody.

3.
Niech A bedzie nieosobliwa macierzą kwadratową wymiaru n. Podać metodę iteracyjną analogiczną do poprzedniej, wyznaczającą przybliżenie A^-1 (tzn. gdy n=1, metoda redukuje się do poprzedniego zadania). Znaleźć wyrażenie na błąd E_k+1 przybliżenia X_k+1 macierzy A^-1,

$\begin{displaymath}E_{k} := I - A_kX_k, \qquad k = 0,1,\ldots, \end{displaymath}$

w terminach błędu na kroku poprzednim. Podać warunek na przybliżenie początkowe, gwarantujący zbieżność metody.
Eksport z Octave do Vigie

Na podstawie ćwiczeń laboratoryjnych, napisać M-funkcję w Octave, eksportującą 2-wymiarowe pole wektorowe zadane przez M-funkcję postaci
```
function z = VecField2D(x)
z(1) = ... funkcja zmiennych x(1) i x(2);
z(2) = ... funkcja zmiennych x(1) i x(2);
endfunction
```
Parametrami funkcji eksportującej mają być

1.
nazwa funkcji (u nas: VecField2D)
2.
zakres x₁ i x₂, czyli obszar, gdzie chcemy zwizualizować nasze pole (zakładamy, że jest to prostokąt).
3.
liczba węzłów siatki w każdym z kierunków
4.
nazwa bazowa pliku wynikowego (przypuśćmy, że będzie nią outvig).

Zatem jeśli zechcemy obejrzeć VecField2D w prostokącie [0,10]x[0,5] na siatce o 20 węzłach w każdym kierunku, to spodziewamy się wywoływać funkcję postaci:
```
octave:11> exportvigie_vf2d('VecField2D', 0, 10, 0, 5, 20, 20, 'outvig');
```
oraz
```
octave:12> help exportvigie_vf2d
```
W wyniku wywołania, nasza funkcja powinna wyświetlić w Octave obie składowe wizualizowanego pola na jednym wykresie oraz - to najważniejsze! - stworzyć dwa pliki: outvig.dsc oraz outvig.ascii2d. Przypominam strukturę tych plików:
outvig.dsc
```
	ascii2d
	outvig.ascii2d
```
outvig.ascii2d
```
	% komentarze
	points [liczba wezlow]
	.
	.
	[lista wspolrzednych wezlow (x,y)]
	.
	.
	quadrangles [liczba prostokatow siatki]
	.
	.
	[lista prostokatow, definiowanych jako sekwencja 4 indeksow wezlow zgodnie z
	ruchem wskazowek zegara (przypominam, ze wezly indeksujemy od zera!!)]
	.
	.
	vectors pole2d pole_wsp1 pole_wsp2 0.1
	.
	.
	[wartosci obu wspolrzednych pola w kolejnych wezlach]
	.
	.
```
Jak to oglądać w Vigie? (przypomnienie). Wywołujemy Vigie komendą
```
vigie -installcm
```
Następnie wybieramy Files -> Read datas i wczytujemy plik outvig.dsc. Dalej, w głównym panelu wybieramy 2D -> Muliti 2D plot; pojawi się wtedy panel roboczy 2D multiple solutions. Na tym panelu zaznaczamy Drawing areas, a następnie mesh oraz Ivalue. Odświeżamy obraz klikając redraw. Dalszy ciąg postępowania był przedstawiony na wykładzie; zob. też sekcję z linkami .

Dodatkowe opcje metody Newtona

zmodyfikować procedurę metody Newtona tak, by, oprócz rozwiązania, liczby iteracji i historii zbieżności, podawała także kod wykonania (0 - zbieżność, 1 - niezbieżność) oraz czas wykonania. Ponadto, w zależności od wartości dodatkowego parametru, będzie (lub nie) wyświetlany wykres historii zbieżności w skali półlogarytmicznej (semilogy()). Uwaga: pamiętajmy o odpowiedniej modyfikacji opisu funkcji wyświetlanego w helpie!
Wskazówka. Do liczenia czasu należy użyć funkcji tic i toc. Polecenie tic ``włącza stoper'', natomiast toc - wyłącza go, podając czas, który upłynął. Tak więc np. fragment kodu
```
tic;
...start obliczen
.
...obliczenia...
.
...koniec obliczen
czas = toc;
```
spowoduje, że zmienna czas będzie zawierać czas (w sekundach) wykonania bloku [start obliczeń] - [koniec obliczeń]).

Wariant metody cięciw

Przy jakich założeniach na $A\approx F'(x^*)$ , metoda (wariant metody cięciw)

x_n+1 = x_n - A^-1F(x_n)

jest lokalnie zbieżna liniowo do x^*? Podać oszacowanie szybkości zbieżności: ||e_n+1|| w zależności od ||e_n|| i ||A-F'(x^*)||.

Metoda siecznych

1.
Znaleźć x₀>0 takie, że metoda Newtona dla f(x) = arctg(x) będzie oscylować cyklicznie między x₀ a -x₀.

2.
Zaimplementować metodę siecznych, przez modyfikację ostatniej wersji funkcji realizującej metodę Newtona.

Metoda gradientów sprzężonych

Niech x_k będą iteracjami metody CG. Wykazać, że residua r_k = b - Ax_k spełniają:

$\begin{displaymath}r_l \in \mathcal{K}_k \quad \forall l < k, \end{displaymath}$

gdzie - k-ta przestrzeń Kryłowa.

Metoda Jacobiego i Richardsona

1.
Udowodnić, że jeśli macierz A jest silnie diagonalnie dominująca,

$\begin{displaymath}\vert a_{ii}\vert> \sum_{j\neq i}\vert a_{ij}\vert> 0, \end{displaymath}$

to metoda Jacobiego jest zbieżna w normie $\vert\vert\cdot\vert\vert _\infty$ . Ponadto, macierz A jest wtedy nieosobliwa.

2.
Dla jakich parametrów $\tau$ , metoda Richardsona jest zbieżna. Znaleźć optymalną wartość parametru $\tau$ .