Mathematical Engineering - Matematica Numerica

Exercise - 08 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 8 - Soluzione Metodi iterativi per sistemi lineari (2) Esercizio 1 Si consideri la seguente matrice : A =2 6 6 410 1 2 3 1 101 2 2 3 201 3 2 1 203 7 7 5 e siaAx=bil sistema lineare associato, dove il termine noto è scelto in modo tale che la soluzione esatta siax ex= [1 ;1;1;1]T . a)Vericare se la matriceAè simmetrica, denita positiva e/o a dominanza diagonale stretta. Stabilire, inoltre, se i metodi di Jacobi e Gauss-Seidel applicati al sistema convergano. b)Si calcoli la norma 2 della matrice di iterazioneB Jassociata al metodo di Jacobi, per il sistema dato. Stimare il numero minimo di iterazionikminnecessarie anchè il metodo di Jacobi riduca l'errore iniziale in norma 2 di un fattore"= 10 10 . c)Vericare numericamente il risultato ottenuto al punto precedente, implementando la funzionemyjacobicon rma[xn,iter]=myjacobi(A,b,x0,nmax,toll)e partendo dal vettore iniziale x0 = zeros(4,1). La funzione deve arrestarsi qualora sia raggiunto il numero massimo di iterazioni Nmaxo sia vericato il criterio di arresto : kr( k) k2k bk 2< toll; (1) essendor( k) =bAx( k) ilresiduoassociato allak-esima iterazione. d)Si consideri ora il metodo SOR, x( k+1) i=!a ii0 @bii 1 X j=1a ijx( k+1) jn X j=i+1a ijx( k) j1 A+ (1!)x( k) i; i = 1; : : : ; n; k0:(2) Modicare la funzionemyjacobiper ottenere una funzionemysorche implementi il metodo SOR per la risoluzione del generico sistemaAx=b. Si utilizzi poi la funzione per risolvere lo stesso problema con il metodo di Gauss-Seidel, con una tolleranzatoll= 10 10 : quante iterazioni vengono eettuate ? L'errore relativo è minore del fattore"= 10 10 considerato al punto b) ? Soluzione 1a)Verichiamo i criteri richiesti e che il raggio spettrale delle corrispondenti matrici di iterazione è minoredi uno.1 A = [10 -1 2 3; 1 10 -1 2; 2 3 20 -1; 3 2 1 20]; x_ex = ones(size(A,1),1); b = A *x_ex;if isequal(A,A)disp(Aesimmetrica) else disp(Anonesimmetrica) end if all(eig(A)>0)disp(Aedefinitapositiva) else disp(Anonedefinitapositiva:) disp(eig(A)=) eig(A) end dom_diag = [true true]; n = size(A,1); for i = 1:nif abs(A(i,i)) < sum(abs(A(i,[1:i-1, i+1:n])))dom_diag(1) = false; end if abs(A(i,i)) < sum(abs(A([1:i-1, i+1:n],i)))dom_diag(2) = false; end end if dom_diag(1)disp(Aeadominanzadiagonaleperrighe) else disp(Anoneadominanzadiagonaleperrighe) end if dom_diag(2)disp(Aeadominanzadiagonalepercolonne) else disp(Anoneadominanzadiagonalepercolonne) end D = diag(diag(A)); E = -tril(A,-1); F = -triu(A,1); BJ = D(E+F); rhoJ = max(abs(eig(BJ))) BGS = (D-E)\F; rhoGS = max(abs(eig(BGS))) Notiamo che Aè denita positiva ma non simmetrica. Tuttavia, è a dominanza diagonale stretta per righe (e anche per colonne), dunque possiamo aermare che entrambi i metodi convergono. Ciò è confermato dal valore dei raggi spettrali delle matrici di iterazione :(B J) = 0 :2950; (B GS) = 0 :0821. b)Abbiamonorm(BJ)=jjB Jjj 2= 0 :4034 0)disp(Matricesimmetricadefinitapositiva) else disp(Matricesimmetricamanondefinitapositiva) end else disp(Matricenonsimmetrica) end disp(Numerodicondizionamentodellamatrice:) max(v)/min(v) Matlab restituisce a schermo : Matrice A: Matrice simmetrica definita positiva Numero di condizionamento della matrice: ans = 336.2412 c)Si riporta il listato della funzione richiesta :7 function [x, k] = richardson(A, b, P, x0, tol, nmax, alpha)% %MetododiRichardsonstazionarioprecondizionato %odinamicoprecondizionato(gradienteprecondizionato) % %Parametridiingresso: % %A:matricedisistema %b:vettoreterminenoto %P:precondizionatore %x:guessiniziale %tol:tolleranzacriteriodarresto %nmax:numeromassimodiiterazioniammesse %alpha:parametrodiaccelerazione;senonassegnatosiconsidera %ilmetododinamico(gradienteprecondizionato) % %Parametridiuscita: % %x:soluzione %k:numerodiiterazioni % n = length(b); if ((size(A,1) ~= n) || (size(A,2) ~= n) || (length(x0) ~= n))error(Dimensioniincompatibili) end x = x0; k = 0; r = b - A *x;res_norm = norm(r) / norm(b); while ((res_norm > tol) && (k < nmax))z = P \ r; x = x + alpha *z;r = b - A *x; %equivalentea:r=r-alpha *A *z;res_norm = norm(r) / norm(b); k = k + 1; end if (res_norm < tol)fprintf(Richardsonconvergein%diterazioni\n, k); else fprintf(Richardsonnonconvergein%diterazioni.\n, k) end d)Si può usare la funzione richardson.mper risolvere il sistema lineare del caso in esame come segue :nmax = 1e4; P = eye(n); %Precondizionatoredisp(Richardson:P=I,alpha=0.20) alpha = 0.2; B_alpha = eye(n) - alpha *A; %inv(P)=Idisp(Raggiospettrale:) max(abs(eig(B_alpha))) [x02, k02] = richardson(A, b, P, x0, tol, nmax, alpha); disp(Richardson:P=I,alpha=0.33) alpha = 0.33; B_alpha = eye(n) - alpha *A; %inv(P)=Idisp(Raggiospettrale:) max(abs(eig(B_alpha))) [x033, k033] = richardson(A, b, P, x0, tol, nmax, alpha); 8 disp(Richardson:P=I,alpha=OPT) alpha = 2/(max(v)+min(v)); B_alpha = eye(n) - alpha *A; %inv(P)=Idisp(Raggiospettrale:) max(abs(eig(B_alpha))) [xopt, kopt] = richardson(A, b, P, x0, tol, nmax, alpha); Dai risultati si può dedurre che : seP = Iealpha = 0.20il raggio spettrale della matrice di iterazione è minore di uno (0.9963) ed il metodo converge dopo 3074 iterazioni ; seP = Iealpha = 0.33il raggio spettrale della matrice di iterazione è maggiore di uno (1.0581) ed il metodo non converge ; seP = Iealpha = 2/(max(v)+min(v))(valore ottimo) il raggio spettrale della matrice di iterazione è 0.9941, leggermente inferiore a quello del primo caso ed il numero di iterazioni richieste scende a 1921. e)Si può dimostrare la seguente stima : kx( k) x exk A Ck kx(0) x exk A; doveC=K 2( A)1K 2( A) + 1: (5) Possiamo vericarne la validità dopokpassi con i seguenti comandi :x_ex = A \ b; [xopt, kopt] = richardson(A, b, eye(n), x0, tol, nmax, alpha); err0_normA = sqrt((x0-x_ex) *(A *(x0-x_ex)));C = (cond(A)-1) / (cond(A)+1); errK_normA = sqrt((xopt-x_ex) *(A *(xopt-x_ex)))errK_normA