Na tych laboratoriach skupimy się na rozwiązaniu układu równań: \[
\mathbf{A} \mathbf{x} = \mathbf{b}
\] W tym celu rozpatrzymy kilka metod, których działanie sprawdzimy na układzie równań otrzymanym na poprzednich zajęciach dzięki metodzie MES.
Naszym celem będzie napisanie funkcji Solve, która zastąpi funkcję Gauss i wyznaczy wartość wektora \(\mathbf{x}\). Nie będziemy jednak tego układu rozwiązywać metodą bezpośrednią, taką jak eliminacja Gaussa, ale metodą iteracyjną. Skonstruujemy ciąg \(\left(\mathbf{x}^{k}\right)\), którego elementy \(\mathbf{x}^{k}\) będą dążyły do rozwiązania dokładnego \(\mathbf{x}\).
Niezbędne definicje
Łatwo zauważyć, że zachodzi zależność \(\mathbf{r}^k = \mathbf{A} \mathbf{e}^k\).
Widać także, że skoro \(\mathbf{x}^{k}\) dąży do \(\mathbf{x}\) to \(\mathbf{r}^{k}\) dąży do zera.
Wyznacz residuum dla zadania z poprzednich zajęć. Następnie oblicz i wyświetl jego normę: \(\|\mathbf{r}\| = \sqrt{\mathbf{r}^T\mathbf{r}}\) (napisz funkcję liczącą normę wektora norm(int, double *)). Ile wynosi ta norma przed i po rozwiązaniu układu metodą eliminacji Gaussa?
Weźmy dowolny wektor \(\mathbf{x}^0\) i obliczmy odpowiadający mu wektor prawych stron \(\mathbf{b}^0 = \mathbf{A} \mathbf{x}^0\). Różnica między ,,prawdziwym’’ wektorem \(\mathbf{b}\) a przybliżeniem jest wtedy równa \[ \mathbf{r}^0 = \mathbf{b} -\mathbf{b}^0 = \mathbf{A} \mathbf{x} - \mathbf{A} \mathbf{x}^0 = \mathbf{A} \mathbf{e}^0 \] Zatem różnica między ,,prawdziwym’’ rozwiązaniem a przybliżonym \(\mathbf{e}^0 = \mathbf{A}^{-1} \mathbf{r}^0\). Co ostatecznie pozwala nam zapisać: \(\mathbf{x} = \mathbf{x}^0 + \mathbf{A}^{-1}\mathbf{r}^0\). Nie mamy jednak \(\mathbf{A}^{-1}\) (w tym rzecz). Zamiast niej użyjemy macierzy \(\mathbf{M}^{-1}\). Wtedy jednak nie dostaniemy dokładnej wartości1 \(\mathbf{x}\) a jedynie przybliżenie. Prowadzi to nas do wzoru: \[ \mathbf{x}^{k+1} = \mathbf{x}^k + \mathbf{p}^k = \mathbf{x}^k + \mathbf{M}^{-1}\mathbf{r}^k \] Wektor \(\mathbf{p}^k\) jest ,,poprawką’’ w k-tej iteracji a macierz \(\mathbf{M}\) nazywamy preconditioner’em. Najlepiej byłoby gdyby macierz \(\mathbf{M}\) była ,,podobna’’ do macierzy \(\mathbf{A}\) a jednocześnie łatwo odwracalna.
Rozpatrzmy układ \(\mathbf{r}^k = \mathbf{A} \mathbf{p}^k\). Widać, że jeśli pominiemy większość jego elementów: \[ \left\{ \begin{array}{ccccccccccc} A_{11} p_1^k &+& \color{Gray}{A_{12} p_2^k} &+& \color{Gray}{A_{13} p_3^k} &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{1N} p_N^k} &=& r_1^k \\ \color{Gray}{A_{21} p_1^k} &+& A_{22} p_2^k &+& \color{Gray}{A_{23} p_3^k} &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{2N} p_N^k} &=& r_2^k \\ \color{Gray}{A_{31} p_1^k} &+& \color{Gray}{A_{32} p_2^k} &+& A_{33} p_3^k &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{3N} p_N^k} &=& r_3^k \\ \vdots & & \vdots & & \vdots & & \ddots & & \vdots & & \vdots \\ \color{Gray}{A_{N1} p_1^k} &+& \color{Gray}{A_{N2} p_2^k} &+& \color{Gray}{A_{N3} p_3^k} &+& \color{Gray}{\ldots} &+& A_{NN} p_N^k &=& r_N^k \\ \end{array} \right. \] dostaniemy prosty wzór na \(\mathbf{p}^k\): \[ p^k_i = \frac{1}{A_{ii}}{r^k_i} \]
Jest to równoważne z wzięciem za macierz \(\mathbf{M}\) diagonalnej części macierzy \(\mathbf{A}\). Ten prosty schemat iteracji z powyższą poprawką nazywamy metodą Jacobiego.
res_draw(double), która posłuży do wykonania wykresu zbieżności.Solve część odpowiedzialną za mnożenie przez \(\mathbf{A}\): Mult(int N, double **A, double *x, double *r) i preconditioner: Precond(int N, double **A, double *x, double *p).Spróbujmy poprawić nasz schemat biorąc lepszy preconditioner. Zauważmy, że obliczając \(p_2^k\) mamy już obliczone \(p_1^k\) i możemy go użyć. Nie musimy zatem pomijać elementów układu ,,pod diagonalą’’: \[ \left\{ \begin{array}{ccccccccccc} A_{11} p_1^k &+& \color{Gray}{A_{12} p_2^k} &+& \color{Gray}{A_{13} p_3^k} &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{1N} p_N^k} &=& r_1^k \\ A_{21} p_1^k &+& A_{22} p_2^k &+& \color{Gray}{A_{23} p_3^k} &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{2N} p_N^k} &=& r_2^k \\ A_{31} p_1^k &+& A_{32} p_2^k &+& A_{33} p_3^k &+& \color{Gray}{\ldots} &+& \color{Gray}{A_{3N} p_N^k} &=& r_3^k \\ \vdots & & \vdots & & \vdots & & \ddots & & \vdots & & \vdots \\ A_{N1} p_1^k &+& A_{N2} p_2^k &+& A_{N3} p_3^k &+& \ldots &+& A_{NN} p_N^k &=& r_N^k \\ \end{array} \right. \] Daje nam to prosty wzór na \(\mathbf{p}^k\): \[ p_i^k = \frac{1}{A_{ii}}\left( r_i - \sum_{j=1}^{i-1} A_{ij} p_j^k \right) \] Gdy \(\alpha = 1\) schemat taki nazywamy metodą Gaussa-Seidla.
Schematy, dla których \(\alpha > 1\) nazywamy metodami Successive Over-Relaxation (SOR).
Widać wyraźnie, że zbieżność bardzo zależy od wartości współczynnika \(\alpha\) i jasnym jest, że najlepiej byłoby dobierać ten współczynnik w każdej iteracji: \[ \mathbf{x}^{k+1} = \mathbf{x}^{k} + \alpha^k \mathbf{p}^k \] Zastanówmy się teraz jak będzie się zmieniało residuum w zależności od kroku. Jeśli pomnożymy powyższy wzór przez \(-\mathbf{A}\) a następnie dodamy \(\mathbf{b}\) i skorzystamy z definicji residuum otrzymamy: \[ \mathbf{r}^{k+1} = \mathbf{r}^{k} - \alpha^k \mathbf{A} \mathbf{p}^k \] Kwadrat normy tego residuum jest równy: \[ \|\mathbf{r}^{k+1}\| = (\mathbf{r}^{k+1})^T \mathbf{r}^{k+1} = (\mathbf{r}^{k} - \alpha^k \mathbf{A} \mathbf{p}^k)^T (\mathbf{r}^{k} - \alpha^k \mathbf{A} \mathbf{p}^k) = \] \[ (\mathbf{r}^{k})^T \mathbf{r}^{k} - 2 \alpha^k (\mathbf{r}^{k})^T \mathbf{A} \mathbf{p}^k + (\alpha^k)^2(\mathbf{A}\mathbf{p}^k)^T \mathbf{A}\mathbf{p}^k \] Widać, że kwadrat normy jest kwadratową funkcją \(\alpha^k\) a współczynnik przed \((\alpha^k)^2\) jest dodatni. Oznacza to, że funkcja ta ma minimum. Obliczamy pochodną po \(\alpha^k\): \[ \frac{d}{d\alpha^k} \left( \|\mathbf{r}^{k+1}\| \right) = -2(\mathbf{r}^{k})^T \mathbf{A} \mathbf{p}^k + 2 \alpha^k (\mathbf{A} \mathbf{p}^k)^T \mathbf{A} \mathbf{p}^k \] i przyrównujemy do zera co ostatecznie daje wartość: \[ \alpha^k = \frac{(\mathbf{r}^{k})^T \mathbf{A} \mathbf{p}^k}{(\mathbf{A} \mathbf{p}^k)^T \mathbf{A} \mathbf{p}^k} \] Schemat z taką wartością \(\alpha^k\) nazywamy metodą najmniejszych residuów (Minimal Residual Method — MINRES).
skal(int, double *, double *) liczącą iloczyn skalarny i oblicz \(\alpha^k\) z powyższego wzoru.Do tej pory ignorowaliśmy informację o poprawkach z poprzednich iteracji i poprawkę w k-tym kroku obliczaliśmy ze wzoru \[ \mathbf{p}^k = \mathbf{M}^{-1}\mathbf{r}^k \] Zmienimy to przez wykorzystanie informacji o poprawce z \(k-1\) kroku: \[ \mathbf{p}^k = \mathbf{p}^k - \beta^k \mathbf{p}^{k-1} \] Teraz wzór na nowe residuum będzie miał postać: \[ \mathbf{r}^{k+1} = \mathbf{r}^{k} - \alpha^k \mathbf{A} \left(\mathbf{p}^k - \beta^k \mathbf{p}^{k-1}\right) \] Musimy jeszcze wyznaczyć wartość nowego współczynnika \(\beta^k\).
Schemat, w którym współczynniki \(\alpha^k\) i \(\beta^k\) obliczane są po przez minimalizację residuów nazywamy uogólnioną metodą najmniejszych residuów (Generalized Minimal Residual Method — GMRES).
W przypadku naszego zadania MES możemy wykorzystać fakt, że macierz \(\mathbf{A}\) jest symetryczna i dodatnio określona. Wtedy zamiast minimalizować \((\mathbf{r}^{k+1})^T \mathbf{r}^{k+1}\) możemy zminimalizować pewien specjalny funkcjonał: \[ f(\mathbf{x}) = \frac{1}{2} \mathbf{x}^T \mathbf{A} \mathbf{x} - \mathbf{b}^T \mathbf{x} \]
Zastosuj identyczny schemat iteracji jak w poprzednim punkcie ale zmień \(\alpha^k\) i \(\beta^k\). Zbadaj zbieżność.
Taki schemat nazywamy metodą gradientu sprzężonego (Conjugate Gradient Method — CG).
Uwaga: Aktualnie zbieżność jest bardzo słaba. Wynika to z faktu, że choć \(\mathbf{A}\) jest symetryczna to preconditioner z metody Gaussa-Seidla \(\mathbf{M}^{-1}\) już nie jest.
Zbadaj zbieżność z preconditionerem diagonalnym, lub wyrażeniem \(\mathbf{p}^k = \mathbf{r}^k\) (brakiem preconditionera).
Uwaga: Metodę Conjugate Gradient można zaimplementować w bardziej ,,zwartej’’ formie. Taki schemat można znaleźć na Wikipedii, bądź w notatkach z wykładu.
I tak nie dostalibyśmy dokładnej wartości ze względu na błędy numeryczne.↩︎