2.3. Iterativne metode rješavanja sustava linearnih jednadžbi

Gaussove eliminacije i LU faktorizacija nisu jedini algoritmi za numeričko riješavanje linearnih sustava. Postoje još mnogi sofisticiraniji algoritmi, a u ovoj cjelini opisat ćemo još dva jednostavnija algoritma koja pripadaju familiji iterativnih metoda za rješavanje sustava.

Potreba za razvijanjem i korištenjem iterativnih metoda dolazi iz linearnih sustava kod kojih su matrice velike i rijetko popunjene. Takve matrice se u memoriji ne spremaju na klasičan način budući da imaju veliki broj elemenata jednakih nuli. Konkretno, umjesto da pohranjujemo matricu reda \(n\) kao dvodimenzionalno polje koje zauzima \(n^2\) realnih brojeva u memoriji, pamte se tri vektora \(v^{(i)}\), \(v^{(j)}\) i \(v^{(c)}\) duljine \(N\), gdje je \(N\) broj nenul elemenata matrice. Elementi na indexu \(k\) ovih vektora označavaju da se na poziciji \((v^{(i)}_k, v^{(j)}_k)\) u matrici nalazi broj \(v^{(c)}_k\), a na svim ostalim pozicijama matrice koje nisu navedene u vektorima nalaze se nule. Ako je u svakom retku i stupcu samo nekoliko (na primjer, uvijek 10 ili manje) elemenata matrice različito od nule, ovakav način pohrane velike matrice reda \(n\) koristi samo \(\bigO(n)\) memorije, što je mnogo manje od klasičnog spremanja u dvodimenzionalnom polju. Tada se i neke operacije s tako pohranjenim matricama mogu izvesti efikasnije: na primjer, algoritam množenja matrice vektorom tada ima složenost \(\bigO(n)\), a ne \(\bigO(n^2)\) kao inače. Ova činjenica se može iskoristiti za konstrukciju iterativnih metoda za rješavanje sustava. S druge strane, ako je matrica sustava rijetko popunjena i bez jasne strukture (nije vrpčasta; takav slučaj se pojavljuje kod matrica koje izlaze iz numeričkog rješavanja parcijalnih diferencijalnih jednadžbi), direktne metode poput LU faktorizacije neće iskoristiti tu rijetku popunjenost, što će uzrokovati korištenje previše memorije (ponovno \(n^2\) realnih brojeva za pohranjivanje \(L\) i \(U\)) i previše vremena za rješavanje.

Opišimo ideju iza iterativnih metoda iz ove cjeline. Zapišimo matricu sustava \(A\) u obliku \(A=M-N\), gdje su \(M\) i \(N\) neke matrice isto reda kao \(A\), te neka se \(M\) „jednostavno” invertira (bit će nam kasnije jasno što to znači). Pogledajmo niz jednakosti

\[\begin{split}\begin{align*} Ax &= b \\ (M-N)x&=b \\ Mx &= Nx + b \\ x &= M^{-1}Nx + M^{-1}b \end{align*}\end{split}\]

Zadnja jednakost motivira nas da uvedemo niz iteracija

\[ x^{(k+1)} = M^{-1}Nx^{(k)} + M^{-1}b, \]

uz neku početnu vrijednost \(x^{(0)}\). Primijetimo da ako postoji limes niza \(\hat{x} := \lim_{k} x^{(k)}\), puštanjem \(k\to \infty\) u gornjoj jednakosti dobivamo \(\hat{x} = M^{-1}N \hat{x} + M^{-1}b\), odnosno, možemo se vratiti unatrag do \(A \hat{x}=b\). Drugim riječima, ako niz iteracija konvergira, onda mu je limes upravo rješenje polaznog sustava, tj. \(\lim_k x^{(k)} = x\).

Uz oznake \(B:= M^{-1}N\), \(c:=M^{-1}b\), dobili smo rekurzivnu relaciju

(2.10)\[x^{(k+1)} = Bx^{(k)} + c.\]

koju ćemo proučavati u nastavku. Uočimo da za rješenje linearnog sustava \(Ax = b\) vrijedi \(x = Bx + c\), kao i za limes niza iteracija (ako niz konvergira).

Jacobijeva metoda

