Show code cell source Hide code cell source

import numpy as np;
import scipy as sc;
import matplotlib.pyplot as plt;

%matplotlib inline

2.2. Uvjetovanost i stabilnost rješavanja sustava linearnih jednadžbi

Uvjetovanost rješavanja sustava \(Ax=b\)

Kao za sve probleme za koje se susrećemo u numeričkoj matematici, možemo pričati o uvjetovanosti. Sjetimo se, ideja uvjetovanosti je za preslikavanje \(F : X \to Y\), koje neki \(x\in X\) šalje u \(y = F(x) \in Y\), odrediti koliko griješimo ukoliko se u domeni malo odmaknemo od \(x\). Ako su \(X\) i \(Y\) vektorski prostori s normama \(\| \cdot \|_X\) i \(\| \cdot \|_Y\), za relativnu uvjetovanost tražimo broj \(\kappa_F(x)\) takav da vrijedi

\[\frac{\| \Delta y \|_Y}{\| y \|_Y} \leq \kappa^{\text{rel}}_F (x) \frac{\| \Delta x \|_X}{\| x \|_X},\]

pri čemu je \(F(x+\Delta x) = y+ \Delta y\).

U rezultatima ovog poglavlja smatramo da je matrična norma \(\|\cdot\|\) operatorska norma koja je inducirana vektorskom normom \(\|\cdot\|\).

Promotrimo problem u kojem za neku regularnu kvadratnu matricu \(A\) reda \(n\) računamo rješenje sustava \(Ax=b\) za razne desne strane \(b\). Radi se o preslikavanju koje vektoru \(b \in \mathbb R^n\) pridružuje \(x = A^{-1}b \in \mathbb R^n\).

Teorem 2.16

Neka je \(A \in \R^{n \times n}\) regularna kvadratna matrica, te \(b\), \(\Delta b\), \(x\) i \(\Delta x\) iz \(\mathbb R^n\) takvi da vrijedi

\[Ax=b \quad \text{i} \quad A(x+\Delta x) = b+\Delta b.\]

Ako \(x\not =0\), tada vrijedi

\[\frac{\| \Delta x \|}{\| x \|} \leq \kappa (A) \frac{\| \Delta b \|}{\| b \|},\]

pri čemu je \(\kappa (A) = \| A \| \cdot\| A^{-1} \|\).

Dokaz.

Oduzmimo dvije jednakosti iz pretpostavke. Dobivamo \(A \Delta x = \Delta b\), odakle množenjem slijeva s \(A^{-1}\) dobivamo

\[ \Delta x = A^{-1} \Delta b.\]

Uvedimo pokratu \(\varepsilon:= \dfrac{ \| \Delta b \|}{ \| b \|}\). Primjenom norme na svaku stranu gornje jednakosti i ograničavanjem dobivamo

\[\| \Delta x \| = \|A^{-1} \Delta b \| \leq \|A^{-1}\| \cdot \| \Delta b \| = \varepsilon \|A^{-1}\| \cdot \| b \| = \varepsilon \|A^{-1}\| \cdot \| Ax \| \leq \varepsilon \|A^{-1}\| \cdot \| A \| \cdot \|x \|.\]

Ovdje smo za obje nejednakosti iskoristili činjenicu da je \(\|\cdot\|\) operatorska matrična norma. Dijeljenjem svake strane s \(\|x \|\) dobivamo tvrdnju.

Broj \(\kappa (A) = \| A \| \cdot\| A^{-1} \|\) zovemo uvjetovanost matrice \(A\). Nije rijetkost da se taj broj pojavljuje u promatranju uvjetovanosti problema koji uključuju matrice, i on mjeri koliko ta matrica ima potencijal da numerički algoritam koji se na njoj provodi bude pogrešno izveden. Što je uvjetovanost matrice veća, to je ograda u gornjem rezultatu dana s većim faktorom što znači da mala perturbacija desne strane sustava \(Ax=b\) može značiti veliku grešku rješenja \(x\).

Uvjetovanost matrice očito ovisi o normi u kojoj se računa. Neovisno o tome, uvjetovanost matrice uvijek je veća ili jednaka od \(1\). Naime \(1 = \| I \| = \| A \cdot A^{-1} \|\leq \| A \| \cdot\| A^{-1} \|\). Jednakost se u \(2\)-normi postiže za matrice oblika \(\lambda U\), gdje je \(\lambda \not =0 \) skalar, te \(U\) unitarna matrica.

