Metodo di Jacobi

In analisi numerica il metodo di Jacobi è un metodo iterativo per la risoluzione di sistemi lineari, un metodo cioè che calcola la soluzione di un sistema di equazioni lineari dopo un numero teoricamente infinito di passi. Per calcolare tale risultato il metodo utilizza una successione $x^{(k)}$ che converge verso la soluzione esatta del sistema lineare e ne calcola progressivamente i valori arrestandosi quando la soluzione ottenuta è sufficientemente vicina a quella esatta.

Idea

Rappresentando il sistema lineare in forma matriciale con $Ax=b$, a partire da un $x^{(0)}$ iniziale dato, la successione $x^{(k+1)}=F(x^{(k)})$ con $k \in \N$ è definita come

$A=M-N$ dove M è una matrice non singolare (cioè con determinante non nullo e quindi invertibile)

$\begin{matrix}Ax=b \Leftrightarrow Mx=Nx+b \Leftrightarrow & x &=& M^{-1}Nx+M^{-1}b \\ & &=& F(x)\end{matrix}$
dove F è una funzione affine.

Algoritmo

$\left\{ \begin{matrix} x^{(0)} \\ x^{(k+1)} = M^{-1}Nx^{(k)}+M^{-1}b \end{matrix} \right.$
Se x è la soluzione di $Ax=b$ allora $x = M^{-1}Nx+M^{-1}b$ .

Stima dell'errore

Sia $e^{(k)}$ il vettore errore

$e^{(k+1)}=x^{(k+1)}-x^{(k)}=M^{-1}N(x^{(k)}-x^{(k-1)})=M^{-1}Ne^{(k)}$
Poniamo $B = M^{-1}N$, da cui $e^{(k+1)}=Be^{(k)}=B^{(k+1)}e^{(0)}$.

Convergenza

L'algoritmo converge se $\lim_{k \to \infty} \| e^{(k)} \| = 0 \Leftrightarrow \lim_{k \to \infty} \| B^k \| = 0$.

Vi sono inoltre due teoremi, il primo afferma che condizione necessaria e sufficiente affinché $\lim_{k \to \infty} \| B^k \| = 0$ è che il raggio spettrale (cioè il più grande autovalore in modulo) di B sia strettamente inferiore a 1; il secondo che il metodo converge per tutti i sistemi lineari con matrice a diagonale dominante $\forall x^{(0)}$.

Metodo di Jacobi

Si decompone la matrice A come segue: A = D-L-U dove D è la matrice diagonale, -L la parte triangolare inferiore e -U la parte triangolare superiore. Nel metodo di Jacobi si pongono M = D e N = L+U (mentre nel metodo di Gauss-Seidel, M = D-L e N = U).

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

Osservando che la i-esima riga della matrice $D^{-1}(L+U)$ si può scrivere come

$-(\frac{a_{i,1}}{a_{i,i}},\cdots, \frac{a_{i,i-1}}{a_{i,i}},0,\frac{a_{i,i+1}}{a_{i,i}},\cdots,\frac{a_{i,n}}{a_{i,i}})$

si ricava facilmente una espressione che rappresenta in termini più "algoritmici" il passo di iterazione del metodo:

$x^{(k+1)}_i= -\frac{1}{a_{i,i}} \sum_{j=1,j \ne i}^n a_{i,j}x^{(k)}_j + \frac{b_i}{a_{i,i}}$

Residuo

Sia $r^{(k)}=D e^{(k)}$ il vettore residuo, in cui $e$ è il vettore di errore (vedi sopra). Possiamo scrivere $x_i^{(k+1)}=\frac{r_i^{(k)}}{a_{i,i}} + x_i^{(k)}$ con $r_i^{(k)}$ che calcoliamo come segue: $r_i^{(k+1)}=-\sum_{j=1,j \ne i}^n a_{ij} \frac{r_i^{(k)}}{a_{j,j}}$.

Test d'arresto

Per il test d'arresto utilizziamo il vettore residuo. Data una tolleranza prefissata $\epsilon$, abbiamo: $\frac{\|r^{(k)} \|}{\|b\|} < \epsilon$

Conclusione

Questo metodo ha un costo computazionale dell'ordine di 3n²+2n per ogni iterazione, è facilmente parallelizzabile contrariamente al metodo di Gauss-Seidel che, tuttavia, converge più velocemente (nella maggior parte dei casi). Verificare che la complessità computazionale non sia: 2n²+2n!!

Bibliografia

• Alfio Quarteroni; Riccardo Sacco; Fausto Saleri, Risoluzione di sistemi lineari con metodi iterativi in Matematica numerica, 3^a edizione, Milano, Springer Verlag, gennaio 2008, Pagg.111-158. ISBN 978-88-470-0782-6

Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica