Metodi iterativi semplici

G. Puppo

Riassunto

Problema del fill-inMemorizzazione di matrici sparseMetodo di Jacobi

Problema del fill-in

Se risolvo il problema della membrana elastica usando un numero elevato di nodi interni, la soluzione diventa lentissima. Per esempio il comando:

>> u=membrana(f,100);

blocca il mio PC, anche se 100 nodi interni per ogni lato è una richiesta di dettaglio abbastanza frequente.

La function membrana risolve il sistema lineare che fornisce la soluzione con il comando:.

u=a\b';

Con questo comando, Matlab calcola la fattorizzazione LU della matrice A, e poi risolve due sistemi triangolari.Se la matrice A è M per M ed è piena, la fattorizzazione LU richiede circa M 3 operazioni.

Se costruiamo il problema della membrana con 100 nodi per lato, otteniamo una matrice circa 10000 per 10000. Se trattiamo la matrice A come matrice “piena”, la fattorizzazione LU richiede circa 10000 3 = 10 12 operazioni.

Se calcoliamo la fattorizzazione LU di una matrice A sparsa,cioè con un numero elevato di elementi nulli, otteniamo chei fattori L ed U sono molto più “pieni”

E’ possibile rendere il programma più efficiente utilizzando il fatto che la matrice A è sparsa, cioè la maggior parte degli elementi di A sono zero, e inoltre sappiamo dove questi zeri sono collocati.

function spyLa 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.Inoltre SPY stampa anche il numero di elementi diversi da zero.

Esempio.Studiamo la sparsità della fattorizzazione LU della matrice A della membrana elastica.Per far questo, applichiamo SPY sia ad A che alla sua fattorizzazione LU, per N=10 e per N=20. Per visualizzare sia L che U nella figura che segue, la function SPY è stata applicata alla matrice L+U.

N = 10

N = 20

La struttura a banda di A rimane anche nei fattori L ed U, ma all’interno della banda, la fattorizzazione LU crea elementi diversi da zero

Si verifica un comportamento molto simile anche per il problema della trave elastica:

N = 25

N = 50

Commenti

Matlab si accorge dalla struttura a banda della matrice A, e calcola solo gli elementi della fattorizzazione LU all’interno della banda. Questo riduce il costo del calcolo della fattorizzazione notevolmente: per la matrice della mebrana elastica, il costo diventa di ordine N 4, invece di N 6

Anche tenendo conto di questo fatto, se risolvo un sistema lineare sparso usando la fattorizzazione LU devo:- calcolare un elevato numero di elementi;- memorizzare tutti gli elementi calcolati.

Per risolvere i sitemi lineari sparsi, di grandi dimensioni, sono indicati i metodi iterativi

Implementazione dei metodi iterativi

Per implementare i metodi iterativi in modo efficace, è necessario sfruttare le tecniche del calcolo vettoriale e la struttura della matrice.

Nella tabella compaiono i tempi di CPU per tre diverse implementazioni del metodo di Jacobi per il problema della membrana elastica, con 4 valori diversi si N

Implementare il metodo di Jacobi per componenti significa applicare direttamente la formula:

for i=1:n sum=b(i); for j=1:n sum=sum-a(i,j)*x(j); end xnew(i,1)=x(i) + sum/a(i,i);end

Per sfruttare le proprietà vettoriali, uso il metodo di Jacobi in forma matriciale:

Ho bisogno di un’unica istruzione:

xnew = xold + dd\( b - a*xold );

dove dd è la matrice diagonale che contiene la diagonale principale di A

Matrici sparseMatlab memorizza le matrici sparse in un array che contiene gli elementi diversi da zero di A e il loro indirizzo, ordinati per colonne. Per esempio:

a = (1,1) 4 (2,1) -1 (5,1) -1 (1,2) -1 (2,2) 4 (3,2) -1 (6,2) -1 (2,3) -1 (3,3) 4

Contiene le prime colonne della matrice della membrana elastica memorizzata in forma sparsa

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 posizioneindicate dal vettore d:

Esempio:

>> 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);

Crea una matrice 10X10, con 5 diagonali non nulle

b = (1,1) 6 (2,1) -1 (6,1) 1 (1,2) -1 (2,2) 6 (3,2) -1 (7,2) 1 ………….

Inoltre, spdiags permette anche di estrarre diagonali da una matrice:

>>a=laplaciano(3);>>dd=spdiags(a,0)