Gornji rezultat može biti poopćen i u slučaju kada se i matrica \(A\) perturbira (a ne samo desna strana \(b\) sustava).

Teorem 2.17

Neka je \(A \in \R^{n \times n}\) regularna kvadratna matrica, \(\Delta A \in \R^{n \times n}\) te \(b\), \(\Delta b\), \(x\) i \(\Delta x\) iz \(\mathbb R^n\) takvi da vrijedi

\[Ax=b \quad \text{i} \quad (A+ \Delta A) (x+\Delta x) = b+\Delta b.\]

Neka je \(\varepsilon >0\) takav da je

\[\| \Delta A \| \leq \varepsilon \| A \| \quad \text{i} \quad \| \Delta b \| \leq \varepsilon \| b \|.\]

Pretpostavimo dodatno da je \(x\not =0\), te \(\varepsilon \kappa (A) < 1\). Tada vrijedi

\[\frac{\| \Delta x \|}{ \| x \|} \leq \frac{2\varepsilon\kappa (A)}{1-\varepsilon\kappa (A)}.\]

Dokaz.

Oduzmimo dvije jednakosti iz pretpostavke. Dobivamo

\[A \Delta x = \Delta b - \Delta A (x+\Delta x)\]

odakle množenjem slijeva s \(A^{-1}\) dobivamo

\[\Delta x = A^{-1} \Delta b - A^{-1} \Delta A (x+\Delta x).\]

Primjenom norme i ograničavanjem, dobivamo

\[\|\Delta x\| \leq \| A^{-1}\|\cdot \| \Delta b\| + \|A^{-1}\| \cdot \|\Delta A\| \cdot \| x+\Delta x\|.\]

Za \(\Delta b\) ograničavamo slično kao u prošlom dokazu:

\[\| \Delta b \| \leq \varepsilon \| b \| = \varepsilon \| Ax \| \leq \varepsilon \| A \| \cdot \|x \|.\]

Kad to uvrstimo nazad, dobivamo

\[\begin{split}\begin{align*} \|\Delta x\| & \leq \| A^{-1}\|\cdot \| \Delta b\| + \|A^{-1}\| \cdot \|\Delta A\| \cdot \| x+\Delta x\| \\ & \leq \varepsilon \| A^{-1}\|\cdot \| A \| \cdot \|x \|+ \varepsilon \|A^{-1}\| \cdot \| A\| \cdot \left( \| x \| +\| \Delta x\| \right) \\ & = 2\varepsilon \kappa(A) \|x \|+ \varepsilon \kappa(A) \| \Delta x\|. \end{align*}\end{split}\]

Prebacimo li zadnji sumand na lijevu stranu i podijelimo s \((1-\varepsilon\kappa (A)) \|x \|\) (za što znamo da je pozitivno iz pretpostavke teorema), dobivamo tvrdnju.

Pokažimo još jedan rezultat koji se veže uz uvjetovanost matrice i sustave \(Ax=b\).

Definicija 2.18 (Rezidual)

Neka je \(A \in \R^{n \times n}\) regularna kvadratna matrica, te neka su dani vektori \(\hat{x},b \in \mathbb R^n\). Vektor

\[r(\hat{x}) = A\hat{x} - b\]

zove se rezidual sustava \(Ax=b\) za vektor \(\hat{x}\).

Ako o vektoru \(\hat{x}\) razmišljamo kao o aproksimaciji rješenja \(x\) sustava \(Ax=b\), rezidual je vektor koji mjeri koliko je ta aproksimacija pogrešna. Naravno, rezidual \(r(\hat{x})\) je jednak nulvektoru ako i samo ako je vektor \(\hat{x}\) upravo rješenje \(x\).

Teorem 2.19

Neka je \(A \in \mathbb M_n\) regularna matrica i \(b \in \mathbb R^n\). Neka je \(x\in \mathbb R^n\) rješenje sustava \(Ax=b\), te neka je \(\hat{x}\in \mathbb R^n\). Vrijedi

\[\frac{\| x-\hat{x}\|}{\| x\|} \leq \kappa (A) \frac{\| r(\hat{x})\|}{\| b\|}.\]

Dokaz.

Vrijedi

