Calcolo numerico/Metodi iterativi

Descrizione e osservazioni introduttive

Nei metodi iterativi si crea una successione di vettori ${\boldsymbol {x}}_{k}=\phi {\boldsymbol {x}}_{k-1}+{\boldsymbol {c}}$, avendo assegnato un vettore ${\boldsymbol {x}}_{0}$ iniziale, e una matrice $\phi$. Il metodo viene iterato finché ${\boldsymbol {x}}_{k}$ è un'approssimazione sufficientemente buona della soluzione.

I metodi iterativi sono convenienti quando $A$ è una matrice sparsa, cioè quando il numero degli elementi diversi da 0 è trascurabile.

Applicando un metodo diretto a matrici di questo tipo, ho il fenomeno del fill in.

I metodi iterativi non alterano la struttura della matrice, e la loro operazione di base è quella di un prodotto matrice per vettore, che costa $o(n^{2})$ operazioni se la matrice è piena, e $o(n)$ se la matrice è sparsa (vantaggio computazionale).

Per sua natura un metodo iterativo è meno sensibile all'errore algoritmico di un metodo diretto. Supponiamo che le operazioni iterative convergano. A causa dell'errore algoritmico la ${\bar {\boldsymbol {x}}}_{\ast }$ può variare, ma questo non importa perché con questi metodi non si cerca di calcolare una soluzione esatta, ma solo un'approssimazione.

Il fatto di troncare la successione causa la presenza di un errore analitico, che era assente nei metodi diretti.

Si dice che ${\boldsymbol {x}}_{\ast }$ converge a ${\boldsymbol {x}}_{k}$ se questo vale puntualmente, cioè se

${\boldsymbol {x}}_{\ast ,i}=\lim _{k\to \infty }{\boldsymbol {x}}_{k,i}$

Esempio 2.9 (metodo di Jacobi)

Consideriamo la $i$-esima equazione di un sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$, cioè

$\sum _{j=1}^{n}A_{ij}x_{j}=b_{i}$

ed esplicitando rispetto a $x_{i}$ otteniamo

$x_{i}={\frac {b_{i}-\sum _{i\neq j}A_{ij}x_{j}}{a_{ii}}}$

e al passo successivo otteniamo:

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

In forma compatta, questo metodo si scrive come:

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

dove

$D=\mathrm {diag} (A_{1})$

$L=\mathrm {tril} (A_{1})\quad {\hbox{triangolare dalla prima sottodiagonale a scendere}}$

$U=\mathrm {triu} (A_{1})$

Si ha quindi una scrittura della forma

${\boldsymbol {x}}_{k+1}=\phi {\boldsymbol {x}}_{k}+{\boldsymbol {c}}$

Generalità sulle tecniche di splitting

Dato il sistema lineare $A{\boldsymbol {x}}={\boldsymbol {b}}$, $A\in \mathbb {C} _{n,n}$, chiamiamo ${\boldsymbol {x}}_{\ast }$ la soluzione esatta. Consideriamo uno splitting della matrice $A$ che viene scritta come $M-N$, dove $M,N$ sono matrici quadrate e $M$ è invertibile.

Applichiamo lo splitting alla matrice $A$ nel sistema lineare, cioè

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

e posso riscrivere

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

allora poniamo $P=M^{-1}N$ e ${\boldsymbol {c}}=M^{-1}{\boldsymbol {b}}$ quindi abbiamo il sistema

${\boldsymbol {x}}_{k}=P{\boldsymbol {x}}+{\boldsymbol {c}}$

Il metodo iterativo è

${\boldsymbol {x}}_{k+1}=P{\boldsymbol {x}}_{k}+{\boldsymbol {c}}\quad {\hbox{relazione}}\,\triangle$

e $P$ viene chiamata matrice di iterazione.

Supponiamo che la successione dei vetori ${\boldsymbol {x}}_{k}$ converga, allora

$\lim _{k\to \infty }{\boldsymbol {x}}_{k}={\boldsymbol {x}}_{\ast }$

Passando al limite nella relazione $\triangle$:

${\boldsymbol {x}}_{\ast }=P{\boldsymbol {x}}_{\ast }+{\boldsymbol {c}},\quad {\hbox{relazione }}\circ$

Sottraendo membro a membro le relazioni $\triangle$ e $\circ$ otteniamo:

${\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k+1}=P{\boldsymbol {x}}_{\ast }+{\boldsymbol {c}}-(P{\boldsymbol {x}}_{k}+{\boldsymbol {c}})=P({\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k})$

