Algorytmy Jacobiego i Gaussa-Seidla

Algorytmy Jacobiego i Gaussa-Seidla, które w niniejszej pracy zawsze rozwa-żane są w wersjach blokowych [154], są reprezentantami klasycznych metod iteracji prostej. Pojedyncza iteracja blokowej metody Jacobiego lub Gaussa--Seidla zastosowana do rozwiązania układu równań liniowych (2.41) może zo-stać zapisana jako sekwencja problemów lokalnych postaci:

A_ll· u^m,l_l = −

Nbl

X

k=1,k6=l

A_lk· ˜u_k+ b_l (2.56)

Powyżej u_loznacza podwektor stopni swobody, będący częścią globalnego wek-tora u, związany z pojedynczym obiektem siatki, dowolną grupą obiektów lub grupą obiektów tworzących klasyczny podobszar Ω_lobszaru obliczeniowego Ω.

Podobszary Ω_l mogą zachodzić na siebie (jeśli nie prowadzi to do zachodze-nia na siebie podwektorów u_l). N_bl oznacza liczbę podwektorów u_l. Bloki A_lk odpowiadają parom podwektorów u_l, a co za tym idzie parom obiektów siatki, grup obiektów siatki lub podobszarów, zależnie od definicji u_l. Podwek-tory b_l, fragmenty globalnego wektora b, odpowiadają podwektorom u_l. u^m,l_l oznacza podwektor u_l w trakcie wykonywania m-tej iteracji, po rozwiązaniu l problemów lokalnych (2.56). ˜u jest aktualnym rozwiązaniem, uzyskiwanym różnie dla różnych metod. Dla iteracji Jacobiego uaktualnienie rozwiązania następuje po rozwiązaniu wszystkich problemów lokalnych, stąd ˜u_k = u^m,0_k . Dla iteracji Gaussa-Seidla uaktualnienia dokonuje się po rozwiązaniu każdego problemu lokalnego, tak więc ˜u_k = u^m,l−1_k . W obu wypadkach u^m,0 = u^m−1 jest wektorem początkowym dla m-tej iteracji.

W klasycznej analizie metod iteracji prostej [85] zapisywane są one w po-staci:

u^m= S⁻¹Q · u^m−1− S⁻¹b (2.57) gdzie S i Q uzyskiwane są w efekcie rozkładu macierzy A:

A = S − Q (2.58)

Dla metod Jacobiego i Gaussa-Seidla wykorzystuje się rozkład na części: trój-kątną dolną L, diagonalną D i trójtrój-kątną górną U. Pojedyncza iteracja metody

Jacobiego zapisywana jest jako:

u^m = −D⁻¹(L + U) · u^m−1+ D⁻¹b (2.59) natomiast pojedyncza iteracja metody Gaussa-Seidla przybiera postać:

u^m= −(D + L)⁻¹U · u^m−1+ (D + L)⁻¹b (2.60) Metody iteracji prostej są zbieżne, jeśli promień spektralny macierzy S⁻¹Q jest mniejszy od 1.

W klasycznych sformułowaniach blokowych wariantów metod Jacobiego i Gaussa-Seidla podwektory u_li odpowiadające im bloki macierzy A nie zacho-dzą na siebie. Jeśli zachodzenie występuje, nadal można stosować sekwencje rozwiązań problemów (2.56) dla iteracyjnego rozwiązania układu (2.41). Wy-nikające z tej techniki metody są uogólnieniami iteracji prostych, należącymi do metod Schwarza, czyli metod dekompozycji obszaru na nakładające się podobszary [162].

W analizie metod Schwarza wprowadza się dwa rodzaje operatorów: opera-tory zawężenia (restrykcji) R_l i operatory rozszerzenia (prolongacji) R^T_l . Ope-rator zawężenia R_l wybiera spośród wszystkich składowych wektora u pod-wektor stopni swobody u_l (w klasycznych metodach Schwarza odpowiadający podobszarowi Ω_l). Operator rozszerzenia R^T_l uzupełnia podwektor stopni swo-body u_l zerami, do pełnego wymiaru wektora u. W zapisie macierzowym R_l jest tablicą złożoną z jedynek i zer, a R^T_l jest, zgodnie z oznaczeniem, trans-pozycją R_l, przy czym R_lR^T_l = I.

Wykorzystując powyższą notację, lokalne problemy (2.56) dają się zapisać jako:

R_lu^m,l = −A⁻¹_ll





Nbl

X

k=1,k6=l

A_lkR_ku − R˜ _lb



= (2.61)

= R_lu − A˜ ⁻¹_ll

Nbl

X

k=1

A_lkR_ku + A˜ ⁻¹_ll R_lb gdzie w celu uzyskania ostatniej równości zastosowano tożsamość:

R_lu = A˜ ⁻¹_ll A_llR_lu˜ (2.62) Ostateczny wektorowy zapis problemu (2.56) w notacji metod dekompozycji obszaru wygląda następująco:

u^m,l= ˜u + R^T_l A⁻¹_ll Rl(b − A˜u) (2.63)