\[\begin{split}\begin{align*} \frac{\| x-\hat{x}\|}{\| x\|} &=\frac{\| A^{-1}A(x-\hat{x})\|}{\| x\|}\leq \| A^{-1}\|\frac{\| Ax-A\hat{x}\|}{\| x\|}\\ &=\| A^{-1}\|\frac{\| b-A\hat{x}\|}{\| x\|}\frac{\| A\|}{\| A\|} = \kappa (A) \frac{\| r(\hat{x})\|}{\| A\| \| x\|} \\ &\leq \kappa (A) \frac{\| r(\hat{x})\|}{\| Ax\|}= \kappa (A) \frac{\| r(\hat{x})\|}{\| b\|}. \end{align*}\end{split}\]

Ovaj rezultat potvrđuje našu hipotezu: što je rezidual manji, to je aproksimacija rješenja točnija. Ipak, mali rezidual ne mora (gotovo) ništa značiti ako je matrica loše uvjetovana.

Hilbertova matrica

Definicija 2.20 (Hilbertova matrica)

Za prirodan broj \(n\) definiramo Hilbertovu matricu \(H_n\) reda \(n\) na sljedeći način: na mjestu \((i,j)\), \(i,j \in \{ 1,2,\dots,n\}\), nalazi se broj

\[\frac{1}{i+j-1}.\]

Recimo, za \(n=5\) imamo

