Zagadnienia

8.1 Inne metody Kryłowa
- 8.1.1 Uogólnienia metody CG
- 8.1.2 Metody nie oparte na minimalizacji
8.2 Metody strukturalne
8.3 Metody wielosiatkowe
- Pełna metoda wielosiatkowa

8. Przegląd innych metod rozwiązywania wielkich układów równań liniowych

Na szczęście, nie ma jednej najlepszej metody rozwiązywania wielkich układów równań liniowych; czyni to oczywiście życie znacznie ciekawszym. Doprowadziło to do skonstruowania dużej liczby konkurencyjnych metod: zarówno bezpośrednich, jak i iteracyjnych. Celem nieniejszego rozdziału jest pobieżne zaznajomienie Czytelnika z niewielką ich reprezentacją. Dobór tej właściwej dla konkretnego zadania nie musi być oczywisty! Szerszy przegląd takich metod znajdziemy m.in. w monografii [13], [].

8.1. Inne metody Kryłowa

Na początek, na podstawie [1] przedstawimy w pewien systematyczny sposób uogólnienie metody CG na szerszą klasę problemów. Następnie zrobimy krótki przegląd wybranych metod Kryłowa, nie opartych na zasadzie minimalizacji.

8.1.1. Uogólnienia metody CG

Wielką zaletą metody PCG jest to, że nie wymaga ona pamiętania całej bazy przestrzeni Kryłowa $K_{k}(Mr_{0},MA)=\text{span}\{ Mr_{0},MAMr_{0},\ldots,(MA)^{{k-1}}Mr_{0}\}$ , a iteracja ma prostą formułę $x_{{k+1}}=x_{k}+d_{k}$ , gdzie $d_{k}\in K_{k}(Mr_{0},MA)$ jest tak dobrany by minimalizować pewną normę błędu. Można więc pójść krok dalej i określić całą klasę metod o podobnej strukturze.

Przyjmijmy więc, że

macierz $A$ jest nieosobliwa (na tym poziomie ogólności nie zakładamy symetrii ani dodatniej określoności)

oraz są dane dwie dodatkowe macierze:

nieosobliwa macierz lewostronnie ściskająca $M$ ,
symetryczna, dodatnio określona macierz $B$ , która będzie definiować normę, w której będziemy minimalizować błąd:

$\| x\| _{B}=(x^{T}Bx)^{{1/2}}.$

Będziemy zakładali, że macierz ściśnięta $MA$ jest normalna względem iloczynu skalarnego indukowanego przez macierz $B$ , co możemy zapisać w formie warunku

$GG^{T}=G^{T}G,$

gdzie $G=B^{{1/2}}MAB^{{-1/2}}$ .

W uogólnionej metodzie CG, którą będziemy oznaczali GCG( $B$ , $M$ , $A$ ), kolejna iteracja $x_{{k+1}}\in x_{0}+K_{k}(Mr_{0},MA)$ będzie określona tak, by

$\displaystyle\| x_{{k+1}}-x^{*}\| _{B}\leq\| x-x^{*}\| _{B}\qquad\forall x\in x_{0}+K_{k}(Mr_{0},MA)$

(8.1)

Twierdzenie 8.1

Przy powyższych założeniach, metoda GCG( $B$ , $M$ , $A$ ) jest dobrze określona i w dokładnej arytmetyce osiąga rozwiązanie dokładne po co najwyżej $N$ krokach.

Dowód

Kładąc $\tilde{r}_{0}=Mr_{0}$ oraz $\tilde{A}=MA$ , możemy zastosować twierdzenie 6.1 do przestrzeni Kryłowa $K_{k}=K_{k}(\tilde{r}_{0},\tilde{A})$ , które zagwarantuje nam spełnienie tezy twierdzenia.

∎

Twierdzenie 8.2

Przy powyższych założeniach, błąd na $k$ -tym kroku metody GCG( $B$ , $M$ , $A$ ) spełnia oszacowanie

$\| x_{k}-x^{*}\| _{B}\leq\min _{{p\in\tilde{P}_{k}}}\max _{{\lambda\in\sigma(MA)}}|p(\lambda)|\| x_{0}-x^{*}\| _{B}.$

Dowód

Na mocy wniosku 6.1,

$\| x_{k}-x^{*}\| _{B}=\min _{{p\in\tilde{P}_{k}}}\| p(MA)(x_{0}-x^{*})\| _{B}.$