W konsekwencji, pojedyncza iteracja addytywnej metody Schwarza, będącej uogólnieniem blokowej metody Jacobiego, uzyskuje formę:

u^m = u^m−1− natomiast pojedyncza iteracja multiplikatywnej metody Schwarza, będącej uogólnieniem blokowej metody Gaussa-Seidla, przybiera postać:

u^m = u^m−1− Kluczem do analizy metod dekompozycji obszaru jest pojęcie rzutowania błędu aktualnego rozwiązania na lokalną przestrzeń (wektorową lub funkcyjną) związaną z podobszarem (podwektorem stopni swobody u_l). Operator projek-cji na lokalną przestrzeń wektorową definiuje się wzorem:

P_l= R^T_l A⁻¹_ll R_lA (2.66) Jeśli przez u^∗ oznaczone zostanie rozwiązanie dokładne układu równań, tak że b = A · u^∗, to lokalny problem (2.63) okazuje się być korektą aktualnego rozwiązania ˜u o zrzutowany na lokalną przestrzeń globalny błąd rozwiązania

˜ e:

u^m,l= ˜u + Pl(u^∗− ˜u) = ˜u + Pl˜e (2.67) W dalszych rozważaniach proces iteracyjny zapisywany jest w języku ana-lizy funkcjonalnej. Słabe sformułowanie (2.30) lub (2.31) związane z rozwią-zywanym układem równań, po uproszczeniu zapisu przez pominięcie indeksów konkretnych metod, przybiera formę:

Znajdź taką funkcję u ∈ V_h, aby spełnione było:

a(u, v) = l(v) ∀v ∈ V_h (2.68)

Rozważa się dekompozycję obszaru obliczeniowego Ω na nakładające się podobszary Ω_l i lokalne problemy projekcji na przestrzenie funkcyjne V_l⊂ V_h, związane z podobszarami Ω_l. Dla każdego podobszaru przestrzeń V_l jest roz-pięta na lokalnych dla Ω_l funkcjach bazowych zerujących się na brzegu we-wnętrznym podobszaru²⁰. Dla dowolnej funkcji u_l ∈ V_l określa się jej roz-szerzenie u⁰_l do pełnej przestrzeni V_h (poprzez funkcję zerową). Dla dowolnej

20Brzegiem wewnętrznym podobszaru jest ta część jego brzegu, która nie pokrywa się z brzegiem obszaru obliczeniowego.

funkcji g problem rzutowania związany ze sformułowaniem (2.68) i rozwiązy-wanym układem równań ma postać:

Znajdź taką funkcję g_l∈ V_l, aby spełnione było:

a_l(g_l, v_l) = a(g, v_l⁰) ∀v_l∈ V_l (2.69) Lokalną formę dwuliniową a_l(., .), odpowiadającą lokalnym problemom (2.56), uzyskuje się z formy globalnej poprzez zastąpienie obszaru Ω podobszarem Ω_l. Jeśli funkcję z przestrzeni V_l odpowiadającą aktualnemu rozwiązaniu ˜u oznaczymy przez ˜u, a funkcję odpowiadającą rozwiązaniu układu u^∗ przez u^∗, to dla funkcji błędu ˜e = u^∗− ˜u zachodzi równość:

a(˜e, v_l⁰) = a(u^∗− ˜u, v_l⁰) = l(v⁰_l) − a(˜u, v⁰_l) (2.70) Zagadnienie lokalnej projekcji globalnego błędu rozwiązania, zapisywane dla wektorów jako:

˜

e_l= P_l˜e (2.71)

dla odpowiadających funkcji ma ostatecznie postać:

Znajdź taką funkcję ˜e_l∈ V_l, aby spełnione było:

a_l(˜e_l, v_l) = l(v_l⁰) − a(˜u, v_l⁰) ∀v_l∈ V_l (2.72) Podobnie jak dla wektora niewiadomych, także dla odpowiadającej mu funkcji lokalne problemy przyjmują postać korekty aktualnych wartości za pomocą błędu zrzutowanego na lokalną przestrzeń:

u^m,l= ˜u + ˜e_l (2.73)

Zinterpretowane jako metody dekompozycji obszaru klasyczne metody blo-kowe Jacobiego i Gaussa-Seidla ukazują się jako serie lokalnych projekcji błędu i następujących (w różnych momentach dla różnych metod) korekt aktual-nego rozwiązania. Istotnym elementem w rozważanych metodach dekompozy-cji (metodach Schwarza) jest nakładanie się podobszarów, z którymi związane są lokalne podprzestrzenie. Zbieżność iteracji jest w tych metodach zależna od rozmiaru podobszarów, ich kształtu, rozmiaru strefy nakładania się podobsza-rów i, oczywiście, od własności problemu [67, 66, 68, 145]. W standardowych przypadkach, im większe podobszary i im większa strefa nakładania się podob-szarów, tym zbieżność metod szybsza. Zbieżność ta jest jednak na tyle wolna, że w praktycznych zastosowaniach zawsze stosuje się techniki jej przyspiesza-nia.