Metodi iterativi G Puppo Metodi iterativi Problema del

Metodi iterativi G. Puppo

Metodi iterativi • • Problema del fill-in Metodo di Jacobi Metodo di Gauss Seidel Studio della convergenza

Problema del Fill-in Se calcoliamo la fattorizzazione LU di una matrice A sparsa cioè con un numero elevato di elementi nulli, otteniamo che i fattori L ed U sono molto più “pieni”

Esempio Consideriamo per esempio la matrice: >> a=[4 -1 0 0 1; -1 4 -1 0 0; 0 -1 4 -1 0; 0 0 -1 4 -1; 1 0 0 -1 4] a= 4 -1 0 0 1 -1 4 -1 0 0 0 -1 4 -1 1 0 0 -1 4 Questa matrice è tridiagonale, con solo 2 elementi diversi da 0 fuori della struttura tridiagonale

Calcolando la fattorizzazione LU otteniamo >> [l, u]=lu(a) l= 1. 0000 0 0 -0. 2500 1. 0000 0 0 -0. 2667 1. 0000 0 0 -0. 2679 1. 0000 0 0. 2500 0. 0667 0. 0179 -0. 2632 1. 0000 u= 4. 0000 -1. 0000 0 0 1. 0000 0 3. 7500 -1. 0000 0 0. 2500 0 0 3. 7333 -1. 0000 0. 0667 0 0 0 3. 7321 -0. 9821 0 0 3. 4737

Continuo lo studio del fill-in Per approfondire lo studio del fill in, devo creare matrici sparse più grandi e analizzarne la struttura. Servono dei nuovi comandi • Function SPDIAGS • Function FULL • Function SPY

Function SPDIAGS Il comando A=spdiags(B, d, m, n) crea una matrice A m per n, con diagonali uguali alle colonne di B, disposte nella posizione indicate dal vettore d: b= Esempio: (1, 1) 6 (2, 1) -1 >> n=10; (6, 1) 1 >> e=ones(n, 1); (1, 2) -1 >> b=[e, -e, 6*e, -e, 2*e]; (2, 2) 6 >> d=[-n/2 -1 0 1 n/2]; (3, 2) -1 >> a=spdiags(b, d, n, n); (7, 2) 1 Crea una matrice 10 X 10, con 5 diagonali …………. non nulle

Function FULL La matrice creata da SPDIAGS è memorizzata considerando solo gli elementi diversi da zero, in modo da conservarne la struttura sparsa. Per avere la matrice completa, devo espanderla con FULL.

>> n=10; >> e=ones(n, 1); >> b=[e, -e, 6*e, -e, 2*e]; >> d=[-n/2 -1 0 1 n/2]; >> a=spdiags(b, d, n, n); >> full(a) ans = 6 -1 0 0 0 2 0 0 0 0 -1 6 -1 0 0 0 2 1 0 0 0 -1 6 -1 0 0 0 0 1 0 0 0 -1 6

Function SPY La function spy(a) permette di visualizzare la sparsità di una matrice. Il comando SPY(A) genera un grafico, nel quale sono evidenziati con un punto solo gli elementi di A che sono diversi da zero. Esempio. Studiamo la sparsità della fattorizzazione LU della matrice A. Per far questo, applichiamo SPY sia ad A che alla fattorizzazione LU, per N=10 e per N=100. Per visualizzare sia L che U nella figura che segue, la function SPY è stata applicata alla matrice L+U.

Risultati ottenuti con il programma fill_in. m

Listato dello script fill_in. m % Questo programma studia il fill-in della fattorizzazione LU % di una matrice sparsa A k=0; for n=[10, 100] e=ones(n, 1); b=[e, -e, 6*e, -e, 2*e]; d=[-n/2 -1 0 1 n/2]; a=spdiags(b, d, n, n); full(a); k=k+1; subplot(2, 2, k) spy(a) title('Matrice A') k=k+1; [l, u]=lu(a); subplot(2, 2, k) spy(l+u) title('Fattorizzazione LU') end