Ponieważ dla dowolnego $y$

$\displaystyle\|(MA)^{k}y\| _{B}=\|(B^{{1/2}}MAB^{{-1/2}})^{k}\, B^{{1/2}}y\| _{2},$

(8.2)

to w konsekwencji $\| p(MA)y\| _{B}=\| p(B^{{1/2}}MAB^{{-1/2}})\, B^{{1/2}}y\| _{2}\leq\| p(B^{{1/2}}MAB^{{-1/2}})\| _{2}\,\| y\| _{B}$ . Z założenia macierz $G=B^{{1/2}}MAB^{{-1/2}}$ jest normalna, jest więc diagonalizowalna i jej wektory własne są ortogonalne:

$G=X\Lambda X^{T},\qquad\text{ gdzie }X^{T}X=I.$

Stąd zaś wynika, że $\| p(G)\| _{2}=\| Xp(\Lambda)X^{T}\| _{2}=\| p(\Lambda)\| _{2}=\max _{{\lambda\in\sigma(G)}}|p(\lambda)|$ , co kończy dowód, gdyż macierze $G$ i $MA$ są podobne, więc mają to samo spektrum.

∎

Ćwiczenie 8.1

Udowodnij równość (8.2).

Można, postępując analogicznie jak w przypadku klasycznej metody CG wykazać, że GCG( $B$ , $M$ , $A$ ) ma własności podobne jak jej protoplastka:

Stwierdzenie 8.1

W metodzie GCG( $B$ , $M$ , $A$ ) zachodzi:

$x_{k}=x_{{k-1}}+\alpha _{k}p_{k}$ , gdzie $p_{k}$ tworzą bazę $B$ -ortogonalną przestrzeni Kryłowa $K_{k}(Mr_{0},MA)$ .
$w^{T}B(x_{k}-x^{*})=0$ dla każdego $w\in K_{k}(Mr_{0},MA)$ .

Dowód

W oczywisty sposób $x_{k}=x_{{k-1}}+d_{k}$ , gdzie $d_{k}\in K_{k}(Mr_{0},MA)$ . Jeśli przez $V_{k}$ oznaczyć bazę w $K_{k}(Mr_{0},MA)$ , to $x_{k}=x_{0}+V_{k}a_{k}$ i na mocy (6.2) $a_{k}$ spełnia układ równań normalnych,

$V_{k}^{T}B(V_{k}a_{k}-(x^{*}-x_{0}))=0.$

Ponieważ $V_{k}a_{k}=x_{k}-x_{0}$ , dostajemy

$V_{k}^{T}B(x_{k}-x^{*})=0,$

co dowodzi drugiego punktu.

∎

Metodę GCG( $B$ , $M$ , $A$ ) możnaby zrealizować następującym uogólnieniem algorytmu CG, który nazwiemy za [1] algorytmem Odir( $B$ , $M$ , $A$ ):

Metoda Odir dla GCG, wersja bazowa

$p_{0}=Mr_{0}$ ; $k=0$ ;

while not stop begin

$\alpha _{k}=\dfrac{p_{k}^{T}B(x^{*}-x)}{p_{k}^{T}Bp_{k}}$

$x_{k}=x_{{k-1}}+\alpha _{k}p_{k}$

$r_{k}=r_{{k-1}}-\alpha _{k}Ap_{k}$

$\gamma _{k}=\dfrac{p_{k}^{T}BMAp_{k}}{p_{k}^{T}Bp_{k}}$

$\sigma _{k}=\dfrac{p_{k}^{T}BMAp_{{k-1}}}{p_{{k-1}}^{T}Bp_{{k-1}}}$

$p_{{k+1}}=MAp_{k}-\gamma _{k}p_{k}-\sigma _{k}p_{{k-1}}$

$k=k+1$

end

Algorytm ten jeszcze nie jest gotowy do użycia, albowiem wymaga obliczania $B(x^{*}-x)$ — a więc potencjalnie wymaga odwołania się do nieznanego a priori wektora błędu! Możliwość skutecznego obliczania współczynnika $\alpha _{k}$ ogranicza więc zbiór $B$ , których możemy użyć do zdefiniowania konkretnej metody. Z drugiej strony, wybierając konkretne $B$ możemy dalej uprościć powyższą bazową implementację. Na tym etapie nie jest także oczywiste, czy może zdarzyć się $p_{k}=0$ — a więc stagnacja i załamanie się algorytmu (dzielenie przez zero!).

