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𝒙=(M-N)𝒙=𝒃$

$⟶𝒙=M^{-1}N𝒙+M^{-1}𝒃$

$A=M^{-1}N,𝒄=M^{-1}𝒃$

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

$P=-D^{-1}(L+U),𝒄=D^{-1}𝒃$

Indicando con $P_j$ la matrice d'iterazione del metodo di Jacobi si ha:

${𝒙}_k=P_j{𝒙}_{k-1}+{𝒄}_j$

e moltiplicando per $D$:

$D{𝒙}_k=D{P}_j{𝒙}_{k-1}+D{𝒄}_j$

e sostituendo le espressioni ricavate prima:

$D{𝒙}_k=-(L+U){𝒙}_{k-1}+𝒃$

Considerando la $i$-esima equazione si ha:

$x_i=\frac{-\sum _{j\ne 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,𝒄=(D+L{)}^{-1}𝒃$

Come prima, moltiplichiamo per $M$ la relazione

${𝒙}_k=P{𝒙}_{k-1}+𝒃$

e otteniamo

$(D+L){𝒙}_k=-U{𝒙}_{k-1}+𝒃$

cioè

$D{𝒙}_k=-L{𝒙}_k-U{𝒙}_{k-1}+𝒃$

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

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

$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{array}{c}{A}_{11}{x}_1+A_{12}{x}_2=b_1\\ A_{22}{x}_2=b_2\end{array}$

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}\ne 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

$|A{|}^{\mathrm{\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:

$|P|\le |D^{-1}|\ast |L+U|$

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

$|L+U|\le \underset{i}{max}|A_{ii}|$

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

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

quindi

$|P|\le \frac{1}{\underset{i}{max}|A_{ii}|}\ast \underset{i}{max}|A_{ii}|\le 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 {ℂ}_{n,n}$ tridiagonale, siano

$𝒂=(a_1,\dots ,a_n)\text{ diagonale principale}$

$𝒄=(c_1,\dots ,c_{n-1})\text{ sopradiagonale}$

$𝒃=(b_2,\dots ,b_n)\text{ 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 \ne 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:

$|E_k|\le |P^k|\ast |E_{k-1}|$

Allora valgono le seguenti definizioni

$\begin{array}{cc}\text{fattore di convergenza}& =|P^k|\\ \\ \text{fattore di convergenza medio}& =|P^k{|}^{1/k}\\ \\ \text{tasso di convergenza medio}& R_k(P)=-1/k\log |P^k|\\ \\ \text{tasso asintotico di convergenza}& R(P)=\underset{k\to \mathrm{\infty }}{lim}{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

$|{𝒙}_k-{𝒙}_{k-1}|\le \epsilon$

dove $\epsilon$ è una tolleranza assegnata.

$E_k=|{𝒙}_{\ast }-{𝒙}_k|$

$⟶|{𝒙}_k-{𝒙}_{k-1}|=|{𝒙}_k-{𝒙}_{\ast }|+|{𝒙}_{\ast }-{𝒙}_{k-1}|\le E_k+E_{k-1}$

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

Quindi

${𝒙}_k-{𝒙}_{k-1}\le (I+P)E_{k-1}$

$⟶E_{k-1}=(I+P{)}^{-1}({𝒙}_k-{𝒙}_{k-1})$

Passando alla norma:

$|E_{k-1}|\le |(I+P{)}^{-1}|\ast |{𝒙}_k-{𝒙}_{k-1}|\text{relazione }\ast$

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

$|(I+A{)}^{-1}|\le \frac{1}{1-|A|}$

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

$|E_{k-1}|\le \frac{1}{1-|P|}\ast |{𝒙}_k-{𝒙}_{k-1}|$

Quindi

$|E_k|=|P{E}_{k-1}|\le |P|\ast |E_{k-1}|$

$=\frac{|P|}{1-|P|}|{𝒙}_k-{𝒙}_{k-1}|$

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

Test d'arresto Test del residuo

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

Sapendo che

$A{𝒙}_{\ast }=𝒃$

allora

$R_k=𝒃-A{𝒙}_k=A\ast ({𝒙}_{\ast }-{𝒙}_k)=A\ast E_k\text{equazione dell'errore}$

allora

$E_k=A^{-1}{R}_k$

$|E_k|=|A^{-1}{R}_k|\le |A^{-1}|\ast |R_k|$

$=\frac{K(A)}{|A|}\ast |R_k|$

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{|E_k|}{|{𝒙}_{\ast }|}$

$\le \frac{|A^{-1}|\ast |R_k|}{|{𝒙}_{\ast }|}$

Allora, osserviamo che:

$|𝒃|\le |A|\ast |{𝒙}_{\ast }|$

$⟶\frac{1}{|{𝒙}_{\ast }|}\le \frac{|A|}{|𝒃|}$

allora

$\frac{|E_k|}{|{𝒙}_{\ast }|}\le \frac{|A^{-1}|\ast |A|\ast |R_k|}{|𝒃|}$

$=K(A)\ast \frac{|R_k|}{|𝒃|}$

La quantità $\frac{|R_k|}{|𝒃|}$ viene chiamata residuo relativo, se supponiamo $r_0=𝒃$ se ${𝒙}_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)$.