\[\begin{split}H_5 = \begin{bmatrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{bmatrix}.\end{split}\]

Ove matrice nisu „umjetne”, pojavljuju se u aproksimacijama funkcija polinomom, kao što ćemo vidjeti kasnije. Matrica ima naoko lijepa svojstva: simetrična je i pozitivno definitna. No, vidjet ćemo da je s njima teško raditi. Pokažimo na primjeru problem rješavanja sustava \(Ax=b\) Hilbertovom matricom.

Izračunajmo na računalu rješenje sustava Hilbertovom matricom. Rješavamo za dimenziju \(n=1,2,\dots,20\). Matrica sustava je \(H_n\), a vektor desne strane \(b\) namještamo tako da je egzaktno rješenje sustave jednako \(x = (1,\dots,1)^T\). Zatim, sustav rješavamo pomoću Pythonove funkcije solve te tako dobivamo vektor \(\tilde{x}\) kao aproksimaciju za \(x\). Na kraju, ispisujemo relativnu grešku \(\frac{\|x - \tilde{x}\|}{\|x\|}\).

Show code cell source Hide code cell source

def hilbert_matrix(n):
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            H[i, j] = 1 / (i + j + 1)

    return H

for n in range(1, 21):
    x = np.ones(n)
    H = hilbert_matrix(n)
    b = H @ x
    xtilde = np.linalg.solve(H, b)
    relativna_greska = np.linalg.norm(x-xtilde) / np.linalg.norm(x)
    print(f'Relativna greška za n={n}:  \t{relativna_greska:.2e}')

Relativna greška za n=1:  	0.00e+00
Relativna greška za n=2:  	5.66e-16
Relativna greška za n=3:  	8.36e-15
Relativna greška za n=4:  	1.23e-13
Relativna greška za n=5:  	1.57e-11
Relativna greška za n=6:  	2.58e-10
Relativna greška za n=7:  	3.52e-09
Relativna greška za n=8:  	1.21e-07
Relativna greška za n=9:  	8.14e-06
Relativna greška za n=10:  	4.17e-04
Relativna greška za n=11:  	1.19e-02
Relativna greška za n=12:  	1.04e-03
Relativna greška za n=13:  	8.28e-02
Relativna greška za n=14:  	4.69e+00
Relativna greška za n=15:  	8.73e+00
Relativna greška za n=16:  	3.87e+00
Relativna greška za n=17:  	3.51e+00
Relativna greška za n=18:  	7.61e+00
Relativna greška za n=19:  	8.47e+00
Relativna greška za n=20:  	2.03e+01

Uočavamo da rješavanjem sustava s Hilbertovom matricom na računalu dobivamo ogromne relativne greške: već za \(n=14\) greška je veća od \(100\%\)!

Razlog za to je jako velika uvjetovanost Hilbertove matrice. Ispisujemo uvjetovanosti u \(2\)-normi:

Show code cell source Hide code cell source

for n in range (1, 21):
    H = hilbert_matrix(n)
    # Funkcija cond iz numpy.linalg računa uvjetovanost u 2-normi.
    cond = np.linalg.cond(H)
    print(f'Uvjetovanost Hilbertove matrice reda {n}: \t{cond:.2e}')

Uvjetovanost Hilbertove matrice reda 1: 	1.00e+00
Uvjetovanost Hilbertove matrice reda 2: 	1.93e+01
Uvjetovanost Hilbertove matrice reda 3: 	5.24e+02
Uvjetovanost Hilbertove matrice reda 4: 	1.55e+04
Uvjetovanost Hilbertove matrice reda 5: 	4.77e+05
Uvjetovanost Hilbertove matrice reda 6: 	1.50e+07
Uvjetovanost Hilbertove matrice reda 7: 	4.75e+08
Uvjetovanost Hilbertove matrice reda 8: 	1.53e+10
Uvjetovanost Hilbertove matrice reda 9: 	4.93e+11
Uvjetovanost Hilbertove matrice reda 10: 	1.60e+13
Uvjetovanost Hilbertove matrice reda 11: 	5.23e+14
Uvjetovanost Hilbertove matrice reda 12: 	1.62e+16
Uvjetovanost Hilbertove matrice reda 13: 	4.95e+17
Uvjetovanost Hilbertove matrice reda 14: 	2.57e+17
Uvjetovanost Hilbertove matrice reda 15: 	3.19e+17
Uvjetovanost Hilbertove matrice reda 16: 	3.00e+17
Uvjetovanost Hilbertove matrice reda 17: 	6.98e+17
Uvjetovanost Hilbertove matrice reda 18: 	2.52e+18
Uvjetovanost Hilbertove matrice reda 19: 	3.16e+18
Uvjetovanost Hilbertove matrice reda 20: 	3.31e+18

Razlog velike uvjetovanosti matrice \(H_n\) su veliki elementi u inverzu. Pogledajmo \(H_5^{-1}\):

Show code cell source Hide code cell source

# Funkcija inv iz numpy.linalg vraća inverz matrice.
# Općenito je LOŠA PRAKSA u numeričkim algoritmima računati inverze matrica !!!
# Ovdje nam to služi samo za ilustraciju.
print(np.linalg.inv(hilbert_matrix(5)))

[[ 2.500e+01 -3.000e+02  1.050e+03 -1.400e+03  6.300e+02]
 [-3.000e+02  4.800e+03 -1.890e+04  2.688e+04 -1.260e+04]
 [ 1.050e+03 -1.890e+04  7.938e+04 -1.176e+05  5.670e+04]
 [-1.400e+03  2.688e+04 -1.176e+05  1.792e+05 -8.820e+04]
 [ 6.300e+02 -1.260e+04  5.670e+04 -8.820e+04  4.410e+04]]

Općenito za uvjetovanost Hilbertove matrice reda \(n\) u normi \(\|\cdot \|_2\) vrijedi

\[\kappa_2(H_n) \approx \frac{\displaystyle \big( \sqrt{2} + 1 \big)^{4n + 4}} {\displaystyle 2^{15/4} \, \sqrt{\pi n}},\]

kad \(x \to \infty\), pa vidimo da uvjetovanost raste eksponencijalno s dimenzijom matrice.

Stabilnost LU faktorizacije

Na točnost rješavanja sustava \(Ax=b\), pored uvjetovanosti tog problema, odnosno uvjetovanosti matrice \(A\), utječe i stabilnost numeričkog algoritma kojim rješavamo taj sustav. Vidjeli smo da nije isto rješavati sustav LU faktorizacijom i LU faktorizacijom s parcijalnim pivotiranjem. To se može pokazati i nekim metrikama.

Prilikom rješavanja sustava LU faktorizacijom (bez i s pivotiranjem) imamo tri koraka: provođenje same LU faktorizacije, rješavanje jednog donjetrokutastog sustava i rješavanje gornjetrokutastog sustava. Pokazuje se da je rješavanje trokutastih sustava supstitucijom unaprijed ili unatrag stabilno, tako da problem u stabilnosti ovog algoritma je sama LU faktorizacija, odnosno koliko su ti faktori \(L\) i \(U\) točno izračunati.

Sljedeći rezultat nećemo dokazivati (vidi Teorem 9.3 u [Hig02]), a on govori o grešci unatrag algoritma za računanje LU-faktorizacije u uvjetima aritmetike konačne preciznosti na računalu.

Teorem 2.21

U aritmetici računala računamo LU faktorizaciju zadane kvadratne regularne matrice \(A\) reda \(n\). Pretpostavimo da je algoritam uspješno završio bez pojave brojeva koji nisu prikazivi u računalu (jer su preveliki ili premali) i bez pokušaja djeljenja s nulom.

Pretpostavimo da su umjesto točnih faktora \(L\) i \(U\) takvih da je \(A=LU\) izračunati faktori \(\hat{L}\) i \(\hat{U}\). Neka je \(\Delta A\) matrica takva da je \(\hat{L} \hat{U} = A + \Delta A\). Tada je

(2.8)\[|\Delta A| \leq \gamma_n |\hat{L}| |\hat{U}|,\]

pri čemu je \(\displaystyle \gamma_n := \frac{n u}{1 - n u}\), \(u\) jedinična greška zaokruživanja, te \(|X|\) matrica dobivena od \(X\) tako da je svakom elementu primijenjena apsolutna vrijednost. Gornja nejednakost vrijedi za svaki element matrica s lijeve i desne strane.

Kada ovakvu \(\Delta A\) uvrstimo u rezultat Teorema 2.17, možemo dobiti ocjenu za relativnu grešku za \(x\). Dakle, točnost rješavanja sustava \(Ax=b\) ovisi o broju \(|\hat{L}| |\hat{U}|\) za koji želimo da je što manji.

U idealnom slučaju kada je \(|\hat{L}| |\hat{U}| = |\hat L \hat U|\), prema gornjem teoremu bismo imali

\[ |\hat{L}| |\hat{U}| = |\hat{L} \hat{U}| = |A + \Delta A| \leq |A| + |\Delta A| \leq |A| + \gamma_n |\hat{L}| |\hat{U}|,\]

pa kad prebacimo članove na lijevu stranu i ponovno iskoristimo gornji teorem

\[|\Delta A| \leq \gamma_n |\hat{L}| |\hat{U}| \leq \frac{\gamma_n}{1 - \gamma_n} |A|.\]

Problem je što općenito ne mora vrijediti \(|\hat{L}| |\hat{U}| = |\hat L \hat U|\). Ova jednakost vrijedi u nekim važnim situacijama (recimo kod splajn interpolacije, što ćemo vidjeti kasnije), kada kod oba faktora \(L\) i \(U\) imamo elemente koji su svi nenegativni. No u općenitom slučaju, taj odnos brojeva \(|\hat{L}| |\hat{U}| \) i \( |\hat L \hat U| = |A|\) može biti velik, pogotovo kod \(LU\) faktorizacije bez pivotiranja, što rezultira velikom greškom. Da bi se kontrolirala veličina te greške, u analizi se često računa omjer

\[\dfrac{\| |\hat{L}| |\hat{U}| \| }{\| A \|}.\]

Što je taj omjer manji, to očekujemo veću točnost rješenja.

Nadalje, u ogradi (2.8) možemo računati \(\| \cdot \|_{\infty}\) normu matrica koje se pojavljuju u ocjeni. Tako broj \(\| A \|_{\infty}\) ocjenjujemo brojevima \(\| \hat{L} \|_{\infty}\) i \(\| \hat{U} \|_{\infty}\). Kod LU faktorizacije s parcijalnim pivotiranjem znamo da su svi elementi matrice \(L\) po apsolutnoj vrijednosti manji od \(1\), pa vrijedi \(\| \hat{L} \|_{\infty} \leq n\). Preostalo je ocijeniti veličinu elemenata matrice \(\hat{U}\). To motivira uvođenje broja

\[ \rho_n (A) = \frac{\max_{i, j, k} |a_{i,j}^{(k)}|}{\max_{i, j} |a_{i,j}|},\]

gdje su \(a_{i,j}^{(k)}\) elementi matrice \(A^{(k)}\) u \(k\)-tom koraku LU faktorizacije. Taj broj naziva se pivotni rast. Kako je \(A^{(n)} = U\), vrijedi \(|u_{i,j}| \leq \rho_n (A) \max_{i, j} |a_{i,j}| \). Ponovno, cilj je da je taj broj što manji. Pivotni rast, osim o matrici \(A\), ovisi i o metodi nalaženja faktora u LU faktorizaciji (tj. provodimo li pivotiranje ili ne).

Odnos pivotnog rasta i povratne greške iskazan je u sljedećem teoremu (navodimo bez dokaza; za dokaz vidjeti [Hig02]).

Teorem 2.22 (Wilkinson)

Neka je \(A\) kvadratna regularna matrica reda \(n\) i neka je \(\hat{x}\) izračunato rješenje sustava \(Ax=b\) Gaussovim eliminacijama s parcijalnim pivotiranjem u aritmetici računala. Neka je \(\Delta A\) matrica takva da je \((A + \Delta A) \hat{x} = b\). Tada vrijedi:

(2.9)\[\| \Delta A \|_\infty \leq n^2 \gamma_{3n} \rho_n^{(p)} \| A \|_\infty,\]

pri čemu je \(\displaystyle \gamma_n := \frac{n u}{1 - n u}\), a \(\rho_n^{(p)}\) pivotni rast u prilikom parcijalnog pivotiranja matrice \(A\).

Primjer 2.23

Promotrimo ponovno Primjer 2.2. Matrica sustava je oblika

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

za neki mali \(\varepsilon>0\). Uvjetovanost ove matrice u \(\infty\)-normi je \(\kappa_\infty (A) = \dfrac{4}{1-\varepsilon}\), što nije veliki broj pa bi se sustav s ovom matricom trebao točno riješiti.

LU faktorizacija bez pivotiranja iznosi

\[\begin{split}\begin{bmatrix} \varepsilon & 1 \\ 1 & 1 \end{bmatrix} = A = LU = \begin{bmatrix} 1 & 0 \\ \dfrac{1}{\varepsilon} & 1 \end{bmatrix} \cdot \begin{bmatrix} \varepsilon & 1 \\ 0 & 1-\dfrac{1}{\varepsilon} \end{bmatrix}\end{split}\]

Umnožak \(|L| \cdot |U|\) iznosi

\[\begin{split} |L| |U| = \begin{bmatrix} 1 & 0 \\ \dfrac{1}{\varepsilon} & 1 \end{bmatrix} \cdot \begin{bmatrix} \varepsilon & 1 \\ 0 & \dfrac{1}{\varepsilon}-1 \end{bmatrix} = \begin{bmatrix} \varepsilon & 1 \\ 1 & \dfrac{2}{\varepsilon}-1 \end{bmatrix}.\end{split}\]

Što je \(\varepsilon\) manji, to je element u donjem desnom kutu ove matrice veći, te prema Teoremu 2.21 kvari točnost LU faktorizacije u računalu. Primijetimo da smo ovaj račun provodili za točne faktore \(L\) i \(U\), a ne \(\hat{L}\) i \(\hat{U}\) koji bi bili dobiveni aritmetikom računala. Pivotni rast u ovom slučaju je \(\rho(A) = \dfrac{1}{\varepsilon}-1\), što je također nekontrolirano veliko za mali \(\varepsilon\).

Promotrimo sada istu tu matricu kada bismo ju rješavali s parcijalnim pivotiranjem. Potrebno je promatrati matricu \(PA\), gdje \(P = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}\) mijenja retke. LU faktorizacija s parcijalnim pivotiranjem iznosi