Wśród metod, które dają się skutecznie zaimplementować na podstawie algorytmu Odir( $B$ , $M$ , $A$ ) — to znaczy: mają obliczalne $\alpha _{k}$ — znajdują się metody wymienione w tabeli 8.1, którą przytaczamy za [1, Table 5.2]. Z tabeli wynika m.in., że PCG wcale nie wymaga dodatnio określonej macierzy ściskającej!

Metoda	$B$	$M$	Ograniczenia
CG	$A$	$I$	$A=A^{T}>0$
CR	$A^{2}$	$I$	$A=A^{T}$
PCG	$A$	$M$	$A=A^{T}>0$ , $M=M^{T}$
PCR	$AMA$	$M$	$A=A^{T}$ , $M=M^{T}>0$
CGNR	$A^{T}A$	$A^{T}$
CGNE	$I$	$A^{T}$
D'yakonov	$B$	$B^{{-1}}A^{T}B^{{-1}}$	$B=B^{T}>0$

Tabela 8.1. Wybrane metody Kryłowa dające się zinterpretować jako metoda GCG. W tabeli zamieszczono założenia na macierze $A$ i $M$ gwarantujące poprawność określenia metody.

Aby otrzymać dobrą implementację konkretnej metody opartej na algorytmie Odir, należy zawsze ją dopracować, wykorzystując zależności specyficzne dla konkretnej metody.

Warta uwagi jest metoda PCR, która działa w przypadku, gdy macierz $A$ jest jedynie symetryczna — nie musi być dodatnio określona. Metoda D'aykonova [6] co prawda działa dla dowolnej nieosobliwej, niesymetrycznej macierzy $A$ , ale ze względu na jej podobieństwo do równań normalnych dla $B^{{-1}}A$ , w wielu wypadkach skuteczniejsza (pod względem szybkości zbieżności, a nie pod względem oszczędności pamięci) będzie metoda typu GMRES.

Uwaga 8.1

Czasami wymaganie, by można było wykonywać mnożenie przez $A^{T}$ może być trudne do spełnienia — na przykład wtedy, gdy $A$ jest zadana jako operator, tzn. wyłącznie przez procedurę mnożenia przez zadany wektor $x$ , tzn. obliczania $A\cdot x$ .

Ćwiczenie 8.2

Porównaj zbieżność metod: PCG i PCR w przypadku, gdy $A$ oraz $M$ są macierzami symetrycznymi i dodatnio określonymi. Przeprowadź weryfikację eksperymentalną swoich przypuszczeń.

Wskazówka:

W MATLABie dyponujesz gotową implementacją zarówno metody pcg, jak i pcr.

8.1.2. Metody nie oparte na minimalizacji

Zamiast określać kolejną iterację metody, jak dotychczas, przez pewien warunek minimalizacji, możemy inaczej położyć warunek na $x_{k}$ . Na przykład, możemy wymagać by oprócz $x_{k}\in x_{0}+K_{k}(r_{0},A)$ zachodził warunek ortogonalności residuów:

$v^{T}r_{k}=0\qquad\forall v\in K_{k}(r_{0},A^{T}).$

(Zwróć uwagę na to, że residua mają być ortogonalne nie do $K_{k}(r_{0},A)$ , tylko do $K_{k}(r_{0},A^{T})$ .)

Tego typu warunek prowadzi m.in. do metody BiCG (ang. Bi-Conjugate Gradient) która wymaga mnożenia przez $A^{T}$ i nie jest też zbyt stabilna. Jej modyfikacje: CGS i BiCG-stab — nie wymagają już mnożenia przez $A^{T}$ , BiCG-stab jest też bardziej stabilna. Niestety, jak na razie dla BiCG-stab nie ma pełnej i satysfakcjonującej teorii zbieżności.

W przeciwieństwie do metod minimalizacji, kolejna iteracja może nie być realizowalna mimo, że nie osiągnięto jeszcze rozwiązania — mówimy wtedy potocznie o załamaniu się metody (ang. breakdown).

W innej metodzie, QMR (ang. Quasi-Minimal Residual) — i jej wariancie bez mnożenia przez macierz transponowaną: TFQMR (ang. Transpose-Free QMR) — kolejną iterację wybiera się tak, by zminimalizować (stosunkowo łatwą do obliczenia) wielkość, która ma tylko trochę wspólnego z normą residuum.