Da bi metoda zaista bila lakša za implementaciju od ostalih metoda za rješavanja sustava \(Ax=b\), matrica \(M\) iz rastava \(A=M-N\) mora biti jednostavna za invertiranje, odnosno, račun članova \(M^{-1}N\) i \(M^{-1}b\) ne smije biti kompliciraniji od računa \(A^{-1}b\), s kojim smo krenuli. Uzmemo li da je \(M\) dijagonalna, koja jest jednostavna za invertiranje, dobit ćemo Jacobijevu metodu za rješavanje sustava. Dakle, \(M:=D\), gdje je \(D\) dijagonalni dio matrice \(A\), te je \(N := M-A = D-A\) (odnosno vandijagonalni dio matrice \(A\) s minus predznakom).

Definicija 2.24 (Jacobijeva metoda)

Jacobijeva metoda za rješavanje sustava \(Ax=b\) glasi: za bilo koju početnu iteraciju \(x^{(0)}\), do određenog uvjeta zaustavljanja, ponavljati postupak

\[\begin{split}x^{(k+1)} = \begin{bmatrix} 0 & \dfrac{-a_{1,2}}{a_{1,1}} & \dots & \dfrac{-a_{1,n}}{a_{1,1}} \\ \dfrac{-a_{2,1}}{a_{2,2}}& 0 & \dots & \dfrac{-a_{2,n}}{a_{2,2}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{-a_{n,1}}{a_{n,n}}& \dfrac{-a_{n,2}}{a_{n,n}} & \dots & 0 \end{bmatrix} x^{(k)} + \begin{bmatrix} \dfrac{b_1}{a_{1,1}} \\ \dfrac{b_2}{a_{2,2}} \\ \vdots \\ \dfrac{b_n}{a_{n,n}}. \end{bmatrix}\end{split}\]

Zapisano po komponentama, u \(k\)-tom koraku (\(k=1,2,3,\dots\)) ažuriranje \(i\)-te komponente (\(i=1,\dots,n\)) glasi

(2.11)\[\begin{split}x^{(k+1)}_i = \frac{1}{a_{i,i}} \left( b_{i} - \sum_{\substack{j=1 \\ i\not=j}}^n a_{i,j} x_j^{(k)} \right).\end{split}\]

Gauss-Seidelova metoda

Pogledajmo ponovno (2.11). Prilikom računa \(i\)-te komponente novog vektora rješenja, u sumi s desne strane između ostalog uzimamo prvih \((i-1)\) komponenti starog rješenja iako na tim pozicijama znamo i vrijednosti novog vektora rješenja. Ako želimo metoda konvergira što brže, i ako vjerujemo da su vrijednosti komponenti \(x^{(k)}_i\) sve bliže komponentama rješenja \(x_i\), možemo metodu učiniti boljom tako da tih prvih \((i-1)\) komponenti računamo iz rješenja \(x^{k+1}\). Preciznije, možemo prilagoditi iteracije da u \(k\)-tom koraku (\(k=1,2,3,\dots\)) ažuriranje \(i\)-te komponente (\(i=1,\dots,n\)) glasi

(2.12)\[x^{(k+1)}_i = \frac{1}{a_{i,i}} \left( b_{i} - \sum_{j=1}^{i-1} a_{i,j} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} \right).\]

Time dobivamo Gauss-Seidelovu metodu za rješavanje linearnog sustava.

Kako bi se uvjerili da ona dolazi iz iste familije iterativnih sustava kao i Jacobijeva (tj. da postoje matrice \(M\) i \(N\) s početka lekcije), zapišimo ove iteracije na ekvivalentan način:

\[\begin{split}\begin{align*} x^{(k+1)}_i &= \frac{1}{a_{i,i}} \left( b_{i} - \sum_{j=1}^{i-1} a_{i,j} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} \right) \\ a_{i,i} x^{(k+1)}_i &= b_{i} - \sum_{j=1}^{i-1} a_{i,j} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} \\ \sum_{j=1}^{i-1} a_{i,j} x_j^{(k+1)} + a_{i,i} x^{(k+1)}_i &= b_{i} - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} \\ \sum_{j=1}^{i} a_{i,j} x_j^{(k+1)} &= - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} + b_i. \end{align*}\end{split}\]

