Metodi iterativi - Jacobi, Gauss-Seidel e convergenza

I metodi diretti (Fattorizzazione LU, pivoting e CholeskyRisolvere $A\mathbf x=\mathbf b$ con Gauss equivale a fattorizzare $A=LU$ ($L$ triangolare inferiore con $1$ sulla diagonale e i moltiplicatori $l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ sotto la diagonale, $U$ la matrice a scala finale) e poi risolvere $L\mathbf y=\mathbf b$ (in avanti) e $U\mathbf x=\mathbf y$ (all'indietro). Esiste (ed è unica) se i minori principali di testa $\det A_k$, $k<n$, sono $\ne0$. Costo $\frac23n^3$ flop per la fattorizzazione, $2n^2$ per ogni termine noto; $\det A=\prod u_{kk}$. Con un pivot piccolo l'algoritmo è instabile: il pivoting parziale (scambio di righe che porta in posizione di pivot l'elemento di modulo massimo della colonna) dà $PA=LU$ con $|l_{ik}|\le1$ ed esiste per ogni $A$ invertibile. Se $A$ è simmetrica definita positiva esiste $A=LL^T$ (Cholesky), con $l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}$ e $l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}$; costo $\frac13n^3$, senza pivoting.Fattorizzazione LU, pivoting e Cholesky →) costano $\frac23n^3$ operazioni e, soprattutto, riempiono la matrice: i fattori $L$ e $U$ di una matrice sparsa (pochi elementi non nulli per riga, come in un reticolo di elementi finiti o in una discretizzazione 2D o 3D) hanno molti più elementi non nulli di $A$. Con milioni di incognite memoria e tempo sono proibitivi. I metodi iterativi evitano la fattorizzazione: partono da una stima iniziale e la correggono ripetutamente, usando in ogni iterazione solo prodotti matrice-vettore, che costano quanto il numero di elementi non nulli di $A$. Se bastano poche iterazioni, vincono.

Prerequisiti: Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari → (norme, raggio spettrale), Metodo di bisezione e iterazioni di punto fissoBisezione: se f è continua e f(a)f(b) < 0 si dimezza l'intervallo tenendo la metà con il cambio di segno; dopo n passi |c_n − α| ≤ (b−a)/2^(n+1), cioè un bit di precisione per iterazione (convergenza lineare, sempre garantita). Punto fisso x = g(x): se g manda [a,b] in sé con |g'| ≤ L < 1 (contrazione) il punto fisso è unico, l'iterazione x_{k+1} = g(x_k) converge da ogni x_0 e |x_k − α| ≤ L^k/(1−L)·|x_1−x_0|; l'errore si moltiplica per circa g'(α) a ogni passo. Una stessa equazione ha molte riscritture x = g(x), alcune convergenti e altre no.Metodo di bisezione e iterazioni di punto fisso → (la forma $\mathbf x=B\mathbf x+\mathbf c$ è un'iterazione di punto fisso).

1. Splitting e iterazione di punto fisso

Si scrive $A=M-N$ con $M$ invertibile e facile da invertire (diagonale o triangolare). Allora $$A\mathbf x=\mathbf b\iff M\mathbf x=N\mathbf x+\mathbf b\iff\mathbf x=\underbrace{M^{-1}N}_{B}\,\mathbf x+\underbrace{M^{-1}\mathbf b}_{\mathbf c}$$ e si definisce $\mathbf x^{(k+1)}=B\mathbf x^{(k)}+\mathbf c$, cioè si risolve a ogni passo il sistema facile $M\mathbf x^{(k+1)}=N\mathbf x^{(k)}+\mathbf b$. La soluzione esatta $\mathbf x$ è un punto fisso. Sottraendo $\mathbf x=B\mathbf x+\mathbf c$ da $\mathbf x^{(k+1)}=B\mathbf x^{(k)}+\mathbf c$ si ottiene la regola dell'errore $\mathbf e^{(k)}=\mathbf x-\mathbf x^{(k)}$: $$\mathbf e^{(k+1)}=B\,\mathbf e^{(k)}\quad\Longrightarrow\quad\mathbf e^{(k)}=B^k\mathbf e^{(0)} .$$

Teorema (convergenza, D). L'iterazione converge per ogni $\mathbf x^{(0)}$ se e solo se $\rho(B)<1$. Basta che $\|B\|<1$ per una norma indotta.

Dimostrazione. ($\Leftarrow$ nel caso $\|B\|<1$) $\|\mathbf e^{(k)}\|\le\|B\|^k\|\mathbf e^{(0)}\|\to0$. Nel caso generale $\rho(B)<1$ esiste una norma indotta con $\|B\|<\rho(B)+\varepsilon<1$ (con la forma di Jordan), e si riconduce al caso precedente. ($\Rightarrow$) Se $\rho(B)\ge1$ sia $\lambda$ un autovalore con $|\lambda|\ge1$ e $\mathbf v$ un autovettore: scegliendo $\mathbf x^{(0)}=\mathbf x-\mathbf v$ (per autovalori complessi si lavora nel complesso, o con parte reale e immaginaria) è $\mathbf e^{(k)}=\lambda^k\mathbf v\not\to0$. $\blacksquare$

La convergenza dipende dunque solo da $B$, non da $\mathbf b$ né dalla stima iniziale. Se $\rho(B)\ge1$ nessuna stima iniziale (tranne quelle fortunate) fa convergere il metodo.

2. Jacobi e Gauss-Seidel

Si decompone $A=D+L+U$, con $D$ diagonale di $A$, $L$ la parte strettamente inferiore, $U$ la parte strettamente superiore (non sono i fattori della LU). Si supponga $a_{ii}\ne0$.

Metodo di Jacobi ($M=D$, $N=-(L+U)$): $\mathbf x^{(k+1)}=-D^{-1}(L+U)\mathbf x^{(k)}+D^{-1}\mathbf b$, cioè, per ogni $i$, $$x_i^{(k+1)}=\frac1{a_{ii}}\Big(b_i-\sum_{j\ne i}a_{ij}\,x_j^{(k)}\Big).$$ Metodo di Gauss-Seidel ($M=D+L$, $N=-U$): $(D+L)\mathbf x^{(k+1)}=-U\mathbf x^{(k)}+\mathbf b$, cioè $$x_i^{(k+1)}=\frac1{a_{ii}}\Big(b_i-\sum_{j<i}a_{ij}\,x_j^{(k+1)}-\sum_{j>i}a_{ij}\,x_j^{(k)}\Big).$$

Interpretazione: l'equazione $i$ viene risolta rispetto all'incognita $x_i$, tenendo fissi gli altri valori. In Jacobi si usano tutti i valori della iterazione precedente (aggiornamento «in parallelo»: si può calcolare ogni $x_i$ indipendentemente, ma servono due vettori). In Gauss-Seidel si usano subito i valori appena calcolati ($x_j^{(k+1)}$ con $j<i$): in genere converge più in fretta e basta un solo vettore, ma l'ordine delle equazioni conta e il metodo è sequenziale. Le matrici di iterazione sono $$B_J=-D^{-1}(L+U),\qquad B_{GS}=-(D+L)^{-1}U .$$

Esempio (a mano). $A=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}$, $\mathbf b=(2,4,10)^T$, soluzione $\mathbf x=(1,2,3)^T$, partenza $\mathbf x^{(0)}=\mathbf 0$. Le equazioni risolte per l'incognita diagonale sono $x_1=\frac{2+x_2}4$, $x_2=\frac{4+x_1+x_3}4$, $x_3=\frac{10+x_2}4$.

Jacobi.

• $\mathbf x^{(1)}$: $x_1=\frac{2+0}4=0{,}5$; $x_2=\frac{4+0+0}4=1$; $x_3=\frac{10+0}4=2{,}5$.
• $\mathbf x^{(2)}$ (usando $\mathbf x^{(1)}$): $x_1=\frac{2+1}4=0{,}75$; $x_2=\frac{4+0{,}5+2{,}5}4=1{,}75$; $x_3=\frac{10+1}4=2{,}75$.
• $\mathbf x^{(3)}$: $x_1=\frac{2+1{,}75}4=0{,}9375$; $x_2=\frac{4+0{,}75+2{,}75}4=1{,}875$; $x_3=\frac{10+1{,}75}4=2{,}9375$.

Gauss-Seidel (ogni componente usa quelle già aggiornate).

• $\mathbf x^{(1)}$: $x_1=\frac{2+0}4=0{,}5$; $x_2=\frac{4+0{,}5+0}4=1{,}125$ (usa il nuovo $x_1$); $x_3=\frac{10+1{,}125}4=2{,}78125$.
• $\mathbf x^{(2)}$: $x_1=\frac{2+1{,}125}4=0{,}78125$; $x_2=\frac{4+0{,}78125+2{,}78125}4=1{,}890625$; $x_3=\frac{10+1{,}890625}4=2{,}972656$.
• $\mathbf x^{(3)}$: $x_1=\frac{2+1{,}890625}4=0{,}972656$; $x_2=\frac{4+0{,}972656+2{,}972656}4=1{,}986328$; $x_3=\frac{10+1{,}986328}4=2{,}996582$.

Dopo tre iterazioni l'errore (in norma $\infty$) è $2{,}7\cdot10^{-2}$ per Gauss-Seidel e $0{,}125$ per Jacobi.

3. Quando convergono

Teorema (dominanza diagonale). Se $A$ è a dominanza diagonale stretta per righe, $|a_{ii}|>\sum_{j\ne i}|a_{ij}|$ per ogni $i$, Jacobi e Gauss-Seidel convergono.

Per Jacobi la dimostrazione è immediata: $B_J=-D^{-1}(L+U)$ ha elementi $-\frac{a_{ij}}{a_{ii}}$ ($j\ne i$) e diagonale nulla, quindi $\|B_J\|_\infty=\max_i\sum_{j\ne i}\frac{|a_{ij}|}{|a_{ii}|}<1$. (Per Gauss-Seidel si dimostra che $\|B_{GS}\|_\infty\le\|B_J\|_\infty<1$.) Nell'esempio precedente $\|B_J\|_\infty=\frac24=\frac12<1$.

Teorema (SPD). Se $A$ è simmetrica definita positiva, Gauss-Seidel converge (Jacobi in generale no: richiede anche $2D-A$ definita positiva). Tridiagonali. Se $A$ è tridiagonale con diagonale non nulla, $\rho(B_{GS})=\rho(B_J)^2$ e Jacobi e Gauss-Seidel convergono o divergono insieme.

Per l'esempio, $B_J$ ha autovalori $\pm\frac{\sqrt2}4$ e $0$ (la matrice tridiagonale con zeri sulla diagonale e $\frac14$ fuori), quindi $\rho(B_J)=0{,}3536$ e $\rho(B_{GS})=0{,}3536^2=0{,}125$: Gauss-Seidel è due volte più veloce in numero di cifre guadagnate per iterazione.

La dominanza è solo sufficiente, e dipende dall'ordine delle equazioni. $A=\begin{pmatrix}1&2\\3&1\end{pmatrix}$ non è dominante: $B_J=\begin{pmatrix}0&-2\\-3&0\end{pmatrix}$ ha autovalori $\pm\sqrt6$, $\rho=2{,}45>1$ e Jacobi diverge per ogni $\mathbf x^{(0)}\ne\mathbf x$. Scambiando le equazioni, $\begin{pmatrix}3&1\\1&2\end{pmatrix}$ è dominante: $\rho(B_J)=\sqrt{\frac16}=0{,}408<1$. In pratica si riordinano le equazioni (o si fa il pivoting) per portare il termine grande sulla diagonale.

Nessun metodo domina l'altro. Esistono matrici per cui Jacobi converge e Gauss-Seidel no (e viceversa): $\begin{pmatrix}1&2&-2\\1&1&1\\2&2&1\end{pmatrix}$ ha $B_J$ nilpotente ($\rho=0$, Jacobi converge in 3 passi) e $B_{GS}$ con autovalori $0,2,2$ ($\rho=2$, diverge); $\begin{pmatrix}2&-1&1\\2&2&2\\-1&-1&2\end{pmatrix}$ ha $\rho(B_J)=\frac{\sqrt5}2=1{,}118>1$ (diverge) e $\rho(B_{GS})=\frac12$ (converge).

4. Stime dell'errore e arresto

A priori. Da $\mathbf e^{(k)}=B^k\mathbf e^{(0)}$: $\|\mathbf e^{(k)}\|\le\|B\|^k\|\mathbf e^{(0)}\|$. Per ottenere $\|\mathbf e^{(k)}\|\le\varepsilon\|\mathbf e^{(0)}\|$ basta $k\ge\dfrac{\ln\varepsilon}{\ln\|B\|}$.

A posteriori (dall'incremento). Siccome $\mathbf x^{(k+1)}-\mathbf x^{(k)}=B\big(\mathbf x^{(k)}-\mathbf x^{(k-1)}\big)$, gli incrementi decrescono geometricamente; sommando la serie geometrica degli incrementi successivi $$\|\mathbf x-\mathbf x^{(k)}\|\le\frac{\|B\|}{1-\|B\|}\,\|\mathbf x^{(k)}-\mathbf x^{(k-1)}\| .$$ Quindi si ferma l'iterazione quando l'incremento è abbastanza piccolo: con $\|B\|$ vicino a $1$ l'incremento piccolo non garantisce un errore piccolo (il fattore $\frac{\|B\|}{1-\|B\|}$ è grande). L'altro criterio è il residuo $\|\mathbf b-A\mathbf x^{(k)}\|\le\text{tol}\,\|\mathbf b\|$, che garantisce un errore relativo $\le\kappa(A)\,\text{tol}$ (Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari →).

Velocità asintotica. L'errore si riduce di circa un fattore $\rho(B)$ a ogni passo (rapporti degli errori che tendono a $\rho$, o oscillano attorno, quando ci sono autovalori $\pm\rho$). Per guadagnare $d$ cifre decimali servono $$k\approx\frac{d}{-\log_{10}\rho(B)}$$ iterazioni. Nell'esempio: Jacobi $\rho=0{,}354$, $-\log_{10}\rho=0{,}45$: 2,2 iterazioni per cifra; Gauss-Seidel $\rho=0{,}125$, $0{,}90$: 1,1 iterazioni per cifra.

Esempio meccanico: barra con sorgente di calore. La temperatura $T$ di una barra di lunghezza $L$ con generazione uniforme e estremi a 0\,^\circC soddisfa $-T''=s$. Con $n$ nodi interni e passo $h=\frac L{n+1}$ le differenze finite danno $2T_i-T_{i-1}-T_{i+1}=h^2s$, cioè $K_n\mathbf T=h^2s\,\mathbf 1$ con $K_n=\operatorname{tridiag}(-1,2,-1)$. Qui $B_J=I-\frac12K_n$ ha autovalori $\cos\frac{j\pi}{n+1}$, quindi $$\rho(B_J)=\cos\frac\pi{n+1},\qquad\rho(B_{GS})=\cos^2\frac\pi{n+1}.$$ Per $n=3$: $0{,}707$ e $0{,}5$. Per $n=10$: $0{,}9595$ e $0{,}9206$, cioè $334$ e $167$ iterazioni per guadagnare 6 cifre. Per $n=100$: $0{,}99952$ e $0{,}99903$, 28,554 e 14,277 iterazioni. Più la griglia è fine, più $\rho\to1$ e più il metodo rallenta (legato al cattivo condizionamento $\kappa\approx\frac{4(n+1)^2}{\pi^2}$ della stessa matrice); i metodi moderni (gradiente coniugato, multigriglia) rimediano. Per questa matrice tridiagonale, tuttavia, Cholesky costa $O(n)$: un metodo iterativo conviene per matrici sparse non a banda stretta (griglie 2D e 3D).

5. Implementazione

import numpy as np

def iterazioni(A, b, metodo, x0, nmax):
    """Jacobi o Gauss-Seidel: restituisce la lista delle approssimazioni x^0, x^1, ..."""
    D, L, U = np.diag(np.diag(A)), np.tril(A, -1), np.triu(A, 1)
    M = D if metodo == "jacobi" else D + L           # A = M - N
    N = M - A
    xs = [np.array(x0, dtype=float)]
    for _ in range(nmax):
        xs.append(np.linalg.solve(M, N @ xs[-1] + b))   # M x_new = N x_old + b
    return xs

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]], dtype=float)
b = np.array([2, 4, 10], dtype=float)
x = np.linalg.solve(A, b)                              # soluzione esatta (1, 2, 3)