Można pokazać [9], że metoda TFQMR dla macierzy $N\times N$ (realizowana w dokładnej arytmetyce) w ciągu $\lceil(N+1)/2\rceil$ iteracji albo załamie się, albo osiągnie dokładne rozwiązanie.

8.2. Metody strukturalne

W bardzo wielu przypadkach, w macierzach układów równań jakie przychodzi nam rozwiązywać w praktycznych zastosowaniach, nie tylko niewiadome są powiązane jedynie z nielicznymi innymi niewiadomymi (co właśnie prowadzi do rozrzedzenia macierzy) ale w dość naturalny sposób można wskazać grupy niewiadomych i równań, które między sobą nie mają żadnych powiązań.

W najogólniejszym przypadku, mielibyśmy więc

$A=\begin{pmatrix}A_{{11}}&&A_{{13}}\\ &A_{{22}}&A_{{23}}\\ A_{{31}}&A_{{32}}&A_{{33}},\end{pmatrix}$

(8.3)

gdzie $A_{{ii}}$ są pewnymi macierzami kwadratowymi. Niewiadome i równania z grupy ,, $1$ ” nie mają bezpośredniego związku z grupą ,, $2$ ” (są one powiązane ze sobą jedynie poprzez równania i niewiadome grupy ,, $3$ ”). W dalszym ciągu będziemy zakładać, że macierze $A_{{11}}$ i $A_{{22}}$ są nieosobliwe. Strukturalizacja równania liniowego to nic innego jak wskazanie pewnej permutacji równań i niewiadomych, w wyniku której dostajemy układ postaci (8.3). W rzeczywistości, dla macierzy $A_{{11}}$ i $A_{{22}}$ mogłyby zachodzić podobne relacje¹⁰Taką permutację nazywamy wtedy zagnieżdżonym podziałem grafu macierzy (ang. nested dissection)., jak w poniższym przykładzie.

Przykład 8.1 (Strukturalizacja dyskretyzacji równania różniczkowego)

Rozważmy dyskretyzację równania różniczkowego

$-\Delta u=f\text{ w }\Omega$

z jednorodnymi warunkami brzegowymi Dirichleta. Jeśli obszar $\Omega$ podzielić na kilka mniejszych części, podobszarów, $\bar{\Omega}=\cup _{{i=1}}^{K}\bar{\Omega}_{i}$ i przez $\Gamma$ oznaczyć zestaw krawędzi (ang. interface) je spajających: $\Gamma=\left(\cup _{{i=1}}^{K}\partial\Omega _{i}\right)-\partial\Omega$ , to — ze względu na lokalny charakter działania pochodnej (i jej sensownych aproksymacji) — niewiadome z dwóch różnych podobszarów $\Omega _{i}$ i $\Omega _{j}$ nie będą ze sobą powiązane. Będą za to powiązane z niewiadomymi leżącymi na krawędziach spajających $\Gamma$ , por. rysunek 8.1.

Obszar prostokątny jest podzielony na krzyż na cztery mniejsze prostokąty. Brzegi mniejszych prostokątów (z wyłączeniem brzegu wyjściowego obszaru) tworzą interfejs --- w kształcie krzyża.

Rys. 8.1. Przykładowy podział obszaru na $K=4$ podobszary; węzły wspólne (należące do $\Gamma$ ) zaznaczono kolorem niebieskim.

Numerując niewiadome tak, by $U=(U_{1},U_{2},\ldots,U_{K},U_{\Gamma})^{T}$ , gdzie wektor $U_{i}$ zawiera niewiadome odpowiadające wewnętrznym węzłom w $\Omega _{i}$ , a $U_{\Gamma}$ — niewiadome odpowiadające węzłom na interfejsie $\Gamma$ , otrzymalibyśmy macierz o blokowej strukturze

$D_{h}=\begin{pmatrix}D_{{11}}&&&&D_{{1\Gamma}}\\ &D_{{22}}&&&D_{{2\Gamma}}\\ &&\ddots&&\vdots\\ &&&D_{{KK}}&D_{{K\Gamma}}\\ D_{{\Gamma 1}}&D_{{\Gamma 2}}&\cdots&D_{{\Gamma K}}&D_{{\Gamma\Gamma}}\end{pmatrix}_{{N\times N}}.$

(8.4)

W pustych miejscach $D_{h}$ znajdują się bloki zerowe, a macierze $D_{{ii}}$ faktycznie są nieosobliwe.

