Calcolo numerico/Principali metodi iterativi

Metodo di Jacobi

Dividiamo $A$ in tre matrici $D+L+U$, dove $D$ è la matrice diagonale, $L$ la triangolare inferiore stretta (con 0 sulla diagonale principale), $U$ triangolare superiore stretta.

Nel caso del metodo di Jacobi $A=M-N$ con $M=D$ e $N=-(L+U)$ ($N$ è uguale ad $A$ con zeri sulla diagonale principale).

$A{\boldsymbol {x}}=(M-N){\boldsymbol {x}}={\boldsymbol {b}}$

$\longrightarrow \,{\boldsymbol {x}}=M^{-1}N{\boldsymbol {x}}+M^{-1}{\boldsymbol {b}}$

$A=M^{-1}N,\quad {\boldsymbol {c}}=M^{-1}{\boldsymbol {b}}$

e tenendo conto delle espressioni di $M$ e $N$:

$P=-D^{-1}(L+U),\,{\boldsymbol {c}}=D^{-1}{\boldsymbol {b}}$

Indicando con $P_{j}$ la matrice d'iterazione del metodo di Jacobi si ha:

${\boldsymbol {x}}_{k}=P_{j}{\boldsymbol {x}}_{k-1}+{\boldsymbol {c}}_{j}$

e moltiplicando per $D$:

$D{\boldsymbol {x}}_{k}=DP_{j}{\boldsymbol {x}}_{k-1}+D{\boldsymbol {c}}_{j}$

e sostituendo le espressioni ricavate prima:

$D{\boldsymbol {x}}_{k}=-(L+U){\boldsymbol {x}}_{k-1}+{\boldsymbol {b}}$

Considerando la $i$-esima equazione si ha:

$x_{i}={\frac {-\sum _{j\neq i}A_{ij}x_{j}^{(k-1)}+b_{i}}{A_{ii}}}$

Metodo di Gauss-Seidel

In questo caso si pone $M=D+L$, invece $N=-U$, e si ricava:

$P=-(D+L)^{-1}U,\quad {\boldsymbol {c}}=(D+L)^{-1}{\boldsymbol {b}}$

Come prima, moltiplichiamo per $M$ la relazione

${\boldsymbol {x}}_{k}=P{\boldsymbol {x}}_{k-1}+{\boldsymbol {b}}$

e otteniamo

$(D+L){\boldsymbol {x}}_{k}=-U{\boldsymbol {x}}_{k-1}+{\boldsymbol {b}}$

cioè

$D{\boldsymbol {x}}_{k}=-L{\boldsymbol {x}}_{k}-U{\boldsymbol {x}}_{k-1}+{\boldsymbol {b}}$

Considerando l'equazione $i$-esima,

$x_{i}^{(k)}={\frac {-\sum _{j<i}A_{ij}x_{j}^{(k)}-\sum _{j>i}A_{ij}x_{j}^{(k-1)}}{A_{ii}}}$

Teniamo conto che nel prodotto con la triangolare inferiore si utilizzano le informazioni già aggiornate al passo $k$-esimo, e sopravvivono solo gli elementi con colonne di indice minori di $i$. Nell'altro prodotto invece si utilizzano le informazioni ricavate al passo $k-1$.

Matrici diagonalmente dominanti e irriducibili

Definizione

Data una matrice $A\in \mathbb {C} _{n,n}$ diciamo che:

1. $A$ è diagonalmente dominante in senso stretto per righe se vale

$|A_{ii}|>\sum _{j\neq i}|A_{ij}|\quad i=1,\dots ,n$

2. $A$ è diagonalmente dominante in senso stretto per colonne se vale

$|A_{ii}|>\sum _{i\neq j}|A_{ji}|$

3. $A$ è diagonalmente dominante in senso debole per righe se

$|A_{ii}|\geq \sum _{j\neq i}|A_{ij}|$

e se esiste almeno un indice ${\bar {i}}$ tale che

$|A_{{\bar {i}},{\bar {i}}}|>\sum _{j=1}^{n}|A_{{\bar {i}},j}|$

Definizione

$A\in \mathbb {C} _{n,n}$ si dice riducibile se esiste $\pi$ matrice di permutazione tale che

$B=\pi A\pi ^{t}$

può essere descritta a blocchi come:

${\begin{array}{cc}A_{11}&A_{12}\\0&A_{22}\end{array}}$

con il primo blocco fatto da $k$ righe e l'ultimo da $n-k$ righe.

Se $A,B$ sono come sopra, dato il sistema lineare

$B{\boldsymbol {x}}={\boldsymbol {b}}$

esso si può riscrivere come:

$\left({\begin{array}{cc}A_{11}&A_{12}\\0&A_{22}\end{array}}\right)\times \left({\begin{array}{c}x_{1}\\x_{2}\end{array}}\right)=\left({\begin{array}{c}b_{1}\\b_{2}\end{array}}\right)$

allora si ottengono le equazioni

${\begin{cases}A_{11}x_{1}+A_{12}x_{2}=b_{1}\\A_{22}x_{2}=b_{2}\end{cases}}$