Sada vidimo da ako je \(M=L\) i \(N=L-A\), gdje je \(L\) donjetrokutasti dio matrice \(A\) koji uključuje i dijagonalu (a \(N\) zapravo dio matrice \(A\) strogo iznad dijagonale uz promijenjen predznak svim elementima), dobivamo jednakost \(Mx^{(k+1)} = Nx^{(k)} + b\). To je u skladu s idejom iterativnih metoda iz ove lekcije.

Definicija 2.25 (Gauss-Seidelova metoda)

Gauss-Seidelova metoda za rješavanje sustava \(Ax=b\) glasi: za bilo koju početnu iteraciju \(x^{(0)}\), do određenog uvjeta zaustavljanja, ponavljati postupak

\[x^{(k+1)}_i = \frac{1}{a_{i,i}} \left( b_{i} - \sum_{j=1}^{i-1} a_{i,j} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{i,j} x_j^{(k)} \right)\]

(\(i=1,2,\dots,n\), \(k =1,2,\dots\)).

Konvergencija iterativnih metoda

Reći ćemo da iterativna metoda konvergira ako za svaki početni \(x^{(0)}\) niz iteracija \((x^{(k)})_{k\geq 0}\) definiranih s (2.10) konvergira k rješenju sustava \(Ax = b\). Inače kažemo da divergira.

Da bismo analizirali konvergenciju metode, promotrimo razliku \(e^{(k)} := x^{(k)} - x\). Kako vrijedi \(x = Bx + c\), oduzimanjem od (2.10) dolazimo do

(2.13)\[e^{(k+1)} = B e^{(k)}.\]

Primijetite da je to kao da promatramo iteracije na homogenom sustavu (kao da je \(b=0\)).

Propozicija 2.26

Ako je \(\| B \| < 1\), gdje je \(\| \cdot \|\) operatorska norma, tada iterativna metoda konvergira.

Dokaz.

Promatrajući (2.13), dobivamo

\[\|e^{(k+1)}\| = \|B e^{(k)} \| \leq \|B\| \|e^{(k)} \|.\]

Ovu nejednakost sada možemo „odmotati” skroz do \(k=0\):

\[\|e^{(k)}\| \leq \|B\| \|e^{(k-1)}\| \leq \|B\| \|B\| \|e^{(k-2)}\| \leq \ldots \leq \|B\|^k \|e^{(0)} \|.\]

Ako je \(\|B\|<1\), kad \(k\to \infty\) dobivamo da nužno \(\|e_{k}\| \to 0\), čime je dokaz tvrdnje gotov.

Slično kao u dokazu možemo izvesti uvjet zaustavljanja iterativne metode: ako je cilj s iteracijama stati kada je \(\|x^{(k)}-x\| < \varepsilon\), može se pokazati da je za to dovoljan uvjet stati kada razlika uzastopnih iteracija zadovoljava

\[\frac{\|B\|}{1-\|B\|} \| x^{(k)} - x^{(k-1)} \| < \varepsilon.\]

Propozicija 2.26 nije nužan i dovoljan uvjet za konvergenciju iteratine metode. To ćemo pokazati u sljedećim primjerima.

Primjer 2.27

Promotrimo sustav u kojem je matrica sustava

\[\begin{split}A = \begin{bmatrix} 1 & \alpha \\ \alpha & 1 \end{bmatrix},\end{split}\]

za neki \(\alpha \geq 0\). U slučaju Jacobijeve metode, matrica \(B\) glasi

\[\begin{split}B = \begin{bmatrix} 0 & -\alpha \\ -\alpha & 0 \end{bmatrix}.\end{split}\]

Lako je vidjeti tu matricu vrijedi \(\|B\|_1 = \|B\|_{\infty} = \alpha\). Čak i za \(p=2\), kako su obje svojstvene vrijednosti matrice

\[\begin{split}B^TB = \begin{bmatrix} 0 & -\alpha \\ -\alpha & 0 \end{bmatrix} \begin{bmatrix} 0 & -\alpha \\ -\alpha & 0 \end{bmatrix} = \begin{bmatrix} \alpha^2 & 0 \\0 & \alpha^2 \end{bmatrix}\end{split}\]