Dzięki szczególnej postaci macierzy (8.3) możemy podać algorytm typu ,,dziel i rządź” rozwiązywania układu równań $Ax=b$ z tą macierzą. Wprowadzając podział niewiadomych i wektora prawej strony zgodnie z blokowym podziałem macierzy $A$ zadanym (8.3) mamy do rozwiązania układ

$\begin{pmatrix}A_{{11}}&&A_{{13}}\\ &A_{{22}}&A_{{23}}\\ A_{{31}}&A_{{32}}&A_{{33}}\end{pmatrix}\begin{pmatrix}x_{{1}}\\ x_{{2}}\\ x_{{3}}\end{pmatrix}=\begin{pmatrix}b_{{1}}\\ b_{{2}}\\ b_{{3}}\end{pmatrix}.$

Dokonując blokowej eliminacji Gaussa (najpierw na blokach $A_{{11}}$ i $A_{{22}}$ ) dostajemy następujący algorytm:

Oblicz $q_{i}=A_{{ii}}^{{-1}}b_{i}$ , $i=1,2$ .
Wyznacz $x_{3}$ jako rozwiązanie układu

$Sx_{3}=f,$

gdzie $S=A_{{33}}-A_{{31}}A_{{11}}^{{-1}}A_{{13}}-A_{{32}}A_{{22}}^{{-1}}A_{{23}}$ oraz $f=b_{3}-A_{{31}}q_{1}-A_{{32}}q_{2}$ .
Wyznacz $x_{i}=A_{{ii}}^{{-1}}(b_{i}-A_{{i3}}x_{3})$ , $i=1,2$ .

Oczywiście, wszędzie tam, gdzie występują odwrotności macierzy, w praktyce będziemy rozwiązywali układy równań z tymi macierzami.

Zwróćmy uwagę na kilka charakterystycznych cech tego algorytmu:

Zarówno w kroku 1. jak i w kroku 3. algorytmu, oba układy możemy rozwiązywać niezależnie od siebie (co daje możliwość równoległej implementacji).
Macierze $A_{{ij}}$ są rozrzedzone (bo $A$ taka była), ale w ogólności $S$ jest macierzą gęstszą.

Jeśli wymiar macierzy $S$ jest niezbyt wielki, możemy ją wyznaczyć explicite i rozwiązać metodą bezpośrednią. Jednak jeśli $S$ wciąż jest zbyt wielka (lub zbyt gęsta), by potraktować ją metodą bezpośrednią — wtedy możemy zastosować do niej metodę iteracyjną operatorową (matrix-free). Rzeczywiście, mnożenie macierzy $S$ przez wektor $z$ można zrealizować bez wyznaczania niej samej:

$Sz=A_{{33}}z-A_{{31}}(A_{{11}}^{{-1}}(A_{{13}}z))-A_{{32}}(A_{{22}}^{{-1}}(A_{{23}}z)).$

(Nawiasy wskazują kolejność mnożenia.) Możemy więc wyznaczyć iloczyn $Sz$ kosztem rozwiązania układów równań z macierzami $A_{{11}}$ i $A_{{22}}$ (co i tak musielibyśmy umieć zrobić!) oraz niewielkim dodatkowym kosztem pomnożenia pewnych wektorów przez rozrzedzone macierze trzeciej kolumny i trzeciego wiersza $A$ .

Ćwiczenie 8.3

Przedyskutuj, jak rozwiązywać układ równań postaci (8.4).

8.3. Metody wielosiatkowe

Macierze rozrzedzone, które otrzymuje się z dyskretyzacji równań różniczkowych cząstkowych metodą różnic skończonych lub elementu skończonego, mają wiele specjalnych cech, które pozwalają na konstrukcję optymalnych metod rozwiązywania takich zadań. Jedną z takich klas metod są metody wielosiatkowe (ang. multigrid).

Omówimy jej ideę na przykładzie macierzy $T_{N}$ dyskretyzacji jednowymiarowego laplasjanu (5.5): wszystkie istotne wady metod iteracyjnych ujawniają się już w jej przypadku, podobnie jak zasadnicze cechy metody wielosiatkowej. Oczywiście, w praktyce metodę wielosiatkową stosowalibyśmy dopiero do dyskretyzacji dwu- lub trójwymiarowego laplasjanu.