dove ${\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{k+1}$ è l'errore assoluto $E_{k+1}$ al passo $k+1$-esimo, e si ha

$E_{k+1}=PE_{k}=P^{2}E_{k-1}=\dots =P^{k+1}E_{0}$

Per ogni $k$,

$E_{k}=P^{k}E_{0}$

Convergenza del metodo iterativo

Definizione

Diciamo che un metodo iterativo è convergente se

$\lim _{k\to \infty }{\boldsymbol {x}}_{k}={\boldsymbol {x}}_{\ast }$

per ogni vettore d'innesco ${\boldsymbol {x}}_{0}$.

Il metodo iterativo $\triangle$ è convergente se e solo se

$\lim _{k\to \infty }E_{k}=0$

per ogni $E_{0}$, e quindi se

$\lim _{k\to \infty }P^{k}=0$

Teorema 2.12 (condizione necessaria e sufficiente per la convergenza)

Il metodo è convergente, cioè

$\lim _{k\to \infty }P^{k}={\hbox{matrice nulla}}$

se e solo se $\rho (P)\leq 1$.

Dimostrazione

$1\longrightarrow 2$: supponiamo che il metodo sia convergente, allora

$\lim _{k\to \infty }E_{k}=0$

Scegliamo ${\boldsymbol {x}}_{0}$ in modo tale che $E_{k}={\boldsymbol {x}}_{\ast }-{\boldsymbol {x}}_{0}$ sia un autovettore rispetto a $\lambda$ che realizza il raggio spettrale, quindi

$E_{k}=P^{k}E_{0}=\lambda ^{k}E_{0}$

(infatti $E_{0}$ è autovettore)

$E_{k}=\vert P^{k}E_{0}\vert =|\lambda |^{k}|\vert E_{0}\vert =\rho (P)^{k}\vert E_{0}\vert$

allora

$\lim _{k\to \infty }\rho (P)^{k}=0$

$|\rho (P)|<1$.

$2\longrightarrow 1$:

Caso di $P$ diagonalizzabile: $P$ può essere scritta come $P=S\lambda S^{-1}$ dove $\lambda$ è la matrice diagonale con entrate uguali agli autovalori (la maggior parte delle matrici sono diagonalizzabili).

Allora

$P^{k}=(S\lambda S^{-1})*(S\lambda S^{-1})(S\lambda S^{-1})\dots (S\lambda S^{-1})=S\lambda ^{k}S^{-1}$

e $P^{k}$ tende a 0 se $\lambda$ tende a 0, cioè se e solo se $|\lambda _{i}|\leq 1\forall i$, quindi se e solo se $\rho (P)\leq 1$.

Caso generale: per dimostrare questo caso serve la nozione di forma normale di Shur.

Teorema 2.13 (forma normale di Shur)

Data una matrice $A\in \mathbb {C} _{n,n}$ con autovalori $\lambda _{1},\dots ,\lambda _{n}$, allora esiste $Q$ unitaria e una matrice $T$ triangolare superiore con elementi diagonali $t_{ii}=\lambda _{i}$ tali che $A=QTQ_{h}$.

Se $A$ è hermitiana, $T$ diventa la matrice $\lambda$ degli autovalori.

Consideriamo la forma di Shur di $P$, cioè

$P=QTQ_{h}$

Allora consideriamo

$P^{k}=(QTQ_{h})(QTQ_{h})\dots (QTQ_{h})$

e siccome $Q$ è unitria, $QQ_{h}=I$, quindi rimane

$P^{k}=QT^{k}Q_{h}$

allora

$\lim _{k\to \infty }P^{k}=0$

se e solo se

$\lim _{k\to \infty }T^{k}=0$

Inoltre, $\rho (P)=\rho (T)$ perché le due matrici sono simili, e si ha

$0\leq |T^{k}|\leq |T|^{k}\leq (|T+E|)^{k}=|F|^{k}$

Definisco la matrice $E$ in modo che:

1. $E$ sia una matrice diagonale, non negativa
2. $|F|=|T+E|\leq 1$,
3. gli elementi diagonali di $T+E$ sono a due a due distinti, cioè richiedo che la matrice $T+E$ sia diagonalizzabile.

Allora, per $F$ diagonalizzabile, ho già dimostrato che

$\lim _{k\to \infty }|F|^{k}=0$

se e solo se

$\rho (|F|^{k})\leq 1$

Allora per confronto

$\lim _{k\to \infty }T^{k}=0$

e segue la tesi.

cvd

Supponiamo di considerare $|T|$ e di ordinarne gli autovalori per modulo decrescente, in modo che $\tau _{1}\geq \tau _{2}\geq \tau _{3}$. Allora la matrice $E$ è la matrice diagonale con entrate $\varepsilon _{1},\dots ,\varepsilon _{t}$ dove

$\varepsilon _{1}=(1-\tau _{1})/2$

$\varepsilon _{2}=\varepsilon _{1}/2$

$\varepsilon _{t}=\varepsilon _{t-1}/2$

Con questa tecnica si riesce a garantire che gli autovalori di $T+E$ sono a due a due distinti.

Altre osservazioni sulla convergenza

Vale la seguente proposizione:

Proposizione 2.4 [condizione sufficiente di convergenza]

Se esiste una norma matriciale indotta tale $\vert P\vert \leq 1$ allora il metodo è convergente.

Dim. Per ogni norma matriciale indotta vale che $\rho (P)\leq \vert P\vert$, allora $\rho (P)\leq 1$.

cvd

Osservazione 2.6

Supponiamo che esista una norma tale che

$\vert P\vert \leq c\leq 1$

(cioè il metodo converge)

Allora, usando la stessa norma, l'errore decresce in modo monotono.

Dim. Consideriamo

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

Inoltre sappiamo che

$\vert {\boldsymbol {x}}_{k}-{\boldsymbol {x}}_{k-1}\vert =\vert P*{\boldsymbol {x}}_{k-1}+{\boldsymbol {c}}-(P*{\boldsymbol {x}}_{k-2}+{\boldsymbol {c}})\vert =$

$=\vert P({\boldsymbol {x}}_{k-1}-{\boldsymbol {x}}_{k-2})\vert \leq \vert P\vert *\vert {\boldsymbol {x}}_{k-1}-{\boldsymbol {x}}_{k-2}\vert$

$\leq c*\vert {\boldsymbol {x}}_{k-1}-{\boldsymbol {x}}_{k-2}\vert$

e si ha una decrescita monotona anche dell'incremento.

cvd

Tuttavia $\rho (P)\leq 1$ non implica che $\vert P\vert \leq 1$ per qualsiasi norma, e l'errore prima di convergere può avere una crescita iniziale.

Esempio 2.10

Consideriamo una matrice d'iterazione

$\left({\begin{array}{cc}0&2\\0&0\end{array}}\right)$

$\rho (P)=0$

invece $\vert P\vert _{2}=2$, infatti

$P_{h}=\left({\begin{array}{cc}0&0\\2&0\end{array}}\right)$

$P_{h}P=\left({\begin{array}{cc}4&0\\0&0\end{array}}\right)$

e il polinomio caratteristico è

$p_{\lambda }=\lambda ^{2}-4\lambda$

quindi ${\sqrt {\max |\rho (A_{h}A)|}}=2$.

$\vert E_{0}\vert _{2}=0,\,\vert E_{1}\vert _{2}=2,\,\vert E_{2}\vert _{2}=0$

Quindi, prima che il metodo converga, l'errore ha una crescita iniziale.

Stiamo considerando metodi che generano la successione

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

tali che $\rho (P)<1$. In uno spazio finitodimensionale, non abbiamo la dipendenza dalla norma.

Infatti

$T({\boldsymbol {x}})=P{\boldsymbol {x}}+{\boldsymbol {c}}$

$T({\boldsymbol {x}})\to {\boldsymbol {x}}_{\ast }$ se e solo se $\rho (P)<1$.

$\rho (P)$ è l'inf delle norme matriciali indotte della norma di $P$. Questo implica che, se $\rho (P)<1$, allora esiste una norma matriciale indotta $\ast$ tale che

$\vert P\vert _{\ast }\leq 1$

quindi

$\vert {\boldsymbol {x}}-{\boldsymbol {y}}\vert =\vert P{\boldsymbol {x}}+{\boldsymbol {c}}-(P{\boldsymbol {y}}+{\boldsymbol {c}})\vert =\vert P({\boldsymbol {x}}-{\boldsymbol {y}})\vert _{\ast }$

$\leq \vert P\vert _{\ast }*\vert {\boldsymbol {x}}-{\boldsymbol {y}}\vert _{\ast }$

allora per il teorema di Banach abbiamo convergenza.

Siccome tutte le norme in uno spazio finitodimensionale sono topologicamente equivalenti, si ha convergenza rispetto a una qualsiasi norma.

Se $\vert E_{k}\vert _{\ast }\to 0$, siccome $\vert E\vert =K^{\ast }\vert E_{k}\vert _{\ast }$, anche la norma a primo membro tende a 0.