D, L, U = np.diag(np.diag(A)), np.tril(A, -1), np.triu(A, 1)
BJ = -np.linalg.solve(D, L + U)
BGS = -np.linalg.solve(D + L, U)
rho = lambda B: max(abs(np.linalg.eigvals(B)))
print(f"rho(Jacobi) = {rho(BJ):.4f}   rho(Gauss-Seidel) = {rho(BGS):.4f}   (cos(pi/4)/2 = {np.cos(np.pi/4)/2:.4f}, quadrato: {(np.sqrt(2)/4)**2:.4f})")

print(" k   Jacobi: errore    rapp.  |  Gauss-Seidel: errore   rapp.")
ej, eg = None, None
for k, (xj, xg) in enumerate(zip(iterazioni(A, b, "jacobi", [0, 0, 0], 8), iterazioni(A, b, "gs", [0, 0, 0], 8))):
    nj, ng = np.linalg.norm(x - xj, np.inf), np.linalg.norm(x - xg, np.inf)
    rj = f"{nj / ej:5.3f}" if ej else "    -"
    rg = f"{ng / eg:5.3f}" if eg else "    -"
    print(f"{k:2d}   {nj:.3e}      {rj}  |  {ng:.3e}        {rg}")
    ej, eg = nj, ng

Output:

rho(Jacobi) = 0.3536   rho(Gauss-Seidel) = 0.1250   (cos(pi/4)/2 = 0.3536, quadrato: 0.1250)
 k   Jacobi: errore    rapp.  |  Gauss-Seidel: errore   rapp.
 0   3.000e+00          -  |  3.000e+00            -
 1   1.000e+00      0.333  |  8.750e-01        0.292
 2   2.500e-01      0.250  |  2.188e-01        0.250
 3   1.250e-01      0.500  |  2.734e-02        0.125
 4   3.125e-02      0.250  |  3.418e-03        0.125
 5   1.562e-02      0.500  |  4.272e-04        0.125
 6   3.906e-03      0.250  |  5.341e-05        0.125
 7   1.953e-03      0.500  |  6.676e-06        0.125
 8   4.883e-04      0.250  |  8.345e-07        0.125

Il rapporto di Gauss-Seidel si assesta su $\rho=0{,}125$. Per Jacobi alterna $0{,}25$ e $0{,}5$ (la media geometrica è $\sqrt{0{,}125}=0{,}354=\rho$): gli autovalori $\pm\rho$ di $B_J$ fanno oscillare l'errore.

In MATLAB non ci sono funzioni predefinite per Jacobi e Gauss-Seidel: si scrivono con un ciclo (x = D\(b-(L+U)*x)); per sistemi sparsi grandi si usano A\b (con matrici sparse), pcg (gradiente coniugato, SPD) o gmres.