dd = 4 4 4 4 4 4 4 4 4

Ottengo il vettore:

Se invece desidero la matrice diagonale che ha sulla diagonale gli elementi della diagonale di A:

>> n=3;>> a=laplaciano(n);>> dd=spdiags(a,0);>> dd=spdiags(dd,0,n^2,n^2)

dd = (1,1) 4 (2,2) 4 (3,3) 4 (4,4) 4 (5,5) 4 (6,6) 4 (7,7) 4 (8,8) 4 (9,9) 4

In questo modo ottengo la matrice diagonale in forma sparsa:

Esempio

function a=spmat_laplaciano(n)% A=SPMAT_LAPLACIANO(N) calcola la matrice del laplaciano% sul quadrato, con N nodi interni su ogni lato, in forma% sparsa

Scrivere una function che costruisca la matrice della membrana elastica in forma sparsa:

Le matrici sparse in Matlab sono costruite per colonne: è necessario prestare attenzione per impostare correttamente le diagonali non costanti

% crea le 5 diagonali della matrice Ab=ones(n^2,5);% diagonale principaleb(:,1)=4; d(1)=0;% sopradiagonaleb(:,2)=-1; b(1:n:n^2,2)=0; d(2)=1;% sottodiagonaleb(:,3)=-1; b(n:n:n^2,3)=0; d(3)=-1;% diagonali lontaneb(:,4)=-1; d(4)=n;b(:,5)=-1; d(5)=-n;a=spdiags(b,d,n^2,n^2);

Qui il primo zero va messo in posizione 1

Qui il primo zero va messo in posizione n

Esercizi

Costruire una function che calcoli la matrice di rigidità in forma sparsa per il problema di convezione diffusione con elementi P1

Costruire una function che calcoli la matrice di rigidità in forma sparsa per il problema della trave elastica

Per il primo esercizio, ottengo:

function a=spmat_cd2d(n,nu,beta)% A=SPMAT_CD2D(N,NU,[A,B]) calcola la matrice FEM % di convezione-diffusione sul quadrato, in forma sparsa,% con N nodi interni su ogni lato.% Nu e' la diffusione, BETA e' il vettore velocita'.betax=beta(1); betay=beta(2);h=1/(n+1); n2=n^2; %ampiezza di grigliapex=betax*h/(2*nu) %numero di Peclet lungo xpey=betay*h/(2*nu) %numero di Peclet lungo ypex=pex/3; pey=pey/3;

b=zeros(n2,7);% diagonale principaleb(:,1)=4; d(1)=0;% sopradiagonaleb(:,2)=-1+pex+2*pey; b(1:n:n2,2)=0; d(2)=1;% sottodiagonaleb(:,3)=-1-pex-2*pey; b(n:n:n2,3)=0; d(3)=-1;% sopra diagonali lontaneb(:,4)=-1+2*pex+pey; d(4)=n;b(:,5)=pex-pey; b(n:n:n2,5)=0; d(5)=n-1;% sotto diagonali lontaneb(:,6)=-1-2*pex-pey; d(6)=-n;b(:,7)=pey-pex; b(1:n:n2,7)=0; d(7)=-n+1;b=nu*b;a=spdiags(b,d,n^2,n^2);

Metodo di Jacobi

Posso ora scrivere una function che implementa il metodo di Jacobi utilizzando il calcolo vettoriale e le matrici sparse.

In input devo avere: i dati (la matrice A e il vettore B) e la tolleranza TOL , per costruire il test di arresto.

In output devo avere: la soluzione, il residuo ed ilnumero di iterazioni (per controllare la convergenza)

function jacobi.m

function [xnew,residuo,nit]=jacobi(a,b,tol)% X=JACOBI(A,B,TOL): Calcola la soluzione X del sistema% A*X=B, usando il metodo iterativo di Jacobi% con tolleranza TOL sul test di arresto% se TOL viene omessa, sceglie TOL=100*eps% [X,NIT]=JACOBI(A,B,TOL) Calcola la soluzione X del sistema% A*X=B e il numero NIT di iterazioni eseguite

if nargin < 3 tol = eps*100;end

% Estrae la diagonale principale di Add=spdiags(a,0);% Costruisce la matrice diagonale come matrice sparsan = length(a);dd = spdiags(dd,0,n,n);% Usa B come stima iniziale X0xold = b;% Stima un tetto al numero massimo di iterazioninmax=length(dd)^2;

