Sustavi nelinearnih jednadžbi

Sadržaj

6.2. Sustavi nelinearnih jednadžbi#

Neke od ideja razvijenih za rješavanje nelinearnih jednadžbi s jednom nepoznanicom ćemo moći poopćiti i na sustave nelinearnih jednadžbi s \(n\) nepoznanica. Ako nepoznanice označimo sa \(x_1, \ldots, x_n\), a postoji \(m\) jednadžbi koje te nepoznanice moraju zadovoljavati, onda svaku od tih jednadžbi možemo zapisati ovako:

\[ f_i(x_1, x_2, \ldots, x_n) = 0, \quad i=1, 2, \ldots, m, \]

za neke funkcije \(f_1, f_2, \ldots, f_m : \R^n \to \R\). Označimo li

\[\begin{split} f(x_1, x_2, \ldots, x_n) = \mb{c} f_1(x_1, x_2, \ldots, x_n) \\ f_2(x_1, x_2, \ldots, x_n) \\ \vdots \\ f_m(x_1, x_2, \ldots, x_n) \me, \end{split}\]

vidimo da nas zapravo zanimaju nultočke funkcije

\[ f : \R^n \to \R^m. \]

Do sada smo u ovom i drugim kolegijima promatrali samo sustave linearnih jednadžbi: za zadanu matricu \(A \in \R^{m \times n}\) i vektor \(b \in \R^m\) tražili smo sve vektore \(x \in \R^n\) koji su nultočke funkcije \(f(x) = Ax - b\). Sada će nas prvenstveno zanimati nelinearne funkcije \(f\).

Radi jednostavnosti, u nastavku uvodimo sljedeće pretpostavke:

Broj jednadžbi je jednak broju nepoznanica, tj. \(m=n\).
Funkcija \(f\) ima samo izolirane nultočke \(\alpha \in \R^n\).
Funkcija \(f\) je dovoljno glatka u okolini nultočke koju želimo odrediti.

Promotrimo sada metode kojima smo rješavali nelinearne jednadžbe s jednom nepoznanicom i pokušajmo ih poopćiti na sustave nelinearnih jednadžbi. Uočimo prvo da je to poopćenje vrlo teško napraviti za metodu bisekcije. Zamislimo da tražimo nultočke funkcije \(f : \R^2 \to \R^2\).

Kako poopćiti uvjet \(f(a) f(b) < 0\) koji jamči da je nultočka \(\alpha \in [a, b]\), a koji smo imali u 1d-slučaju? Kako bi glasio analogni uvjet za 2d-slučaj i pravokutnik \([x_1, x_2] \times [y_1, y_2]\) (ili možda neki drugi ograničeni podskup od \(\R^2\), koji?)
Ako za točku \((x_\ast, y_\ast)\) utvrdimo npr. \(f(x_\ast, y_\ast) < 0\), kako nastaviti proces bisekcije? Nije moguće pronaći manji pravokutnik unutar kojeg se nalazi nultočka!

Ova pitanja se samo još dodatno kompliciraju u općoj situaciji sa \(n\) jednadžbi i \(n\) nepoznanica. S druge strane, Newtonovu metodu i metodu jednostavnih iteracija možemo poopćiti i na funkcije više varijabli.

Newtonova metoda#

Kao i u slučaju \(n=1\), Newtonovu metodu dobivamo linearizacijom funkcije \(f\). Pretpostavimo da u \(k\)-tom koraku imamo aproksimaciju \(x^{(k)} \in \R^n\) za nultočku funkcije \(f\). Promotrimo Taylorov razvoj \(i\)-te komponente funkcije \(f\) oko točke \(x^{(k)}\), te odbacimo sve članove nakon linearnih:

(6.3)#\[ f_i(x) \approx f_i(x^{(k)}) + \sum_{j=1}^n \frac{\partial f_i}{\partial x_j} (x^{(k)}) \cdot (x_j - x_j^{(k)}), \quad i = 1, 2, \ldots n.\]

Ovdje smo označili

\[ x = (x_1, x_2, \ldots, x_n), \quad x^{(k)} = (x_1^{(k)}, x_2^{(k)}, \ldots, x^{(k)}_n). \]

Prisjetimo se definicije Jacobijeve matrice:

\[\begin{split} J_f(x) = \frac{\partial f}{\partial x}(x) = \mb{cccc} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \ldots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \ldots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \ldots & \frac{\partial f_n}{\partial x_n} \\ \me \in \R^{n \times n}. \end{split}\]

Izraz (6.3) možemo kompaktno zapisati kao

\[ f(x) \approx f(x^{(k)}) + J_f(x^{(k)}) \cdot (x - x^{(k)}). \]

Kao i u slučaju nelinearnih jednadžbi s jednom nepoznanicom, novu aproksimaciju \(x^{(k+1)}\) dobivamo iz zahtjeva da je desna strane gornje približne jednakosti jednaka nuli jer tada očekujemo \(f(x^{(k+1)}) \approx 0\):

\[ f(x^{(k)}) + J_f(x^{(k)}) \cdot (x^{(k+1)} - x^{(k)}) = 0. \]

Stavimo \(d^{(k)} := x^{(k+1)} - x^{(k)}\). Ako je \(J_f(x^{(k)})\) regularna matrica, onda je \(d^{(k)}\) jedinstveno rješenje sustava linearnih jednadžbi

(6.4)#\[ J_f(x^{(k)}) \cdot d^{(k)} = -f(x^{(k)}).\]

Slijedi da za iduću iteraciju u Newtonovoj metodi vrijedi

(6.5)#\[ x^{(k+1)} = x^{(k)} - \left( J_f(x^{(k)}) \right)^{-1} f(x^{(k)}).\]

Uočimo da je dobivena formula zaista poopćenje „skalarne” Newtonove metode kod koje smo imali

\[ x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})}, \]