Gdyby przyjrzeć się metodzie Jacobiego dla $T_{N}$ zobaczymy, że niektóre składowe błędu są w niej szybciej redukowane niż inne. Dokładniej, rozważmy metodę Jacobiego z parametrem relaksacyjnym $\omega$ ,

$x_{{k+1}}=x_{k}+\dfrac{\omega}{2}r_{k}$

(w macierzy $T_{N}$ diagonala jest stała równa $2I$ ). Macierz iteracji tej metody ma postać

$G=I-\dfrac{\omega}{2}T_{N}.$

Wyrażając błąd $e_{k}=x_{k}-x^{*}$ w bazie wektorów własnych macierzy $G$ (są dane wzorem $(v_{i})_{j}=\sin(\dfrac{ij\pi}{N+1})$ — por. (5.18)) dostajemy, że

$e_{0}=\sum _{i}a_{i}v_{i}$

i w konsekwencji,

$e_{k}=G^{k}e_{0}=\sum _{i}\lambda _{i}(G)^{k}a_{i}v_{i}.$

Jasne jest więc, że najszybszej redukcji ulegną te składowe błędu, które odpowiadają najbliższym zera wartościom $\lambda _{i}(G)$ , najwolniej zaś będą znikać te, dla których $|\lambda _{i}(G)|$ będzie duży. Rysunek 8.2 przedstawia wykres

$\lambda _{i}(G)=1-\omega(1-\cos(\dfrac{i\pi}{N+1}))$

w zależności od parametru relaksacyjnego $\omega$ .

Gdy omega jest wyraźnie mniejsza od 1, np. równa około 1/2, następuje szybka redukcja szybkozmiennych składowych błędu

Rys. 8.2. Wartości własne macierzy iteracji metody Jacobiego z parametrem $\omega\in\{ 1,\frac{3}{4},\frac{1}{2},\frac{1}{4}\}$ , zastosowanej do macierzy $T_{{100}}$ .

Jak możemy zauważyć, gdy $\omega=1$ (tzn. gdy używamy standardowej metody Jacobiego), najsłabiej są redukowane składowe szybkozmienne i wolnozmienne (jedynie te ,,średniozmienne” są redukowane błyskawicznie).

Ale już dla $\omega=1/2$ , najmocniej są redukowane wszystkie szybkozmienne składowe błędu — zatem po kilku iteracjach takiej metody błąd zostanie wygładzony. Jeśliby więc wspomóc taką prostą iterację poprawką, która redukowałaby pozostałe — wolnozmienne — składowe błędu, moglibyśmy otrzymać szybko zbieżną iterację. Jak pamiętamy z wcześniejszego wykładu, metody projekcji zastępowały idealną poprawkę $A\delta^{*}=r_{k}$ poprawką przybliżoną $\delta$ , wyznaczaną przez rozwiązanie zadania zmniejszonego wymiaru. W metodzie dwusiatkowej poprawkę wyznaczymy podobnie: rozwiązując metodą bezpośrednią (np. eliminacji Gaussa) równanie poprawki na rzadszej siatce — wszak będzie ono tylko gorszą aproksymacją (o mniejszej rozdzielczości!) tego samego rozwiązania równania Poissona, które ma aproksymować $x^{*}$ .

W metodzie wielosiatkowej zadanie na rzadszej siatce (potrzebne do wyznaczenia poprawki na siatce najdrobniejszej) potraktujemy rekurencyjnie w ten sam sposób, otrzymując sekwencję zadań coraz mniejszego rozmiaru, którą w pewnym momencie oczywiście będziemy musieli przerwać i rozwiązać zadanie bardzo małego wyniaru metodą bezpośrednią.

Oto więc kolejne etapy jednej iteracji metody wielosiatkowej, wyznaczającej kolejne przybliżenie $x_{{k+1}}$ na podstawie $x_{k}$ :