Errori da evitare

• Confondere $L$ e $U$ della decomposizione $A=D+L+U$ con i fattori della LU: qui sono le parti strettamente triangolari di $A$.
• Usare i valori nuovi in Jacobi (diventa Gauss-Seidel) o i vecchi in Gauss-Seidel.
• Confondere la dominanza per righe con quella per colonne: la prima corrisponde a $\|B_J\|_\infty<1$, la seconda a $\|B_J\|_1<1$; per Jacobi basta una delle due.
• Fidarsi dell'incremento piccolo quando $\|B\|\approx1$: l'errore può essere $\frac{\|B\|}{1-\|B\|}$ volte più grande.
• Dimenticare di riordinare le equazioni quando la diagonale ha elementi piccoli.
• Credere che la dominanza sia necessaria: è solo sufficiente (condizione necessaria e sufficiente: $\rho(B)<1$).

Esercizi svolti: Esercizio - Jacobi a mano, dominanza diagonale e stima del numero di iterazioni, Esercizio - Gauss-Seidel a mano e confronto con Jacobi, Esercizio - matrici di iterazione e raggio spettrale di sistemi 2x2 e 3x3, Esercizio - conduzione del calore in una barra con Jacobi e Gauss-Seidel.

Domande d'esame

1. Enunciare e dimostrare il teorema di convergenza dei metodi iterativi $\mathbf x^{(k+1)}=B\mathbf x^{(k)}+\mathbf c$. Traccia. $\mathbf e^{(k)}=B^k\mathbf e^{(0)}$; se $\|B\|<1$ tende a zero; se $\rho(B)\ge1$ un autovettore dà $\lambda^k\mathbf v\not\to0$.
2. Scrivere le matrici di iterazione di Jacobi e Gauss-Seidel e dimostrare che Jacobi converge per dominanza diagonale stretta. Traccia. $B_J=-D^{-1}(L+U)$, $B_{GS}=-(D+L)^{-1}U$; $\|B_J\|_\infty=\max_i\sum_{j\ne i}\frac{|a_{ij}|}{|a_{ii}|}<1$.
3. Stimare quante iterazioni servono per ridurre l'errore di $10^{-6}$ con $\rho(B)=0{,}9$. Traccia. $k\approx\frac6{-\log_{10}0{,}9}=\frac6{0{,}0458}=131$.
4. Perché l'incremento $\|\mathbf x^{(k+1)}-\mathbf x^{(k)}\|$ piccolo non basta come criterio d'arresto? Traccia. $\|\mathbf e^{(k)}\|\le\frac{\|B\|}{1-\|B\|}\|\mathbf x^{(k)}-\mathbf x^{(k-1)}\|$: se $\|B\|\to1$ il fattore esplode.