jer se u slučaju \(n=1\) Jacobijeva matrica \(J_f(x^{(k)})\) podudara s derivacijom \(f'(x^{(k)})\).

Upozorenje

Kod implementacije Newtonove metode na računalu, sljedeću iteraciju ne računamo direktnim korištenjem formule (6.5), odnosno, ne računamo inverz Jacobijeve matrice.

Umjesto računanja inverza matrice (što se uvijek izbjegava), rješavamo sustav linearnih jednadžbi (6.4) i onda pomoću dobivene korekcije \(d^{(k)}\) računamo sljedeću iteraciju:

\[ x^{(k+1)} = x^{(k)} + d^{(k)}. \]

Slijedi jednostavna implementacija Newtonove metode, koju zaustavljamo nakon k_max iteracija ili kada niz iteracija postane „dovoljno stacionaran”, tj. kada vrijedi

\[ \max_{i=1, \ldots, n} \frac{|x_i^{(k+1)} - x_i^{(k)}|}{|x_i^{(k)}|} \leq \eps. \]

Ovaj kriterij zaustavljanja podsjeća na kriterij \(|x^{(k+1)} - x^{(k)}| \leq \sqrt{\frac{2m_1}{M_2} \eps}\) koji smo imali kod skalarne Newtonove metode no nećemo ga ovdje detaljnije diskutirati.

import numpy as np;
import numpy.linalg as nla;

def newton_sustav( f, Jf, x0, k_max, eps ):
# Ulaz:
#   - f = funkcija čiju nultočku tražimo.
#   - Jf = funkcija koja vraća Jacobijevu matricu od f u danoj točki.
#   - x0 = početna iteracija
#   - k_max = broj iteracija Newtonove metode
# Izlaz:
#   - x = aproksimacija nultočke od f dobivena nakon k_max iteracija Newtonove metode 
#         ili nakon što niz postane stacionaran.
#   - hist = polje svih izračunatih iteracija (treba samo za tablični prikaz kasnije).

    x = x0; hist = [x0];
    for k in range( 0, k_max ):
        xk = x;
        
        # Riješimo linearni sustav Jf(x_k)*d = -f(x_k).
        dk = nla.solve( Jf(xk), -f(xk) );

        # Iduća Newtonova iteracija = korekcija prethodne iteracije.
        x = xk + dk;

        # Spremi novu iteraciju u hist.
        hist.append( x );

        # Provjeri je li niz postao "dovoljno stacionaran".
        if( np.max( np.abs( dk ) / np.abs( xk ) ) < eps ):
            break;

    return (x, hist);

Rezultati o lokalnoj konvergenciji kao i o brzini konvergencije Newtonove metode su posve analogni rezultatima koje smo izveli u skalarnom slučaju.

Teorem 6.18 (Lokalna konvergencija i brzina konvergencije Newtonove metode)

Neka je \(\alpha\) nultočka funkcije \(f : \R^n \to \R^n\). Pretpostavimo da je \(f\) klase \(C^2\) u nekoj okolini od \(\alpha\) i da je \(J_f(\alpha)\) regularna matrica.

Tada postoji okolina nultočke \(\alpha\) takva da:

Za svaku početnu iteraciju \(x^{(0)}\) iz te okoline, niz iteracija generiran Newtonovom metodom konvergira prema \(\alpha\).
Konvergencija je barem kvadratna, tj. vrijedi

\[ \|\alpha - x^{{(k+1)}}\| \in \bigO\left( \|\alpha - x^{(k)}\|^2 \right). \]

Dokaz.

Dokaz je posve analogan onom u skalarnom slučaju, no treba obratiti pažnju na neke detalje koji postaju netrivijalni jer se umjesto skalara pojavljuju matrice. Skica dokaza:

Taylorova formula u višedimenzionalnom slučaju se može zapisati kao

\[ f(x) = f(x^{(k)}) + J_f(x^{(k)}) \cdot (x - x^{(k)}) + R(x, x^{(k)}), \]

gdje je

\[\begin{split} \begin{align*} \|R(x, x^{(k)})\| &\leq \frac{M_2}{2} \|x - x^{(k)}\|^2, \\ \sum_{i, j, k} |\frac{\partial f_k}{\partial x_i \partial x_j}(z)|^2 &\leq M_2, \text{ za sve } \|z - x^{(k)}\| \leq \|x - x^{(k)}\|. \end{align*}\end{split}\]
Ako je \(J_f(\alpha)\) regularna, onda postoji \(\eps > 0\) takav da je \(J_f(x)\) regularna za svaki \(\|x - \alpha\| < \eps\). To slijedi zato što je \(J_f\) klase \(C^1\), pa je onda i Lipschitzova upravo s konstantom \(M_2\). Stavimo \(\eps = \sigma_n / (2M_2)\), gdje je \(\sigma_n > 0\) najmanja singularna vrijednost od \(J_f(\alpha)\). Tada je \(\|J_f(x) - J_f(\alpha)\| \leq M_2 \|x - \alpha\| \leq \sigma_n/2\). Slijedi da je \(J_f(x)\) punog ranga, tj. regularna, jer u protivnom imamo kontradikciju s Teoremom 4.11 (Eckart, Young, Mirsky).
Štoviše, može se pokazati da iz \(\|J_f(x) - J_f(\alpha)\| \leq \sigma_n/2\) slijedi \(\|J_f(x)^{-1}\| \leq \frac{2}{\sigma_n} = 2 \|J_f(\alpha)^{-1}\|\).

Stoga je \(x^{(k+1)}\) dobro definiran i zbog \(f(\alpha) = 0\) imamo:

\[ x^{(k+1)} - \alpha = x^{(k)} - \alpha - J_f(x^{(k)})^{-1} ( f(x^{(k)}) - f(\alpha) ), \]

iz čega izvodimo ocjenu

\[\begin{split} \begin{align*} \|x^{(k+1)} - \alpha\| &= \|x^{(k)} - \alpha - J_f(x^{(k)})^{-1} ( f(x^{(k)}) - f(\alpha) )\| \\ &= \|J_f(x^{(k)})^{-1} \cdot \left( J_f(x^{(k)}) (x^{(k)} - \alpha) - (f(x^{(k)}) - f) \right) \| \\ &\leq \|J_f(x^{(k)})^{-1}\| \cdot \| J_f(x^{(k)}) (x^{(k)} - \alpha) - (f(x^{(k)}) - f(\alpha)) \| \\ &\leq \frac{2}{\sigma_n} \cdot \| f(\alpha) - f(x^{(k)}) - J_f(x^{(k)}) (\alpha - x^{(k)})\| \\ &= \frac{2}{\sigma_n} \cdot \|R(\alpha, x^{(k)})\| \\ &\leq \frac{M_2}{\sigma_n} \|x^{(k)} - \alpha\|^2. \end{align*} \end{split}\]

Dakle, ako metoda konvergira, onda je konvergencija kvadratična. No metoda konvergira jer gornju nejednakost možemo nastaviti ovako:

\[\begin{split} \begin{align*} \|x^{(k+1)} - \alpha\| &\leq \frac{M_2}{\sigma_n} \|x^{(k)} - \alpha\| \cdot \|x^{(k)} - \alpha\| \\ &\leq \frac{M_2}{\sigma_n} \eps \cdot \|x^{(k)} - \alpha\| \\ &= \frac{M_2}{\sigma_n} \frac{\sigma_n}{2 M_2} \cdot \|x^{(k)} - \alpha\| \\ &= \frac{1}{2} \cdot \|x^{(k)} - \alpha\|, \end{align*} \end{split}\]

pa niz \(x^{(k)}\) konvergira prema \(\alpha\).

Uvjet o regularnosti matrice \(J_f(\alpha)\) iz gornjeg teorema za \(n=1\) glasi \(f'(\alpha) \neq 0\), odnosno, to je poopćenje uvjeta da je \(f\) jednostruka nultočka.

Primjer#

Riješimo sljedeći sustav nelinearnih jednadžbi:

\[\begin{split} \begin{align*} x_1^6 - 5x_1^2 x_2^2 + 136 &= 0, \\ x_2^4 - 3x_1^4 x_2 + 80 &= 0. \end{align*} \end{split}\]

Dakle, funkcija \(f\) i njezina Jacobijeva matrica su dani formulama

\[\begin{split} f(x_1, x_2) = \mb{c} x_1^6 - 5x_1^2 x_2^2 + 136 \\ x_2^4 - 3x_1^4 x_2 + 80 \me, \quad J_f(x_1, x_2) = \mb{cc} 6x_1^5 - 10x_1 x_2^2 & -10x_1^2 x_2 \\ -12 x_1^3 x_2 & 4x_2^3 - 3x_1^4 \me. \end{split}\]

Koristimo Newtonovu metodu s početnom točkom \(x^{(0)} = (1, 2)\).

def f(x):
# Izračunavamo i vraćamo vrijednost funkcije f u točki x.    
    x1 = x[0]; x2 = x[1];
    fx = np.array( [
        x1**6 - 5*(x1**2)*(x2**2) + 136, 
        x2**4 - 3*(x1**4)*x2 + 80
        ] );

    return fx;

def Jf(x):
# Izračunavamo i vraćamo vrijednost Jacobijeve matrice funkcije f u točki x.    
    x1 = x[0]; x2 = x[1];
    Jfx = np.array( [
      [ 6*x1**5 - 10*x1*(x2**2), -10*(x1**2)*x2 ],
      [ -12*(x1**3)*x2,          4*(x2**3) - 3*(x1**4) ]
    ] );

    return Jfx;

x0 = np.array( [1, 2] );

# Iteracije zaustavljamo nakon max 20 koraka 
# ili ako niz postane "stacionaran" uz toleranciju 2.22*10^(-16).
(x, hist) = newton_sustav( f, Jf, x0, 20, 2.22e-16 );

tablica_newton_sustav( hist, f );

 k        x_1              x_2              ||f(x)||  
------------------------------------------------------------
1.000000000000  2.000000000000    147.61097520171052
4.542291950887  1.828103683492   8863.03030110632790
3.775877835532  1.294999404577   2999.01253720699560
3.123284551508  1.027257121946   1034.76915095180448
2.545619794677  1.044378193958    376.16882189724976
1.995630009701  1.565902932840    150.77807034073686
1.697787817020  3.303311264075    116.76086717144095
2.073054677756  2.934926464100     31.42852345303919
2.081716791664  3.178510792324      3.36156552196635
2.088346798860  3.168604854290      0.01691824302317
2.088378995687  3.168732957947      0.00000071691863
2.088378995521  3.168732953137      0.00000000000011
2.088378995521  3.168732953137      0.00000000000004

Vidimo da je Newtonova metoda stala nakon \(12\) iteracija s vrlo točnom aproksimacijom \(x_1 \approx 2.088378995521\), \(x_2 \approx 3.168732953137\); kada uvrstimo ove aproksimacije u sustav onda je desna strana u svakoj od jednadžbi po modulu manja od \(4 \cdot 10^{-14}\).

Napomena

Ako sustav ima mnogo nepoznanica ili je funkcija \(f\) vrlo složena, računanje Jacobijeve matrice može biti računski vrlo zahtjevno. Zbog toga su osmišljene tzv. kvazi-Newtonove metode koje u svakom koraku rade aproksimaciju Jacobijeve matrice. Pokažimo samo neke osnovne ideje jedne ovakve metode. Na primjer, Broydenova metoda umjesto rješavanja sustava

\[ J_f(x^{(k)}) d^{(k)} = -f(x^{(k)}), \]

u svakom koraku korekciju \(d^{(k)}\) računa tako da rješava sustav

(6.6)#\[ Q^{(k)} d^{(k)} = -f(x^{(k)}),\]

gdje je \(Q^{(k)} \approx J_f(x^{(k)})\) aproksimacija Jacobijeve matrice. Kako su elementi Jacobijeve matrice parcijalne derivacije funkcija \(f_i\), možemo ih aproksimirati podijeljenim razlikama. Drugim riječima, matricu \(Q^{(k)}\) možemo odabrati tako da vrijedi

(6.7)#\[ Q^{(k)} (x^{(k)} - x^{(k-1)}) = f(x^{(k)}) - f(x^{(k-1)}).\]

Uočite da ova jednakost ne određuje matricu \(Q^{(k)}\) jedinstveno, nego definira njezino ponašanje za samo jedan vektor. S ciljem da računanje bude što je moguće efikasnije, matrica \(Q^{(k)}\) se određuje tako da se pronađe korekcija ranga jedan matrice \(Q^{(k-1)}\) minimalne norme za koju vrijedi gornja jednakost. To vodi na sljedeću formulu (bez dokaza):

\[ Q^{(k)} = Q^{(k-1)} + \frac{ f(x^{(k)}) - f(x^{(k-1)}) - Q^{(k-1)} d^{(k-1)}}{\|d^{(k-1)}\|^2} (d^{(k-1)})^T. \]

Uvjerite se da ovako definirana matrica \(Q^{(k)}\) zaista zadovoljava (6.7). Dodatno, s ovako definiranom \(Q^{(k)}\) moguće je izbjeći i rješavanje sustava (6.6) uz korištenje Sherman-Morrison-Woodburyjeve formule. Kao početna aproksimacija \(Q^{(0)}\) može se uzeti Jacobijeva matrica funkcije \(f\) u točki \(x^{(0)}\) ili neka njegova aproksimacija.

Metoda fiksne točke#

Kao i kod nelinearnih jednadžbi s jednom nepoznanicom, rješavanje sustava nelinearnih jednadžbi možemo reformulirati kao problem nalaženja fiksne točke \( \alpha\) neke funkcije \(\varphi : \R^n \to \R^n\). Ponovno, nakon odabira pogodne početne iteracije \(x^{(0)}\), promatramo niz \((x^{(k)})_k\) definiran sa

\[ x^{(k+1)} = \varphi(x^{(k)}). \]

Kažemo da je gornji niz iteracija definiran metodom fiksne točke ili metodom jednostavnih iteracija.

Definicija 6.19

Neka je \(\Omega \subseteq \R^n\). Preslikavanje \(\varphi : \Omega \to \R^n\) je kontrakcija ako postoji konstanta \(q < 1\) takva da je

\[ \|\varphi(x) - \varphi(y)\| \leq q \|x - y\|, \]

za sve \(x, y \in \Omega\). Ovdje je \(\|\cdot\|\) neka vektorska norma.

Rezultate dobivene za metodu jednostavnih iteracija sada možemo poopćiti i na sustave. Kontrakcija ima jedinstvenu fiksnu točku ako preslikava zatvoreni skup \(\Omega\) u \(\Omega\).

Teorem 6.20 (Kontrakcije s jedinstvenom fiksnom točkom)

Neka je \(\Omega \subseteq \R^n\) zatvoren skup, te \(\varphi : \Omega \to \R^n\) kontrakcija takva da je \(\varphi(\Omega) \subseteq \Omega\). Tada \(\varphi\) ima jedinstvenu fiksnu točku u \(\Omega\).

Dokaz je posve isti kao i u slučaju \(n=1\); provedite ga za vježbu. Isto vrijedi i za sljedeći teorem, gdje raniji uvjet \(|\varphi'(\alpha)| < 1\) prelazi u \(\|J_{\varphi}(\alpha)\| < 1\).

Teorem 6.21 (Konvergencija metode fiksne točke)

Pretpostavimo da funkcija \(\varphi : \Omega \to \R^n\) ima fiksnu točku \(\alpha\) u interioru od \(\Omega\), te da je \(\varphi\) klase \(C^1\) u nekoj okolini of \(\alpha\). Nadalje, pretpostavimo da za Jacobijevu matricu funkcije \(\varphi\) vrijedi

\[ \|J_{\varphi}(\alpha)\| < 1, \]

gdje je \(\|\cdot\|\) neka (bilo koja!) matrična norma.

Tada postoji okolina \(\hat{\Omega} \subseteq \Omega\) od \(\alpha\) takva da, za svaku \(x^{(0)} \in \hat{\Omega}\), sve iteracije \(x^{(k)}\) generirane metodom fiksne točke leže u \(\Omega\) i konvergiraju ka \(\alpha\).

Primjer 6.22

Promotrimo nelinearnih sustav jednadžbi

\[\begin{split} \begin{align*} x_1^2 + x_2^2 - 1 &= 0, \\ 2x_1 + x_2 - 1 &= 0. \end{align*} \end{split}\]

Ovaj sustav je lako egzaktno riješiti; njegova rješenja su \(\alpha_1 = \mb{c} 0 \\ 1 \me\) i \(\alpha_2 = \mb{c} 4/5 \\ -3/5 \me\). Definirajmo dvije iteracijske funkcije \(\varphi_{1, 2} : \R^2 \to \R^2\) sa

\[\begin{split} \varphi_1(x_1, x_2) = \mb{cc} \frac{1-x_2}{2} \\ \sqrt{1 - x_1^2} \me, \quad \varphi_2(x_1, x_2) = \mb{cc} \frac{1-x_2}{2} \\ -\sqrt{1 - x_1^2} \me. \end{split}\]

Lako se provjeri da je \(\varphi_1 (\alpha_1) = \alpha_1\) i \(\varphi_2(\alpha_2) = \alpha_2\). Pripadne Jacobijeve matrice su

\[\begin{split} J_{\varphi_1}(x_1, x_2) = \mb{cc} 0 & \frac{-1}{2} \\ \frac{-x_1}{\sqrt{1-x_1^2}} & 0 \me, \quad J_{\varphi_2}(x_1, x_2) = \mb{cc} 0 & \frac{-1}{2} \\ \frac{x_1}{\sqrt{1-x_1^2}} & 0 \me, \end{split}\]

pa uvrštavanjem \(\alpha_1\) i \(\alpha_2\) imamo

\[\begin{split} J_{\varphi_1}(\alpha_1) = \mb{cc} 0 & \frac{-1}{2} \\ 0 & 0 \me, \quad J_{\varphi_2}(\alpha_2) = \mb{cc} 0 & \frac{-1}{2} \\ \frac{4}{3} & 0 \me. \end{split}\]

Vidimo da je \(\|J_{\varphi_1}(\alpha_1)\|_1 = \frac{1}{2}\), pa će metoda fiksne točke konvergirati za početnu točku u okolini od \(\alpha_1\).

S druge strane, \(\|J_{\varphi_2}(\alpha_2)\|_1 = \|J_{\varphi_2}(\alpha_2)\|_{\infty} = \|J_{\varphi_2}(\alpha_2)\|_2 = \frac{4}{3} \geq 1\), pa se čini da metoda fiksne točke neće konvergirati ka \(\alpha_2\). Međutim, može se pokazati da je uvjet iz teorema o konvergenciji ekvivalentan zahtjevu

\[ \rho( J_{\varphi}(\alpha) ) < 1, \]

gdje je sa

\[ \rho(A) = \max\{ |\lam| \;:\; \lam \text{ je svojstvena vrijednost od } A \} \]

označen spektralni radijus matrice \(A\). Svojstvene vrijednosti matrice \(J_{\varphi_2}(\alpha_2)\) su \(i \sqrt{2/3}\) i \(-i\sqrt{2/3}\), pa je \(\rho(J_{\varphi_2}(\alpha_2)) = \sqrt{2/3}\) i metoda fiksne točke će konvergirati! Dakle, postoji neka matrična norma takva da je \(\|J_{\varphi_2}(\alpha_2)\| < 1\), no to nije ni \(1\) ni \(2\) ni \(\infty\) norma.

Metoda fiksne točke za linearne sustave

Metodom fiksne točke smo ranije već rješavali sustave linearnih jednadžbi: kako bismo odredili rješenje sustava \(Ax = b\), proučavali smo iterativne metode oblika

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

U Propoziciji 2.26 smo pokazali da takva iterativna metoda konvergira ako je \(\|B\| < 1\) za bilo koju operatorsku matričnu normu. Iskazali smo i Teorem 2.29 koji kaže da to vrijedi ako i samo ako je \(\rho(B) < 1\), gdje je \(\rho(B)\) spektralni radijus matrice \(B\).

Označimo \(\varphi(x) = B x + c\). Promotrimo sada rješenje linearnog sustava kao fiksnu točku funkcije \(\varphi\). Jacobijeva matrica ove funkcije je \(J_\varphi(x) = B\), pa vidimo da je Propozicija 2.26 specijalni slučaj Teorema 6.21. Slično, zaključak o spektralnom radijusu matrice \(J_\varphi(\alpha)\) iz gornjeg primjera poopćava raniji rezultat za sustave linearnih jednadžbi.

Implementirajmo metodu fiksne točke u Pythonu i uvjerimo se da u gornjem primjeru zaista imamo konvergenciju prema nultočki \(\alpha_2\).

def metoda_fiksne_tocke( phi, x0, k_max, eps ):
# Ulaz:
#   - phi = funkcija čiju fiksnu točku tražimo.
#   - x0 = početna iteracija
#   - k_max = broj iteracija metode fiksne točke
# Izlaz:
#   - x = aproksimacija nultočke od f dobivena nakon k_max iteracija metode fiksne točke 
#         ili nakon što niz postane stacionaran.
#   - hist = polje svih izračunatih iteracija (treba samo za tablični prikaz kasnije).

    x = x0; hist = [x0];
    for k in range( 0, k_max ):
        xk = x;
        
        # Iduća iteracija.
        x = phi(xk);

        # Spremi novu iteraciju u hist.
        hist.append( x );

        # Provjeri je li niz postao "dovoljno stacionaran".
        if( np.max( np.abs( x - xk ) / np.abs( xk ) ) < eps ):
            break;

    return (x, hist);

# Primjer: Nađimo nultočku alpha2.
# Ona je fiksna točka funkcije phi2.
def phi2( x ):
    x1 = x[0]; x2 = x[1];

    phi2_x = np.array( 
        [
            (1-x2)/2,
            -np.sqrt( 1 - x1**2 )
        ] );

    return phi2_x;

# Odredimo nultočku alpha2. 
# Početna iteracija = (0.9, 0.9).
# Stajemo nakon max 150 iteracija ili kad je relativna promjena <10^(-10).
x0 = np.array( [0.9, 0.9] )
(x, hist) = metoda_fiksne_tocke( phi2, x0, 150, 1e-10 );

# Potrebno je puno iteracija pa ćemo samo ispisati zadnju.
print( f'Metoda je zaustavljena nakon {len(hist)} iteracija.' );
print( f'Aproksimacija nultočke alpha2: x = ({x[0]:.10f}, {x[1]:.10f})' );

Metoda je zaustavljena nakon 118 iteracija.
Aproksimacija nultočke alpha2: x = (0.8000000000, -0.6000000000)

Na kraju ove cjeline samo spomenimo da se Newtonova metoda kojom rješavamo nelinearni sustav \(f(x) = 0\) ponovno može protumačiti kao metoda fiksne točke za funkciju iteracije

\[ \varphi(x) = x - J_f(x)^{-1} \cdot f(x), \]

na temelju čega se može analizirati konvergencija Newtonove metode i u slučaju nelinearnih sustava.