e risolvendo prima la seconda e sostituendo poi $x_{2}$ nella prima otteniamo un sistema disaccoppiato.

Esempio 2.11

Consideriamo le matrici

${\begin{array}{ccc}1&3&0\\0&2&-1\\0&0&2\end{array}}$

Associamo questa matrice di ordine $n=3$ a un grafo con tre nodi $1,2,3$. Per ogni entrata della matrice, quando $A_{ij}\neq 0$ collego il nodo $i$ con il nodo $j$ con un arco orientato da $i$ a $j$. Ad esempio, $A_{ii}=1$ quindi traccio un arco che entra ed esce nel nodo 1. Nel caso di questa matrice, possiamo andare dal nodo 1 al 2 seguendo l'arco orientato tracciato in corrispondenza di $A_{12}=3$, possiamo andare dal nodo 2 al 3 lungo l'arco tracciato in corrispondenza di $A_{23}=-1$, ma non possiamo andare in nessun modo dal nodo 3 all'1 perché $A_{31}=0$.

Consideriamo invece la matrice

${\begin{array}{ccc}1&3&0\\0&2&-1\\-1&0&2\end{array}}$

Nel caso della seconda matrice, c'è un cammino chiuso che connette tutti i nodi.

Teorema 2.14

Una matrice $A$ è irriducibile se e solo se il grafo associato è fortemente connesso ovvero se e solo se esiste un cammino dal nodo $i$ al nodo $j$ per ogni coppia di nodi $i,j$.

Quindi la prima matrice è riducibile e la seconda non lo è.

Condizioni sufficienti di convergenza

Le condizioni sufficienti per la convergenza del metodo di Jacobi sono le seguenti:

1. $A$ è diagonalmente dominante in senso stretto (per righe o colonne)

Dim. Sappiamo che

$\vert A\vert ^{\infty }=\max _{i=1}^{n}[\sum _{j=1}^{n}|A_{ij}|]$

Nel metodo di Jacobi la matrice di permutazione è

$P=M^{-1}N=-D^{-1}(L+U)$

Se calcolo la norma infinito della matrice di Jacobi $P_{j}$, si ha:

$\vert P\vert \leq \vert D^{-1}\vert *\vert L+U\vert$

Sfruttando il fatto che $A$ è diagonalmente dominante in senso stretto sulle righe, si ha:

$\vert L+U\vert \leq \max _{i}|A_{ii}|$

(infatti la somma di ogni riga privata dell'elemento diagonale è minore dell'elemento diagonale)

$\vert D^{-1}\vert ={\frac {1}{\max _{i}|A_{ii}|}}$

quindi

$\vert P\vert \leq {\frac {1}{\max _{i}|A_{ii}|}}*\max _{i}|A_{ii}|\leq 1$

$\rho (P)<1$ e il metodo converge.

2. $A$ è diagonalmente dominante in senso debole (per righe o per colonne)

3. $A$ è irriducibile.

Nelle stesse ipotesi converge anche il metodo di Gauss-Seidel, però per questo metodo vale anche la seguente proposizione:

Proposizione 2.5

Data una matrice $A$ $m\times n$ hermitiana, non singolare il metodo di Gauss-Seidel converge se e solo se $A$ è definita positiva (la fattorizzazione di Cholesky, che determina se una matrice è definita positiva oppure no, è un buon modo per determinare se il metodo di Gauss-Seidel converge o no).

Teorema 2.15

Sia $A\in \mathbb {C} _{n,n}$ tridiagonale, siano

${\boldsymbol {a}}=(a_{1},\dots ,a_{n}){\hbox{ diagonale principale}}$

${\boldsymbol {c}}=(c_{1},\dots ,c_{n-1}){\hbox{ sopradiagonale}}$

${\boldsymbol {b}}=(b_{2},\dots ,b_{n}){\hbox{ sottodiagonale}}$

Se $\lambda$ è un autovalore della matrice di iterazione di Jacobi, allora $\lambda ^{2}$ è un autovalore per $P_{g}$ (matrice di iterazione di Gauss-Seidel). Inoltre, se $\lambda \neq 0$ è un autovalore di $P_{g}$, allora $\pm {\sqrt {\lambda }}$ è un autovalore di $P_{j}$.

Allora per le matrici tridiagonali, $\rho (P_{g})=\rho (P_{j})^{2}$, allora Gauss-Seidel è convergente se e solo se Jacobi è convergente. Inoltre se Gauss-Seidel converge, ha un tasso di convergenza del doppio rispetto a Jacobi, infatti il tasso asintotico di convergenza è $\log \rho (P)$ e si ha

${\frac {\rho (P_{g})}{\rho (P_{j})}}={\frac {\log \rho (P_{j})}{\log \rho (P_{j})^{2}}}=1/2$

(tuttavia non è detto che sia sempre più conveniente usare Gauss-Seidel).

Tasso di convergenza

Consideriamo la relazione:

$\vert E_{k}\vert \leq \vert P^{k}\vert *\vert E_{k-1}\vert$