Commenti I risultati precedenti dimostrano che la fattorizzazione LU di una matrice sparsa genera un gran numero di elementi diversi da zero. Il numero degli elementi diversi da zero inoltre cresce velocemente all’aumentare delle dimensioni della matrice. Quindi, se risolvo un sistema lineare sparso usando la fattorizzazione LU devo: - calcolare un elevato numero di elementi; - memorizzare tutti gli elementi calcolati. Per questo tipo di sistemi, è conveniente usare i metodi iterativi

Metodi iterativi Per applicare un metodo iterativo ad una matrice sparsa, devo evitare di memorizzare tutta la matrice, altrimenti perdo tutti i possibili vantaggi del metodo iterativo. Infatti, se ho abbastanza memoria per memorizzare A, ho anche abbastanza memoria per memorizzare la fattorizzazione LU. Inoltre, ogni iterazione richiede il prodotto A*x, e diventa quindi molto costosa se non sfrutto la sparsità di A.

Per applicare un metodo iterativo, quindi, devo sfruttare la struttura del sistema lineare che devo risolvere. In particolare, devo usare matrici sparse. Tuttavia, iniziamo per ora ad applicare i metodi di Jacobi e di Gauss Seidel usando la formula generale. . Le functions che otterremo ora non sfrutteranno la struttura di A e quindi non saranno efficienti, ma ci daranno informazioni sulla velocità di convergenza.

Metodo di Jacobi Scriviamo una function che applichi il metodo di Jacobi ad un generico sistema A*x=b. La function richiesta deve: • dare in output il vettore soluzione e il numero delle iterazioni che sono state eseguite; • avere in input la matrice A ed il termine noto B; • contenere un opportuno test di arresto.

Listato per il metodo di Jacobi (function jacobi. m) function [xnew, nit]=jacobi(a, b) % JACOBI(A, B) calcola la soluzione XNEW ottenuta con il metodo di % Jacobi e il numero NIT di iterazioni necessarie % per il sistema lineare A*XNEW=B % Sintassi: [XNEW, NIT]=JACOBI(A, B) % Attenzione: Applicazione naive del metodo di Jacobi, che non % sfrutta la sparsita' di A. [n, m]=size(a); if m ~= n display('A non e'' quadrata') return end m=length(b); if m ~= n display('B non e'' compatibile') return end continua. . .

Calcolo della nuova stima: % Come vettore iniziale usa b: x=b’; kmax=n^2; for k=1: kmax for i=1: n sum=b(i); for j=1: n if j~=i sum=sum-a(i, j)*x(j); end xnew(i)=sum/a(i, i); end continua. . . La formula iterativa per il metodo di Jacobi è