Inizializzazione:

for n = 1:nmax res=b-a*xold; xnew = xold + dd\res; % Test di arresto res = norm(res); residuo(k)=res; diff = norm(xnew-xold); if res <= tol*norm(b) nit=n; return else xold = xnew; endendnit=nmax;

Iterazione principale e test di arresto:

Il metodo di Jacobi converge per la matrice del laplaciano.Il residuo diminuisce ad ogni iterazione secondo questo grafico:

Iterate successive per il metodo di Jacobi applicato al problema della membrana elastica, a partire dalla soluzione iniziale u ≡ 1

Il grafico mostra l'andamento del numero di iterazioni necessario per ottenere la convergenza (curva blu) e il numero di condizionamento della matrice (verde) in funzione di N

EserciziStudiare l'andamento del numero di iterazioni in funzione della tolleranza tol

Modificare il test di arresto: fermare l'iterazione se la differenza fra due iterazioni successive scende sotto la tolleranza specificata. Il comportamento del numero di iterazioni rispetto ad N cambia?

Costruire un problema di cui sia nota la soluzione esatta, per esempio, x=rand(n,1); b=a*x. Confrontare il residuo ottenuto a convergenza con l'errore

Script script_jac.m

% Lancia il metodo iterativo di Jacobi su una matrice assegnata n=20; % Numero nodi internia=spmat_laplaciano(n);b=ones(n^2,1);tol=1e-4;[x1,res,nit]=jacobi(a,b,tol);iterazioni=nitsemilogy(res,'Linewidth',2)

Problema della trave elasticaApplico il metodo di Jacobi al problema della trave elastica

In questo caso il numero di iterazioni cresce molto velocemente con le dimensioni della matrice, perchè il condizionamento della matrice cresce molto velocemente con N.

Infatti: Numero iterazioni Condizionamento N=10 11787 3.05e+3 N=20 148095 3.60e+5

Metodo di Richardson

Il metodo di Richardson è definito dalla seguente relazione:

Il metodo converge per:

mentre la convergenza ottimale si ha per:

function richardson.m

function [xnew,residuo,nit]=richardson(a,b,tol)% X=RICHARDSON(A,B,TOL): Calcola la soluzione X del sistema% A*X=B, usando il metodo iterativo di Richardson% con tolleranza TOL sul test di arresto% se TOL viene omessa, sceglie TOL=100*eps% [X,RESIDUON,NIT]=JACOBI(A,B,TOL) Calcola la soluzione % X del sistema% A*X=B, il residuo, e il numero NIT di iterazioni eseguite

if nargin < 3 tol = eps*100;end% Calcola il parametro alpha di convergenza ottimale, stimando% gli autovalori di modulo massimo e minimolmax=eigs(a,1,'lm');lmin=eigs(a,1,'sm');alpha=2/(lmax+lmin);% Usa B come stima iniziale X0xold = b; normb=norm(b);% Stima un tetto al numero massimo di iterazioninmax=length(b)^4;

Inizializzazione:

for n = 1:nmax res = b-a*xold; xnew = xold + alpha*res; % Test di arresto res=norm(res); residuo(n)=res; diff = norm(xnew-xold); if res <= tol*normb nit=n; return else xold = xnew; endendnit=nmax; display('non converge')

Iterazione principale:

Convergenza per il metodo di Richardson

Metodo di Gauss Seidel

Il metodo di Gauss Seidel è definito dalla relazione

m=tril(a); Per costruire la matrice M:

res = b-a*xold; xnew = xold + m\ res;

L'iterazione principale conterrà:

Confronto fra Jacobi e Gauss Seidel

Jacobi richiede circa il doppio di iterazioni rispetto a Gauss-Seidel

Metodo SOR

Il metodo SOR è definito dalla relazione:

Numero di iterazioni vs. omega

L'efficacia del metodo dipende fortemente da omega.

Esercizi

1) Usare i metodi iterativi appena descritti per il problema di convezione diffusione, variando il numero di Peclet. Che cosa si osserva per numeri di Peclet maggiori di 1?

3) Stimare il valore ottimale di omega per il problema della membrana elastica per diversi valori di N. Usare i risultati ottenuti per dare una “ricetta” per assegnare omega per il metodo SOR applicato al laplaciano.

2) Disegnare un grafico sulla convergenza del residuo per il metodo SOR con parametro ottimale