jednake \(\alpha^2\), vrijedi \(\|B\|_2 = \alpha\).

Pogledajmo za koje \(\alpha\) metoda konvergira. Primijetimo da vrijedi

\[\begin{split}B^2 = \begin{bmatrix} 0 & -\alpha \\ -\alpha & 0 \end{bmatrix} \begin{bmatrix} 0 & -\alpha \\ -\alpha & 0 \end{bmatrix} = \begin{bmatrix} \alpha^2 & 0 \\0 & \alpha^2 \end{bmatrix} = \alpha^2 I.\end{split}\]

Ako uvedemo oznaku \(\displaystyle J:=\begin{bmatrix} 0 & -1 \\ -1 & 0 \end{bmatrix}\), dobivamo

\[\begin{split}e^{(k)} = B^k e^{(0)} = \begin{cases} \alpha^k e^{(0)} & k \text{ paran}, \\ \alpha^k Je^{(0)} & k \text{ neparan}. \end{cases}\end{split}\]

Odavde vidimo da \((e^{(k)})_k\) konvergira k nuli samo ako \(\alpha < 1\).

Primjer 2.28

Promotrimo sustav u kojem je matrica sustava

\[\begin{split}A = \begin{bmatrix} 1 & \alpha \\ 0 & 1 \end{bmatrix},\end{split}\]

za neki \(\alpha \geq 0\). U slučaju Jacobijeve metode, matrica \(B\) glasi

\[\begin{split}B = \begin{bmatrix} 0 & -\alpha \\ 0 & 0 \end{bmatrix}.\end{split}\]

Lako je vidjeti tu matricu vrijedi \(\|B\|_1 = \|B\|_{\infty} = \alpha\). Čak i za \(p=2\), kako su svojstvene vrijednosti matrice

\[\begin{split}B^TB = \begin{bmatrix} 0 & 0 \\ -\alpha & 0 \end{bmatrix} \begin{bmatrix} 0 & -\alpha \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\0 & \alpha^2 \end{bmatrix}\end{split}\]

jednake \(\alpha^2\) i \(0\), vrijedi \(\|B\|_2 = \alpha\).

Pogledajmo za koje \(\alpha\) metoda konvergira. Primijetimo da vrijedi

\[\begin{split}B^2 = \begin{bmatrix} 0 & -\alpha \\ 0 & 0 \end{bmatrix} \begin{bmatrix} 0 & -\alpha \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\0 & 0 \end{bmatrix}.\end{split}\]

Odavde zaključujemo da za svaki \(\alpha\) i za svaki početni \(x^{(0)}\) dobivamo točno rješenje u dva koraka, budući da je \(e^{(2)} = B^2 e^{(0)} = 0\).

Gornji primjeri pokazuju da za nužan i dovoljan uvjet konvergencije iterativne metode ne možemo gledati samo normu matrice \(B\). Razlika u primjerima nije u svojstvenim vrijednosti matrice \(B^TB\), nego u svojstvenim vrijednostima matrice \(B\). To motivira sljedeći teorem koji navodimo bez dokaza.

Teorem 2.29

Iterativna metoda konvergira ako i samo ako je \(\rho(B)<1\).

Definicija 2.30

Za kvadratnu matricu \(A \in \R^{n \times n}\) spektralni radijus \(\rho(A)\) definira se kao najveća apsolutna vrijednost svojstvenih vrijednosti:

\[\rho (A) = \max \{ |\lambda| : \lambda \in \sigma(A)\}.\]

Primijetimo da je u Primjeru 2.27 \(\rho(B) = \alpha\), a u Primjeru 2.28 \(\rho(B) = 0\).

Postoje familije matrica za koje je lakše utvrditi konvergira li neka od iterativnih metoda.

Teorem 2.31

Neka je \(A\) kvadratna matrica koja je strogo dijagonalno dominantna po retcima ili stupcima. Tada i Jacobijeva i Gauss-Seidelova metoda konvergiraju.

Dokaz nećemo raditi, napomenimo da je za Jacobijevu metodu možete dokazati sami tako da dokažete da je 1-norma ili \(\infty\)-norma matrice \(B\) manja od \(1\).