Allora valgono le seguenti definizioni

${\begin{array}{cc}{\hbox{fattore di convergenza}}&=\vert P^{k}\vert \\\hline \\{\hbox{fattore di convergenza medio}}&=\vert P^{k}\vert ^{1/k}\\\hline \\{\hbox{tasso di convergenza medio}}&R_{k}(P)=-1/k\log \vert P^{k}\vert \\\hline \\{\hbox{tasso asintotico di convergenza}}&R(P)=\lim _{k\to \infty }r_{k}(P)=\log \rho (P)\end{array}}$

La definizione di tasso asintotico di convergenza non dipende né dalla norma né dal passo.

Test d'arresto test dell'incremento

Si considera la differenza fra un'iterata e la successiva, misurata con una certa norma, e si richiede che il metodo si arresti quando è verificata la condizione

$\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert \leq \varepsilon$

dove $\varepsilon$ è una tolleranza assegnata.

$E_{k}=\vert {\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k}\vert$

$\longrightarrow \,\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert =\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{\ast }\vert +\vert {\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k-1}\vert \leq E_{k}+E_{k-1}$

$=E_{k-1}+P*E_{k-1}=(I+P)*E_{k-1}$

Quindi

${\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\leq (I+P)E_{k-1}$

$\longrightarrow \,E_{k-1}=(I+P)^{-1}({\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1})$

Passando alla norma:

$\vert E_{k-1}\vert \leq \vert (I+P)^{-1}\vert *\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert \quad {\hbox{relazione }}\ast$

Applichiamo il teorema: se $\vert A\vert <1$ con norma matriciale indotta, allora $I+A$ è non singolare e

$\vert (I+A)^{-1}\vert \leq {\frac {1}{1-\vert A\vert }}$

Allora, sostituendo l'espressione trovata per $\vert (I+P)^{-1}\vert$ nella relazione $\ast$ ottengo:

$\vert E_{k-1}\vert \leq {\frac {1}{1-\vert P\vert }}*\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert$

Quindi

$\vert E_{k}\vert =\vert PE_{k-1}\vert \leq \vert P\vert *\vert E_{k-1}\vert$

$={\frac {\vert P\vert }{1-\vert P\vert }}\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert$

Il test è affidabile tranne nel caso in cui $\vert P\vert \simeq 1$.

Test d'arresto Test del residuo

Il residuo è la quantità ${\boldsymbol {b}}-A{\boldsymbol {x}}$. Si definisce il residuo al passo $k$ come ${\boldsymbol {b}}-A{\boldsymbol {x}}_{k}$. Scelta una norma, si impone che la norma del residuo sia minore di $\varepsilon$, e il metodo si arresta quando questa condizione è verificata.

Sapendo che

$A{\boldsymbol {x}}_{\ast }={\boldsymbol {b}}$

allora

$R_{k}={\boldsymbol {b}}-A{\boldsymbol {x}}_{k}=A*({\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k})=A*E_{k}\quad {\hbox{equazione dell'errore}}$

allora

$E_{k}=A^{-1}R_{k}$

$\vert E_{k}\vert =\vert A^{-1}R_{k}\vert \leq \vert A^{-1}\vert *\vert R_{k}\vert$

$={\frac {K(A)}{\vert A\vert }}*\vert R_{k}\vert$

Allora l'errore assoluto è controllato bene dal residuo, se il numero di condizionamento non è molto grande.

Il test del residuo è affidabile se il problema non è malcondizionato.

Se invece vogliamo stimare l'errore relativo:

${\frac {\vert E_{k}\vert }{\vert {\boldsymbol {x}}_{\ast }\vert }}$

$\leq {\frac {\vert A^{-1}\vert *\vert R_{k}\vert }{\vert {\boldsymbol {x}}_{\ast }\vert }}$

Allora, osserviamo che:

$\vert {\boldsymbol {b}}\vert \leq \vert A\vert *\vert {\boldsymbol {x}}_{\ast }\vert$

$\longrightarrow {\frac {1}{\vert {\boldsymbol {x}}_{\ast }\vert }}\leq {\frac {\vert A\vert }{\vert {\boldsymbol {b}}\vert }}$

allora

${\frac {\vert E_{k}\vert }{\vert {\boldsymbol {x}}_{\ast }\vert }}\leq {\frac {\vert A^{-1}\vert *\vert A\vert *\vert R_{k}\vert }{\vert {\boldsymbol {b}}\vert }}$

$=K(A)*{\frac {\vert R_{k}\vert }{\vert {\boldsymbol {b}}\vert }}$

La quantità ${\frac {\vert R_{k}\vert }{\vert {\boldsymbol {b}}\vert }}$ viene chiamata residuo relativo, se supponiamo $r_{0}={\boldsymbol {b}}$ se ${\boldsymbol {x}}_{0}=0$. L'errore relativo viene controllato dal residuo relativo a meno di condizionamento.

Il costo computazionale dell'iterazione è $o(n^{2})$ operazioni per iterazione, mentre se la matrice è sparsa il costo è $o(n)$.