26 Metodi Iterativi per Sistemi Lineari

Autore/Autrice

Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

I metodi iterativi permettono di risolvere sistemi lineari della forma \bm{A}\bm{x}=\bm{b}, in cui \bm{A} è una matrice e \bm{b} un vettore noto. Sono particolarmente utili nel caso di matrici grandi e sparse, per le quali i metodi diretti, come la fattorizzazione di Gauss, possono richiedere troppo tempo di calcolo e memoria.

Un metodo iterativo aggiorna progressivamente un’approssimazione della soluzione \bm{x} fino a raggiungere una convergenza soddisfacente. Tra i metodi più comuni figurano il metodo di Jacobi, il metodo di Gauss-Seidel e il metodo del gradiente coniugato. Questi procedimenti sfruttano le proprietà della matrice \bm{A} e, in genere, convergono più rapidamente dei metodi diretti quando la matrice è ben condizionata.

Per valutare la qualità di una soluzione approssimata nei metodi iterativi, è fondamentale considerare il residuo, definito come la differenza tra il termine noto e il prodotto tra la matrice e l’approssimazione:

\bm{r} = \bm{b} - \bm{A} \bm{x},

L’obiettivo dei metodi iterativi è ridurre il residuo, avvicinandolo il più possibile a zero. Di solito, la norma del residuo viene utilizzata come criterio di convergenza: il processo iterativo continua finché tale norma non scende sotto una soglia prefissata. Un residuo piccolo indica in genere che l’approssimazione della soluzione è vicina a quella esatta.

Schema iterativo:

Si sceglie un’approssimazione iniziale \bm{x}^0 e si costruisce una successione di approssimazioni \bm{x}^k mediante uno schema iterativo del tipo:

\bm{x}^{k+1} = \bm{G}(\bm{x}^k).

Sotto opportune condizioni, ci aspettiamo che questa successione converga, ovvero:

\bm{x}^* = \lim_{k \to \infty} \bm{x}^k, \qquad \bm{r}^*=\bm{b}-\bm{A}\bm{x}^* = \bm{0}.

In questo contesto, \bm{x}^* rappresenta la soluzione esatta del sistema, e il residuo \bm{r}^* deve annullarsi, indicando che la soluzione trovata soddisfa l’equazione originale.

Questioni da Definire

Quali metodi utilizziamo per costruire gli schemi di iterazione della forma \bm{x}^{k+1} = \bm{G}(\bm{x}^k)?
Quali criteri seguiamo per selezionare la soluzione iniziale \bm{x}^0?
Come determinare l’accuratezza della soluzione approssimata?
Quali strategie adottiamo per arrestare le iterazioni, considerando che non possiamo eseguire un numero infinito di passi?
Come misuriamo l’efficienza di uno schema iterativo?

26.1 Costruzione degli Schemi di Iterazione mediante Splitting

Consideriamo la matrice \bm{A} \in \mathbb{R}^{n \times n} espressa nella forma:

\bm{A} = \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ A_{2,1} & A_{2,2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & A_{n-1,n} \\ A_{n,1} & \cdots & A_{n,n-1} & A_{n,n} \end{pmatrix},

e introduciamo le seguenti matrici:

\begin{aligned} \bm{D} &= \begin{pmatrix} A_{1,1} & 0 & \cdots & 0 \\ 0 & A_{2,2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & A_{n,n} \end{pmatrix}, \\ \bm{L} &= \begin{pmatrix} 0 & 0 & \cdots & 0 \\ A_{2,1} & 0 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ A_{n,1} & \cdots & A_{n,n-1} & 0 \end{pmatrix}, \\ \bm{U} &= \begin{pmatrix} 0 & A_{1,2} & \cdots & A_{1,n} \\ 0 & 0 & \ddots & \vdots \\ \vdots & \ddots & \ddots & A_{n-1,n} \\ 0 & \cdots & 0 & 0 \end{pmatrix}. \end{aligned}

Utilizzando le matrici \bm{D}, \bm{L} e \bm{U}, possiamo decomporre la matrice \bm{A} come segue:

\bm{A} = \bm{D} + \bm{L} + \bm{U}.

Grazie a questa decomposizione o splitting di \bm{A}, il sistema \bm{A}\bm{x} = \bm{b} può essere riscritto in vari modi:

\begin{aligned} (i) &\qquad & \bm{D}\bm{x} &= \bm{b} - (\bm{L} + \bm{U})\bm{x}, \\[1em] (ii) &\qquad & (\bm{D} - \bm{L})\bm{x} &= \bm{b} - \bm{U}\bm{x}, \\[1em] (iii) &\qquad & (\bm{D} - \bm{U})\bm{x} &= \bm{b} - \bm{L}\bm{x}. \end{aligned}

Osserviamo che nella forma (i) la matrice \bm{D} è diagonale, mentre in (ii) la matrice (\bm{D} - \bm{L}) è triangolare inferiore e in (iii) la matrice (\bm{D} - \bm{U}) è triangolare superiore.

26.2 Metodo di Jacobi

La formulazione del sistema \bm{A}\bm{x}=\bm{b} in (i)

\bm{D}\bm{x} = \bm{b} - (\bm{L}+\bm{U})\bm{x}

suggerisce uno schema di iterazione della forma

\boxed{\bm{D} \bm{x}^{k+1} = \bm{b} - (\bm{L}+\bm{U}) \bm{x}^k}

noto come metodo iterativo di Jacobi¹.

26.3 Metodo di Gauss-Seidel

Analogamente la riformulazione in (ii)

(\bm{D}+\bm{L})\bm{x} = \bm{b} - \bm{U}\bm{x}

suggerisce un secondo schema di iterazione della forma

\boxed{(\bm{D}+\bm{L}) \bm{x}^{k+1} = \bm{b} - \bm{U} \bm{x}^k}

noto come metodo iterativo di Gauss²-Seidel³.

Osservazione

A questo punto, con un po’ di creatività, possiamo considerare la riformulazione (iii):

(\bm{D} + \bm{U}) \bm{x} = \bm{b} - \bm{L} \bm{x}

per definire un terzo schema iterativo nel seguente modo:

\boxed{(\bm{D} + \bm{U}) \bm{x}^{k+1} = \bm{b} - \bm{L} \bm{x}^k}.

È importante notare che questo schema può essere interpretato come una variante del metodo di Gauss-Seidel. In effetti, esso deriva dalla forma (ii), ma applicata al sistema lineare con le equazioni e le incognite invertite. Si tratta, quindi, di uno schema di Gauss-Seidel all’indietro.

26.4 Metodo SOR

Separiamo la parte diagonale di \bm{A}, cioè \bm{D}, in due termini pesati con un coefficiente \omega reale e non negativo, da determinare,

\bm{D} = \dfrac{1}{\omega}\bm{D} + \left(1-\dfrac{1}{\omega}\right)\bm{D}.

Introduciamo questi due termini al posto della matrice \bm{D} nella decomposizione utilizzata per il metodo di Gauss-Seidel

\left(\dfrac{1}{\omega}\bm{D}+\bm{L}\right)\bm{x} = \bm{b} - \left[\left(1-\dfrac{1}{\omega}\right)\bm{D}+\bm{U}\right]\bm{x}.

Si ottiene quindi un nuovo schema iterativo

\boxed{ \left(\dfrac{1}{\omega}\bm{D}+\bm{L}\right)\bm{x}^{k+1} = \bm{b} - \left(\dfrac{\omega-1}{\omega}\bm{D}+\bm{U}\right)\bm{x}^k }

noto come Successive Over Relaxation method, da cui l’acronimo SOR. Si noti che il metodo SOR dipende da un parametro \omega che deve essere individuato correttamente per garantire convergenza ed efficienza.

26.5 Derivazione per componenti dei metodi di Jacobi, Gauss–Seidel e SOR

Partiamo dallo splitting \bm{A}=\bm{D}+\bm{L}+\bm{U} e dal sistema componente per componente

\sum_{j=1}^n A_{ij} x_j = b_i,\qquad i=1,\dots,n.

Separiamo il termine diagonale A_{ii}x_i dalla somma:

A_{ii} x_i + \sum_{j\ne i} A_{ij} x_j = b_i.

Risolviamo per x_i ottenendo la forma fondamentale che porta ai metodi iterativi.

Metodo di Jacobi

Nel metodo di Jacobi tutte le componenti del nuovo vettore \bm{x}^{k+1} sono calcolate usando solo i valori vecchi \bm{x}^k. Da

A_{ii} x_i^{k+1} + \sum_{j\ne i} A_{ij} x_j^{k} = b_i,

otteniamo per ogni i

x_i^{k+1} = \dfrac{1}{A_{ii}}\Bigl(b_i - \sum_{j\ne i} A_{ij} x_j^{k}\Bigr) = \dfrac{1}{A_{ii}}\Bigl(b_i - \sum_{j=1}^{i-1} A_{ij} x_j^{k} - \sum_{j=i+1}^{n} A_{ij} x_j^{k}\Bigr).

Tutti gli x_i^{k+1} possono essere calcolati in parallelo.

Jacobi

Input: A[n,n], b[n], x^(0), tol, max_iter
Output: x ≈ soluzione di A x = b

k = 0
x_new = x^(0)

while k < max_iter:
    for i = 1 to n:
        sum = 0
        for j = 1 to n:
            if j != i:
                sum += A[i,j] * x[j]
        end
        x_new[i] = (b[i] - sum) / A[i,i]
    end
    if ||x_new - x|| < tol:
        break
    x = x_new
    k = k + 1
end

Metodo di Gauss–Seidel

Nel metodo di Gauss–Seidel, quando si calcola x_i^{k+1} si usano le componenti già aggiornate per gli indici j<i e le componenti vecchie per j>i. Partendo da

A_{ii} x_i^{k+1} + \sum_{j<i} A_{ij} x_j^{k+1} + \sum_{j>i} A_{ij} x_j^{k} = b_i,

si ottiene

x_i^{k+1} = \dfrac{1}{A_{ii}}\Bigl(b_i - \sum_{j=1}^{i-1} A_{ij} x_j^{k+1} - \sum_{j=i+1}^{n} A_{ij} x_j^{k}\Bigr).

Per i=1 la somma \sum_{j<1} è nulla, e per i=n la somma \sum_{j>n} è nulla.

Gauss–Seidel

Input: A[n,n], b[n], x^(0), tol, max_iter
Output: x ≈ soluzione di A x = b

k = 0
x = x^(0)

while k < max_iter:
    for i = 1 to n:
        sum1 = sum_{j=1}^{i-1} A[i,j] * x[j]   # valori aggiornati
        sum2 = sum_{j=i+1}^{n} A[i,j] * x[j]   # valori vecchi
        x[i] = (b[i] - sum1 - sum2) / A[i,i]
    end
    if ||A*x - b|| < tol:
        break
    k = k + 1
end

Metodo SOR (Successive Over-Relaxation)

Il metodo SOR introduce un parametro di rilassamento \omega. La derivazione componente-wise combina il vecchio valore x_i^k con l’aggiornamento di Gauss–Seidel \tilde x_i^{k+1}:

\tilde x_i^{k+1} = \dfrac{1}{A_{ii}}\Bigl(b_i - \sum_{j=1}^{i-1} A_{ij} x_j^{k+1} - \sum_{j=i+1}^{n} A_{ij} x_j^{k}\Bigr),

poi applichiamo il rilassamento:

x_i^{k+1} = (1-\omega) x_i^{k} + \omega \tilde x_i^{k+1},

ossia, sostituendo \tilde x_i^{k+1}:

x_i^{k+1} = (1-\omega)x_i^{k} + \dfrac{\omega}{A_{ii}}\Bigl(b_i - \sum_{j=1}^{i-1} A_{ij} x_j^{k+1} - \sum_{j=i+1}^{n} A_{ij} x_j^{k}\Bigr).

\omega=1 riduce al metodo sdi Gauss–Seidel.
0<\omega<1 → under-relaxation.
1<\omega<2 → over-relaxation (talvolta accelera la convergenza).

SOR (Successive Over-Relaxation)

Input: A[n,n], b[n], x^(0), omega, tol, max_iter
Output: x ≈ soluzione di A x = b

k = 0
x = x^(0)

while k < max_iter:
    for i = 1 to n:
        sum1 = sum_{j=1}^{i-1} A[i,j] * x[j]   # valori aggiornati
        sum2 = sum_{j=i+1}^{n} A[i,j] * x[j]   # valori vecchi
        x_i_GS = (b[i] - sum1 - sum2) / A[i,i]
        x[i] = (1 - omega) * x[i] + omega * x_i_GS
    end
    if ||A*x - b|| < tol:
        break
    k = k + 1
end

26.6 Generalizzazione

Scriviamo \bm{A}, matrice quadrata non singolare di ordine n, come differenza di due matrici \bm{A}=\bm{P}-\bm{Q} in cui \bm{P} è non singolare e facile da invertire.

(\bm{P}-\bm{Q})\bm{x} = \bm{b}, \quad\Rightarrow\quad \bm{P}\bm{x} = \bm{Q}\bm{x} + \bm{b},

Da questa decomposizione si ricava in maniera naturale lo schema di iterazione

\boxed{ \bm{P}\bm{x}^{k+1} = \bm{Q}\bm{x}^k + \bm{b}, \quad\Rightarrow\quad \bm{x}^{k+1} = \bm{P}^{-1}\bm{Q}\bm{x}^k + \bm{P}^{-1}\bm{b}, }

Definizione 26.1 (Matrice di iterazione) La matrice \bm{P}^{-1}\bm{Q}, che possiamo anche scrivere come segue

\bm{P}^{-1}\bm{Q} = \bm{P}^{-1}\Big(\bm{P}-\bm{A}\Big) = \bm{I}-\bm{P}^{-1}\bm{A}

prende il nome di matrice di iterazione.

Ovviamente, è indifferente scrivere lo schema di iterazione come segue

\bm{x}^{k+1} = (\bm{I}-\bm{P}^{-1}\bm{A}) \bm{x}^k + \bm{P}^{-1} \bm{b},

Osservazione

Tutti i metodi precedenti introdotti finora, Jacobi, Gauss-Seidel, SOR, sono casi particolari del metodo generale.

Tabella 26.1: Quadro riassuntivo

	Jacobi	Gauss-Seidel	SOR
\bm{P}	\bm{D}	\bm{D}+\bm{L}	\dfrac{1}{\omega}\bm{D}+\bm{L}
\bm{Q}	-(\bm{L}+\bm{U})	-\bm{U}	\dfrac{1-\omega}{\omega}\bm{D}-\bm{U}
\bm{I}-\bm{P}^{-1}\bm{A}	-\bm{D}^{-1}(\bm{L}+\bm{U})	-(\bm{D}+\bm{L})^{-1}\bm{U}	\left(\dfrac{\bm{D}}{\omega}+\bm{L}\right)^{-1}\left(\dfrac{1-\omega}{\omega}\bm{D}-\bm{U} \right)

26.7 Convergenza degli schemi iterativi

Definizione 26.2 (Convergenza) Uno schema iterativo della forma

\bm{x}^{k+1} = \bm{G}( \bm{x}^k )

si dice convergente se esiste un vettore \bm{x}^* ed un norma vettoriale \left|\left|\cdot\right|\right| tali che la successione delle iterate costruita a partire da un qualsiasi vettore iniziale \bm{x}^0 converge a \bm{x}^* nella norma considerata

\lim_{k\to\infty} \left|\left|\bm{x}^k-\bm{x}^*\right|\right| = \bm{0}.

Osservazione

Si noti che \bm{x}^* deve soddisfare la condizione

\bm{x}^* = \bm{G}( \bm{x}^*)

Uno schema convergente nel senso della definizione precedente ha convergenza globale. Per problemi non-lineari uno schema iterativo può essere convergente solo per un sottoinsieme di vettori di \mathbb{R}^n. In tal caso si parla di convergenza locale. Per problemi lineari, si assume sempre che la convergenza sia globale.

Questione:

sotto quali condizioni \bm{x}^k converge alla soluzione \bm{x}^*?

Definizione 26.3 (Errore) Definiamo errore \bm{\varepsilon}^k la differenza tra la soluzione esatta \bm{x}^{*} e la soluzione approssimata \bm{x}^k; cioè

\bm{\varepsilon}^k = \bm{x}^{*} - \bm{x}^k,

Quindi la condizione che lo schema iterativo sia convergente si può esprimere anche chiedendo che l’errore sia infinitesimo in norma

\lim_{k\mapsto\infty}\left|\left|\bm{\varepsilon}^k\right|\right| = \bm{0}.

26.8 Come si trasforma l’errore secondo lo schema iterativo generale?

Osserviamo che dalle due espressioni che definiscono \bm{x}^* e \bm{x}^{k+1}

\begin{aligned} \bm{P}\bm{x}^* &= \bm{Q}\bm{x}^* + \bm{b}, \\ \bm{P}\bm{x}^{k+1} &= \bm{Q}\bm{x}^k + \bm{b}, \end{aligned}

si ottiene per differenza

\bm{P}\Big(\bm{x}^*-\bm{x}^{k+1}\Big) = \bm{Q}\Big(\bm{x}^*-\bm{x}^k\Big),

da cui segue iterando

\begin{aligned} \bm{\varepsilon}^{k+1} &= \Big(\bm{P}^{-1}\bm{Q}\Big)\bm{\varepsilon}^k \\ &= \Big(\bm{P}^{-1}\bm{Q}\Big)^2\bm{\varepsilon}^{k-1}\\ &= \Big(\bm{P}^{-1}\bm{Q}\Big)^3\bm{\varepsilon}^{k-2}\\ &= \ldots\\ &= \Big(\bm{P}^{-1}\bm{Q}\Big)^{k+1}\bm{\varepsilon}^0. \end{aligned}

Osservazione

La convergenza di un metodo iterativo può essere definita richiedendo che l’errore sia infinitesimale in una certa norma vettoriale.

Tuttavia, tale definizione dell’errore presuppone la conoscenza della soluzione esatta \bm{x}^*. Si osserva che:

\bm{A}\bm{\varepsilon}^k = \bm{A}\bm{x}^* - \bm{A}\bm{x}^k = \bm{b} - \bm{A}\bm{x}^k = \bm{r}^k,

da cui si deduce immediatamente che:

\bm{x}^k \mapsto \bm{x}^* \qquad \Leftrightarrow \qquad \bm{r}^k \mapsto 0.

Pertanto, studiare la convergenza di un metodo iterativo è equivalente a analizzare quando il residuo tende a zero.

Trasformazione del residuo secondo lo schema iterativo generale

Possiamo scrivere la matrice \bm{A} come:

\bm{A} = \bm{P} - \bm{Q} = \left(\bm{I} - \bm{Q}\bm{P}^{-1}\right)\bm{P}.

Da questa espressione, con alcuni passaggi ovvi, otteniamo:

\begin{aligned} \bm{A}\bm{P}^{-1}\bm{Q} &= \left(\bm{I} - \bm{Q}\bm{P}^{-1}\right)\bm{P} \bm{P}^{-1}\bm{Q}, \\ &= \bm{Q} - \bm{Q}\bm{P}^{-1}\bm{Q}, \\ &= \bm{Q}\bm{P}^{-1}\bm{P} - \bm{Q}\bm{P}^{-1}\bm{Q}, \\ &= \bm{Q}\bm{P}^{-1}\left(\bm{P} - \bm{Q}\right), \\ &= \bm{Q}\bm{P}^{-1}\bm{A}. \end{aligned}

Infine, dalla relazione \bm{A}\bm{P}^{-1}\bm{Q} = \bm{Q}\bm{P}^{-1}\bm{A}, possiamo dedurre:

\begin{aligned} \bm{r}^{k+1} &= \bm{b} - \bm{A}\bm{x}^{k+1}, \\ &= \bm{b} - \bm{A}\bm{P}^{-1}\bm{Q}\bm{x}^k - \bm{A}\bm{P}^{-1}\bm{b}, \\ &= \bm{b} - \bm{Q}\bm{P}^{-1}\bm{A}\bm{x}^k - \left(\bm{P} - \bm{Q}\right)\bm{P}^{-1}\bm{b}, \\ &= \bm{b} - \bm{Q}\bm{P}^{-1}\bm{A}\bm{x}^k - \bm{b} + \bm{Q}\bm{P}^{-1}\bm{b}, \\ &= \bm{Q}\bm{P}^{-1}(\bm{b} - \bm{A}\bm{x}^k), \\ &= \bm{Q}\bm{P}^{-1}\bm{r}^k. \end{aligned}

Teorema 26.1 (Convergenza) Sia \bm{A} una qualsiasi matrice quadrata. Allora le seguenti condizioni sono equivalenti:

\lim_{k\mapsto\infty}\bm{A}^k = \bm{0}.
\lim_{k\mapsto\infty}\bm{A}^k\bm{v} = \bm{0} per ogni vettore \bm{v}.
\varrho(\bm{A}) < 1.
\left|\left|\bm{A}\right|\right|<1 per almeno una norma matriciale.

Teorema 26.2 (Convergenza) Dato un sistema lineare \bm{A}\bm{x}^*=\bm{b} ed una decomposizione \bm{A}=\bm{P}-\bm{Q} con \bm{P} non singolare, allora il metodo iterativo

\bm{P}\bm{x}^{k+1} = \bm{Q}\bm{x}^k + \bm{b},

converge alla soluzione \bm{x}^* se e solo se vale una delle due:

\begin{aligned} (i)&\qquad \boxed{\varrho(\bm{Q}\bm{P}^{-1}) < 1} \\ (ii)&\qquad \boxed{\varrho(\bm{P}^{-1}\bm{Q}) < 1} \end{aligned}

Dimostrazione. Dimostriamo separatamente le due tesi del teorema.

(i) Il residuo k-esimo può essere scritto come \bm{r}^k = (\bm{Q}\bm{P}^{-1})^k \bm{r}^0, per cui dal teorema~1 segue che \varrho(\bm{Q}\bm{P}^{-1}) < 1 \quad\implies\quad \lim_{k\mapsto\infty}\bm{r}^k=\bm{0}.

Viceversa, se il metodo converge, si ha sempre dal teorema Teorema 26.1 che \begin{aligned} \text{ per ogni }\bm{r}^0, &\quad \lim_{k\mapsto\infty}\bm{r}^k = \lim_{k\mapsto\infty}(\bm{Q}\bm{P}^{-1})^k\bm{r}^0=\bm{0} \;\implies\\ &\implies\quad\varrho(\bm{Q}\bm{P}^{-1})<1. \end{aligned}
(ii)
Si dimostra analogamente sfruttando il modo in cui si trasforma l’errore durante il procedimento iterativo. Poiché \bm{\varepsilon}^k = (\bm{P}^{-1}\bm{Q})^k \bm{\varepsilon}^0, dal teorema Teorema 26.1 segue che \varrho(\bm{P}^{-1}\bm{Q}) < 1 \quad\implies\quad \lim_{k\mapsto\infty}\bm{\varepsilon}^k=\bm{0}. Viceversa, se il metodo converge, si ha sempre dal teorema~1 che \begin{aligned} \text{ per ogni }\bm{\varepsilon}^0,&\quad \lim_{k\mapsto\infty}\bm{\varepsilon}^k= \lim_{k\mapsto\infty}(\bm{P}^{-1}\bm{Q})^k\bm{\varepsilon}^0=\bm{0} \;\implies\\ &\implies\quad\varrho(\bm{P}^{-1}\bm{Q})<1. \end{aligned}

Corollario 26.1 (Corollario) Il metodo di Jacobi converge alla soluzione se e solo se

\varrho\left(\bm{D}^{-1}(\bm{L}+\bm{U})\right) < 1.

Corollario 26.2 (Corollario) Il metodo di Gauss-Seidel converge se e solo se

\varrho\left((\bm{D}+\bm{L})^{-1}\bm{U}\right) < 1.

Corollario 26.3 (Corollario) Il metodo SOR converge se e solo se

\varrho\left( \left(\dfrac{\bm{D}}{\omega}+\bm{L}\right)^{-1} \left(\dfrac{1-\omega}{\omega}\bm{D}-\bm{U} \right)\right) < 1.

Definizione 26.4 (Raggio spettrale) Il raggio spettrale della matrice di iterazione del metodo SOR soddisfa

\varrho\left( \left( \dfrac{1-\omega}{\omega}\bm{D} - \bm{U} \right) \left( \dfrac{\bm{D}}{\omega}+\bm{L} \right)^{-1}\right)\geq\left\lvert 1-\omega\right\rvert.

Teorema 26.3 (SOR convergence) Scegliere \omega\in(0,2) è condizione necessaria per la convergenza del metodo SOR.

Dimostrazione. Calcoliamo il determinante della matrice di iterazione del metodo SOR.

\begin{aligned} \left|\begin{matrix} \left[\dfrac{\bm{D}}{\omega}+\bm{L}\right]^{-1} \left[\dfrac{1-\omega}{\omega}\bm{D} - \bm{U}\right] \end{matrix}\right| & = \left|\begin{matrix}\dfrac{\bm{D}}{\omega}+\bm{L}\end{matrix}\right|^{-1}\,\, \left|\begin{matrix}\dfrac{1-\omega}{\omega}\bm{D}-\bm{U}\end{matrix}\right| \\ & = \dfrac { \prod_{i=1}^{n}\dfrac{1-\omega}{\omega} \bm{D}_{i,i} } { \prod_{i=1}^{n}\dfrac{1}{\omega}\bm{D}_{i,i} } \\ & = (1-\omega)^n. \end{aligned}

Osserviamo ora che se \bm{M} è una qualsiasi matrice con autovalori \lambda_i e raggio spettrale \varrho(\bm{M}), allora

\varrho(\bm{M})^n \geq \prod_{i=1}^{n} \left\lvert\lambda_i\right\rvert \geq \left\lvert\det{\bm{M}}\right\rvert.

Prendendo \bm{M} uguale alla matrice di iterazione del metodo SOR si ha

\begin{aligned} \varrho\left( \left( \dfrac{\bm{D}}{\omega}+\bm{L} \right)^{-1} \left( \dfrac{1-\omega}{\omega}\bm{D} - \bm{U} \right) \right)& \geq \\ \left\lvert\det{ \left( \dfrac{\bm{D}}{\omega}+\bm{L} \right)^{-1} \left( \dfrac{1-\omega}{\omega}\bm{D} - \bm{U} \right) } \right\rvert^{\dfrac{1}{n}}&=\left\lvert 1-\omega\right\rvert. \end{aligned}

Indipendentemente dalla matrice \bm{A}=\bm{D}-\bm{L}-\bm{U}, dato che \omega\not\in(0,2) \quad\implies\quad \left\lvert 1-\omega\right\rvert\geq 1.

abbiamo che il raggio spettrale è >1 se \omega\not\in(0,2).

26.9 Controllo della Convergenza

Se \bm{\varepsilon}^{k-1}\neq\bm{0}, la quantità \left|\left|\bm{\varepsilon}^k\right|\right|/\left|\left|\bm{\varepsilon}^{k-1}\right|\right| esprime la riduzione dell’errore al k-esimo passo. La media geometrica delle riduzioni dell’errore sui primi k passi

\sigma_k = \left( \dfrac{\left|\left|\bm{\varepsilon}^1\right|\right|}{\left|\left|\bm{\varepsilon}^0\right|\right|}\, \dfrac{\left|\left|\bm{\varepsilon}^2\right|\right|}{\left|\left|\bm{\varepsilon}^1\right|\right|}\, \ldots \dfrac{\left|\left|\bm{\varepsilon}^k\right|\right|}{\left|\left|\bm{\varepsilon}^{k-1}\right|\right|} \right)^{1/k} = \left( \dfrac{\left|\left|\bm{\varepsilon}^k\right|\right|}{\left|\left|\bm{\varepsilon}^0\right|\right|} \right)^{1/k}

esprime la riduzione media per passo dell’errore.

Dalle proprietà delle norme naturali segue che

\bm{\varepsilon}^k = (\bm{P}^{-1}\bm{Q})^k \bm{\varepsilon}^0 \quad\implies\quad \left|\left|\bm{\varepsilon}^k\right|\right|\leq\left|\left|(\bm{P}^{-1}\bm{Q})^k\right|\right|\left|\left|\bm{\varepsilon}^0\right|\right|

e quindi si ottiene la relazione

\sigma_k\leq \Big(\left|\left|(\bm{P}^{-1}\bm{Q})^k\right|\right|\Big)^{\dfrac{1}{k}}.

Quindi \sigma_k stima la velocità di convergenza sulle prime k iterazioni ed è controllato da una quantità che dipende dalla matrice di iterazione dello schema. Come possiamo stimare questa quantità?

Stima della velocità di convergenza

Teorema 26.4 (Limite potenza e norma) Sia \bm{M}\in\mathbb{K}^{m\times n} e \left|\left|\cdot\right|\right| una qualunque norma indotta. Allora

\lim_{k\mapsto\infty}\left|\left|(\bm{M})^k\right|\right|^{\dfrac{1}{k}}= \varrho(\bm{M}).

Osservazione

La quantità \varrho(\bm{M}) non dipende dalla norma matriciale utilizzata e nemmeno dall’indice di iterazione k. Possiamo assumerla come una stima della velocità di convergenza sui nostri metodi iterativi ponendo \bm{M}=\bm{P}^{-1}\bm{Q}, cioè uguale alla matrice di iterazione.

Definizione 26.5 (Tasso asintotico di convergenza) Si definisce tasso asintotico di convergenza di uno schema iterativo la costante

\mathcal{C} = -\log_{10}{\varrho(\bm{P}^{-1}\bm{Q})}.

La costante \mathcal{C} è all’incirca il numero di iterazioni richieste per ridurre l’errore di un fattore \dfrac{1}{10}, cioè grossomodo per avere ottenere una cifra decimale in più. Infatti

\Big(\varrho(\bm{P}^{-1}\bm{Q})\Big)^k\approx\dfrac{1}{10} \quad\implies\quad k\approx -\dfrac{1}{\log_{10}{\varrho(\bm{P}^{-1}\bm{Q})}}.

Criteri di Arresto

Sia fissata una tolleranza \eta sull’errore. I criteri più comunemente usati sfruttano un controllo sulla variazione assoluta e relativa delle iterate al generico passo k

\begin{aligned} \text{\textbf{IF}} \quad\left|\left|\bm{x}^k-\bm{x}^{k-1}\right|\right| &\leq\eta\quad & \text{\textbf{THEN STOP}} \\ \text{\textbf{IF}} \quad\left|\left|\bm{x}^k-\bm{x}^{k-1}\right|\right| &\leq \eta\left|\left|\bm{x}^k\right|\right|\quad & \text{\textbf{THEN STOP}} \end{aligned}

Tuttavia è importante rendersi conto che queste condizioni non garantiscono che la soluzione sia stata approssimata con accuratezza \eta

Infatti, si ha che

\begin{aligned} \bm{x}^k-\bm{x}^{k-1} &= \bm{x}^k - \bm{x}^* + \bm{x}^* -\bm{x}^{k-1} \\ &= -\bm{\varepsilon}^k + \bm{\varepsilon}^{k-1} \\ &= -\bm{P}^{-1}\bm{Q}\bm{\varepsilon}^{k-1} + \bm{\varepsilon}^{k-1} \\ &= \left(\bm{I}-\bm{P}^{-1}\bm{Q}\right)\bm{\varepsilon}^{k-1} \end{aligned}

per cui possiamo esprimere l’errore al passo k-1 in funzione della variazione assoluta delle iterate al passo k

\bm{\varepsilon}^{k-1} = \left(\bm{I}-\bm{P}^{-1}\bm{Q}\right)^{-1}\left( \bm{x}^k-\bm{x}^{k-1}\right).

Dalla relazione

\bm{\varepsilon}^{k-1} = \left(\bm{I}-\bm{P}^{-1}\bm{Q}\right)^{-1}\left( \bm{x}^k-\bm{x}^{k-1}\right), \qquad\text{(ripetiamo)}

passando alle norme se \left|\left|\bm{P}^{-1}\bm{Q}\right|\right|\leq 1 si ha

\begin{aligned} \left|\left|\bm{\varepsilon}^{k-1}\right|\right| &\leq\left|\left|(\bm{I}-\bm{P}^{-1}\bm{Q})^{-1}\right|\right|\,\left|\left|\bm{x}^k-\bm{x}^{k-1}\right|\right| \\ &\leq\dfrac{\left|\left|\bm{x}^k-\bm{x}^{k-1}\right|\right|}{1-\left|\left|\bm{P}^{-1}\bm{Q}\right|\right|} \end{aligned}

per cui se \left|\left|\bm{P}^{-1}\bm{Q}\right|\right|\approx 1 l’errore in norma \left|\left|\bm{\varepsilon}^{k-1}\right|\right| può essere ancora grande.

26.10 Quale metodo converge più in fretta

I criteri visti finora dicono se un metodo converge, non quale dei due converge più in fretta. Il confronto si fa sul raggio spettrale della matrice di iterazione, e il risultato chiave è il seguente: fra due splitting della stessa matrice, converge più rapidamente quello che lascia meno matrice in \bm{Q}.

Confronto fra due splitting

Siano \bm{A}=\bm{P}_1-\bm{Q}_1=\bm{P}_2-\bm{Q}_2 due splitting debolmente non negativi della stessa matrice. Se \bm{A}^{-1}\geq\bm{0} e \bm{Q}_2\geq\bm{Q}_1 allora

\varrho(\bm{P}_1^{-1}\bm{Q}_1)\leq\varrho(\bm{P}_2^{-1}\bm{Q}_2)<1.

La dimostrazione è in appendice, Teorema C.4.

Applichiamolo a Jacobi e Gauss-Seidel. Con \bm{A}=\bm{D}-\bm{L}-\bm{U} si ha

\bm{Q}_{\text{J}}=\bm{L}+\bm{U}, \qquad \bm{Q}_{\text{GS}}=\bm{U}.

Se \bm{A} è una M-matrice (Definizione 25.1) i suoi elementi fuori diagonale sono non positivi, quindi \bm{L}\geq\bm{0} e \bm{U}\geq\bm{0}, e inoltre \bm{A}^{-1}\geq\bm{0} per definizione. Allora

\bm{Q}_{\text{GS}}=\bm{U}\;\leq\;\bm{L}+\bm{U}=\bm{Q}_{\text{J}} \qquad\Longrightarrow\qquad \varrho\big(\bm{P}_{\text{GS}}^{-1}\bm{Q}_{\text{GS}}\big) \leq \varrho\big(\bm{P}_{\text{J}}^{-1}\bm{Q}_{\text{J}}\big)<1,

cioè su una M-matrice Gauss-Seidel non è mai più lento di Jacobi, e i due convergono o divergono insieme. Si noti che l’ipotesi è necessaria: senza di essa esistono matrici per cui vale il contrario.

Lo stesso criterio ordina famiglie più ampie di metodi. L’appendice Appendice C inquadra Jacobi, Gauss-Seidel e EWA come casi particolari di un’unica costruzione per splitting, riassunta in Tabella C.1, e ne confronta le velocità di convergenza.

In sintesi

Ogni metodo iterativo nasce da uno splitting \bm{A}=\bm{M}-\bm{N} con \bm{M} facile da invertire; la matrice di iterazione è \bm{B}=\bm{M}^{-1}\bm{N}.
Jacobi, Gauss-Seidel e SOR differiscono solo per la scelta di \bm{M}: diagonale, triangolare inferiore, combinazione pesata da \omega.
La condizione necessaria e sufficiente di convergenza è \rho(\bm{B})<1, per ogni dato iniziale; \|\bm{B}\|<1 in una qualsiasi norma è solo sufficiente.
L’errore e il residuo si trasformano con la stessa matrice \bm{B}: il residuo è osservabile, l’errore no, ed è questo che si usa come criterio d’arresto.
Il parametro \omega di SOR va scelto, non indovinato: per matrici consistentemente ordinate esiste un \omega_{\mathrm{opt}} esplicito in funzione di \rho(\bm{B}_{\mathrm{Jacobi}}).