\[\begin{split}\begin{bmatrix} 1 & 1 \\ \varepsilon & 1 \end{bmatrix} = PA = LU = \begin{bmatrix} 1 & 0 \\ \varepsilon & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 1 \\ 0 & 1-\varepsilon \end{bmatrix}\end{split}\]

Matrice \(|L|\) i \(|U|\) upravo su jednake \(L\) i \(U\), i svi elementi tih matrica i njihovog umnoška su manji ili jednaki \(1\) po apsolutnoj vrijednosti, što sugerira da je LU faktorizacija u ovom slučaju stabilna. To potvrđuje i pivotni rast koji iznosi \(1\).

Dok ovaj primjer pokazuje da u slučaju kada ne provodimo pivotiranje pivotni rast može biti nekontrolirano velik, u slučaju parcijalnog pivotiranja postoji ograda \(\rho_n (A) \leq 2^{n-1}\) koja ovisi o dimenziji matrice \(A\). J. H. Wilkinson je pokazao da se ta ograda može postići za matrice oblika

\[\begin{split}\begin{bmatrix} \ \ 1 & & & & 1 \\ -1 & \ \ 1 & & & 1 \\ -1 & -1 & \ddots & & 1 \\ -1 & -1 & \ddots & \ \ 1 & 1 \\ -1 & -1 & -1 & -1 & 1 \end{bmatrix}\end{split}\]

Ipak, u praksi se pokazuje da su ovakvi primjeri rijetki. Pivotni rast kod parcijalnog pivotiranja je u praksi mnogo manji, pa je i rješavanje sustava LU faktorizacijom s parcijalnim pivotiranjem stabilno.

Literatura

[Hig02] (1,2)

Nicholas Higham. Accuracy and Stability of Numerical Algorithms. SIAM, 2002.