Wygładzenie (ang. pre-smoothing) — zaczynając od $x_{k}$ , wykonujemy kilka iteracji (zwykle jedną lub dwie) prostej metody typu Jacobiego lub Gaussa–Seidela z właściwie dobranym parametrem $\omega$ . Dostajemy przybliżenie pośrednie $x_{{k+1}}^{{(1)}}$ , ze zredukowaną szybkozmienną składową błędu.
Restrykcja — sformułowanie na rzadszej siatce zadania na poprawkę $\delta _{c}$ dla $x_{{k+1}}^{{(1)}}$ .
Gruba poprawka (ang. coarse grid correction) — wyznaczenie poprawki $\delta _{c}$ na rzadszej siatce — przez iterację taką, jak ta! (Jako przybliżenie początkowe dla $\delta _{c}$ sensownie wziąć zero.)
Przedłużenie — wyrażenie poprawki obliczonej na rzadszej siatce w terminach wartości na gęstej siatki. Uwzględnienie poprawki i aktualizacja przybliżenia: $x_{{k+1}}^{{(2)}}=x_{{k+1}}^{{(1)}}+\delta$ .
Dodatkowe wygładzenie (ang. post-smoothing) — ewentualnie. Znów kilka iteracji jak na pierwszym etapie, z przybliżeniem początkowym równym $x_{{k+1}}^{{(2)}}$ , zakończone wyznaczeniem $x_{{k+1}}$ . Gdy ten etap jest pominięty, za $x_{{k+1}}$ bierze się $x_{{k+1}}^{{(2)}}$ .

Okazuje się, że w ten sposób uzyskujemy metodę, której szybkość zbieżności na modelowym zadaniu nie zależy od parametru dyskretyzacji $h$ (a więc i od rozmiaru zadania). Przypomnijmy, że było to zmorą wszystkich dotychczas rozważanych metod iteracyjnych; aby to przezwyciężyć, musieliśmy tam użyć spektralnie równoważnej macierzy ściskającej. Metoda wielosiatkowa nie potrzebuje tego, co więcej, można pokazać, że koszt jednej iteracji metody jest rzędu $O(N)$ — a więc (z dokładnością do stałej…) taki sam, jak koszt np. jednej iteracji metody CG.

Pełna metoda wielosiatkowa

Idąc dalej tym tropem możemy zauważyć, że jako dobre przybliżenie początkowe dla (iteracyjnej) metody wielosiatkowej można byłoby wybrać przybliżone rozwiązanie na rzadszej siatce: wtedy jedynymi istotnymi składowymi błędu na gęstej siatce byłyby składowe szybkozmienne — a te wyeliminowalibyśmy w kilku iteracjach! Oczywiście, wymagane przybliżenie na rzadszej siatce uzyskalibyśmy w ten sam sposób z jeszcze rzadszej siatki; na siatce o największych oczkach moglibyśmy użyć jakiejś metody bezpośredniej.

Pamiętając, że rozwiązywane przez nas zagadnienie — układ równań liniowych — pochodzi z dyskretyzacji pewnego równania różniczkowego musimy zdać sobie sprawę z tego, że nie ma sensu rozwiązywać go z dokładnością większą niż dokładność aproksymacji samej dyskretyzacji! Dlatego iterację można przerwać po osiągnięciu z góry zadanego poziomu błędu aproksymacji rozwiązania zadania różniczkowego. Gdy będziemy startować z przybliżenia wyznaczonego metodą opisaną powyżej, okazuje się, że często wystarczy tylko niewielka, niezależna od $h$ , liczba iteracji by osiągnąć zadany poziom błędu.

Taka metoda wyznaczania rozwiązania zadania powstałego po dyskretyzacji równania różniczkowego — którą tutaj tylko zasygnalizowaliśmy w uproszczeniu — nazywa się pełną metodą wielosiatkową (ang. full multigrid). Choć wciąż jest to metoda iteracyjna, ma ona cechy metody bezpośredniej. Okazuje się, że koszt wyznaczenia rozwiązania tą metodą jest rzędu $O(N)$ (kilku iteracji metody wielosiatkowej na finalnej, najdrobniejszej siatce) — a więc optymalny co do rzędu, gdyż samych niewiadomych w zadaniu jest $N$ .

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Matematyka obliczeniowa II

Zagadnienia

8. Przegląd innych metod rozwiązywania wielkich układów równań liniowych

8.1. Inne metody Kryłowa

8.1.1. Uogólnienia metody CG

Quiz 8.1

Twierdzenie 8.1

Dowód

Twierdzenie 8.2

Dowód

Ćwiczenie 8.1

Stwierdzenie 8.1

Dowód

Metoda Odir dla GCG, wersja bazowa

Uwaga 8.1

Ćwiczenie 8.2

8.1.2. Metody nie oparte na minimalizacji

8.2. Metody strukturalne

Przykład 8.1 (Strukturalizacja dyskretyzacji równania różniczkowego)

Ćwiczenie 8.3

8.3. Metody wielosiatkowe

Pełna metoda wielosiatkowa