Test di arresto: % Test di arresto: res = norm(a*xnew' -b); diff = norm(x-xnew); if res <= eps*norm(b) | diff<=eps*norm(x) nit=k; return else x=xnew; end nit=kmax;

Risolvo, usando il metodo di Jacobi, il sistema lineare A*X=B, dove A è la stessa matrice N X N costruita dalla function SPARSA(N) e B è il vettore ONES(N, 1): function a=sparsa(n) % Genera la matrice sparsa n*n a diagonale dominante % usata come esempio in questo capitolo e=ones(n, 1); % per avere una matrice a diagonale dominante, diag>=5 diag=6; b=[e, -e, diag*e, -e, 2*e]; d=[-n/2, -1, 0, 1, n/2]; a=spdiags(b, d, n, n);

Per applicare il metodo di Jacobi, devo dare i seguenti comandi: Considero N=10 >> n=10; >> a=sparsa(n); >> afull=full(a); >> b=ones(n, 1); >> [x, nit]=jacobi(afull, b); Ottengo: >> x x= Columns 1 through 8 0. 1197 0. 1393 0. 1440 Columns 9 through 10 0. 2059 0. 1728 >> nit = 54 0. 1503 0. 1692 0. 2106 0. 2140 0. 2126

Esercizio Modificare la function per il metodo di Jacobi, facendo in modo che la function immagazzini in un vettore res(k) la norma del residuo all’iterazione k: res(k)= norm(b - A*xk), e aggiungere questo vettore alle variabili in output. Fare un grafico del residuo in funzione di k, per il metodo di Jacobi applicato al sistema Ax=b, dove A è la matrice costruita dalla function tridiag(n), per N=10, N=20, N=40, e b è il vettore b = 1.

Metodo di Gauss Seidel Scriviamo una function che applichi il metodo di Gauss Seidel ad un generico sistema A*x=b. La function richiesta deve: • dare in output il vettore soluzione e il numero delle iterazioni che sono state eseguite; • avere in input la matrice A ed il termine noto B; • contenere un opportuno test di arresto.

Listato per il metodo di Gauss Seidel (function gs. m) function [x, nit]=gs(a, b) % GS(A, B) calcola la soluzione XNEW ottenuta con il metodo di % Gauss Seidel e il numero NIT di iterazioni necessarie % per il sistema lineare A*XNEW=B % Sintassi: [XNEW, NIT]=GS(A, B) % Attenzione: Applicazione naive del metodo di Gauss Seidel, che non % sfrutta la sparsita' di A. [n, m]=size(a); if m ~= n display('A non e'' quadrata') return end m=length(b); if m ~= n display('B non e'' compatibile') return end

Calcolo della nuova stima: % Come vettore iniziale usa b: x=b'; kmax=n^2; for k=1: kmax xold=x; %immagazzina il vecchio vettore X for i=1: n sum=b(i); for j=1: n if j~=i sum=sum-a(i, j)*x(j); end x(i)=sum/a(i, i); %Riscrive su X end

Test di arresto: % Test di arresto: res = norm(a*x' -b); diff = norm(x-xold); if res <= eps*norm(b) | diff<=eps*norm(x) nit=k; return end nit=kmax;

Per applicare il metodo di Gauss-Seidel, devo dare i seguenti comandi: >> n=10; >> a=sparsa(n); >> afull=full(a); >> b=ones(n, 1); >> [x, nit]=gs(afull, b); Ottengo: >> x x= Columns 1 through 8 0. 1197 0. 1393 0. 1440 Columns 9 through 10 0. 2059 0. 1728 >> nit = 30 0. 1503 0. 1692 0. 2106 0. 2140 0. 2126

Commenti Ottengo circa lo stesso vettore soluzione, X, ma il numero di iterazioni è molto più basso (circa la metà). Si può dimostrare infatti che per matrici a diagonale dominante: • convergono sia il metodo di Jacobi che il metodo di Gauss Seidel. • Il metodo di Gauss Seidel richiede circa la metà delle iterazioni effettuate dal metodo di Jacobi.

Esercizio Ripetere l’esercizio precedente, costruendo questa volta un grafico per il residuo ottenuto con il metodo di Gauss Seidel, in funzione di k. Confrontare i risultati del metodo di Gauss Seidel con quelli ottenuti con il metodo di Jacobi

Metodo efficiente di Jacobi Tutti i programmi scritti finora non utilizzano le funzionalità vettoriali di Matlab. E’ possibile ottenere una versione molto più efficiente. Per ottenere qualcosa di meglio, devo cercare di evitare di scrivere i conti per componenti. L’algoritmo del metodo di Jacobi in forma vettoriale è: D*x n+1 = -(A-D)*x n + b Quindi devo estrarre la diagonale D di A e risolvere il sistema: x n+1 = -D ((A-D) * x n + b) Notare che in questo modo A può essere una matrice memorizzata in forma sparsa

Nuova versione della function per il metodo di Jacobi: function [xnew, nit]=jacobi(a, b) % X=JACOBI(A, B): Calcola la soluzione X del sistema % A*X=B, usando il metodo iterativo di Jacobi % [X, NIT]=JACOBI(A, B) Calcola la soluzione X del sistema % A*X=B e il numero NIT di iterazioni eseguite % Estrae la diagonale principale di A dd=diag(a, 0); % Costruisce la matrice diagonale come matrice sparsa n = length(a); dd = spdiags(dd, 0, n, n); % Usa B come stima iniziale X 0 xold = b; % Stima un tetto al numero massimo di iterazioni nmax=length(dd)^2; La diagonale dd deve essere scritta come matrice diagonale per poter essere sommata ad A

…continua. . . for n = 1: nmax xnew = dd( (dd-a)*xold +b); % Test di arresto res = norm(a*xnew -b); diff = norm(xnew-xold); if res <= eps*norm(b) | diff<=eps*norm(xold) nit=n; return else xold = xnew; end nit=nmax; N. B. questa function può essere usata assegnando in input sia una matrice A scritta in forma sparsa, che una matrice A piena: la velocità di esecuzione cambia drasticamente

Questi sono i tempi di esecuzione ottenuti con matrici piene…

…e questi sono i tempi di esecuzione ottenuti con le stesse matrici in forma sparsa Se vi sembra che i due grafici siano simili, confrontate le scale verticali

Convergenza dei metodi iterativi Sappiamo che un metodo iterativo converge se e solo se il raggio spettrale della matrice di iterazione è minore di 1. Quindi un metodo per stabilire la convergenza di un metodo iterativo è il seguente: • Calcolo la matrice di iterazione. • Calcolo gli autovalori della matrice di iterazione. • Prendo l’autovalore di modulo massimo e ne studio il modulo

Autovalori di una matrice Per calcolare gli autovalori di una matrice, Matlab dispone della function EIG: >>x=eig(a) crea un vettore x che contiene una stima degli autovalori di a. >> [x, d]=eig(a) Crea una matrice X che contiene gli autovettori di A e una matrice D diagonale che contiene gli autovalori

Esempio: >> a=[1 2 3; 4 5 6; 7 8 9]; >> x=eig(a) x= 16. 1168 -1. 1168 -0. 0000 Verifico che gli autovalori trovati sono una stima degli autovalori esatti, calcolando il determinante di A - l*eye(3): >> for i=1: 3 res(i)=det(a-x(i)*eye(3)); end >> res = 1. 0 e-011 * 0. 1848 0. 0002 -0. 0018

Convergenza del metodo di Jacobi La matrice di iterazione per il metodo di Jacobi è B = D-1*(A-D), dove D contiene gli elementi sulla diagonale di A Costruisco un programma che calcoli il raggio spettrale della matrice di iterazione B per il metodo di Jacobi.

Function RHO=CONV_JAC(A) function rho=conv_jac(a) % Calcola il raggio spettrale RHO per la matrice di iterazione % del metodo di Jacobi applicato alla matrice A % Sintassi RHO=CONV_JAC(A) [n, m]=size(a); if m ~= n display('A non e'' quadrata') return end for i=1: n d(i, i)=a(i, i); end b=d(a-d); x=eig(b); rho=max( abs(x));

Convergenza del metodo di Gauss-Seidel La matrice di iterazione per il metodo di Gauss-Seidel è B = E-1*(A-E), dove E è formata dagli elementi della parte triangolare inferiore di A Costruisco un programma che calcoli il raggio spettrale della matrice di iterazione B per il metodo di Gauss-Seidel.

Function CONV_GS(A) function rho=conv_gs(a) % Calcola il raggio spettrale RHO per la matrice di iterazione % del metodo di Jacobi applicato alla matrice A % Sintassi RHO=CONV_JAC(A) [n, m]=size(a); if m ~= n display('A non e'' quadrata') return end for i=1: n for j=1: i d(i, j)=a(i, j); end b=d(a-d); x=eig(b); rho=max( abs(x));

Esercizio Scrivere una function che calcoli una matrice 5 per 5 con tutti gli elementi ai, j = -1, tranne che sulla diagonale principale, dove ai, j = 4 + δ. Studiare l’andamento del raggio spettrale della matrice di iterazione del metodo di Jacobi per valori piccoli di δ, per esempio δ compreso fra 0 e 0. 1. Per quali valori di δ il metodo di Jacobi converge più velocemente?