Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Matematyka obliczeniowa II – 8. Przegląd innych metod rozwiązywania wielkich układów równań liniowych – MIM UW

Zagadnienia

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)
Quiz 8.1
  1. Czy metoda CG to nic innego jak GCG(I, I, A)?



    Dobrze Oczywiście, CG to nic innego jak GCG(A, I, A)

    Źle Sprawdź, w jakiej normie CG minimalizuje błąd

  2. Czy metoda PCG z macierzą ściskającą M to nic innego jak GCG(A, M, A)?



    Dobrze Oczywiście

    Źle Zob. [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 relacje10Taką 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:

  1. Oblicz q_{i}=A_{{ii}}^{{-1}}b_{i}, i=1,2.

  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}.

  3. 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}:

  1. 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.

  2. Restrykcja — sformułowanie na rzadszej siatce zadania na poprawkę \delta _{c} dla x_{{k+1}}^{{(1)}}.

  3. 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.)

  4. 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.

  5. 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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.