Mathematical Engineering - Matematica Numerica

Exercise - 10 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 10 - Soluzione Metodi iterativi per sistemi lineari (3) Esercizio 1 Si consideri il sistemaAx=bdove la matriceAed il temine notobsi ottengono con i seguenti comandiA = delsq( numgrid(B,50) ); b = ones( size(A,1),1 ); La matrice A2R1814 1814 è simmetrica e denita positiva con struttura sparsa, come si può visualizzare con il comandospy(A). Essa deriva dalla discretizzazione mediante dierenze nite del problema di Laplace sulla griglia di nodi visualizzabile conspy(numgrid (B, 50)). Per ciascuno dei seguenti punti, usare il comando Matlabpcgche implementa il metodo del gradiente coniugato precondizionato, imponendo una tolleranza pari a10 12 sul residuo, indi memorizzare il vettore dei residui ottenuti in funzione del numero di iterazioni (si rimanda all'help del comando). Visualizzare quindi le tre curve di abbattimento del residuo in scala semilogaritmica consemilogy. a)Calcolare il numero di condizionamento diA(usare il comandocondest) e risolvere il sistemaAx=b tramitepcg(con precondizionatore identico). b)Vericare che la parte tridiagonalePdiAè simmetrica e denita positiva. Calcolare il condizionamento del sistema precondizionato con la matricePe risolvere il sistema tramitepcg. c)Costruire algebricamente il sistema con precondizionatore centratoA x =b , doveA =^ H T A^ H 1 , b =^ H T b,x =^ Hx. La matrice^ Hè la matrice triangolare superiore ottenuta tramite la fattorizza- zione di Cholesky incompleta (senza ll-in) della matriceA, calcolata tramite il seguente comando,H=ichol(A); spy(H); Calcolare il numero di condizionamento associato al sistema precondizionato e risolvere tale sistema conpcg. Cosa emerge dal confronto ? Soluzione 1a)Costruiamo la matriceAe il vettorebcome segue :%CreiamoAeb G = numgrid(B,50); A = delsq( G ); b = ones( size(A,1),1 ); %Nota:sesivuolevisualizzarelagrigliausataper %costruireA... figure spy(numgrid(B,50)); %Nota:sesivuolevisualizzarelasoluzionedellequazionedi %Laplacesullagrigliadata...(vedereanchelademodelsqdemo) 1 x = A\b; U = G; U(G>0) = full(x(G(G>0))); figure(11); colormap(cool); mesh(U) axis square ij colorbar Usiamo innanzitutto il metodo del gradiente coniugato non precondizionato (ovvero precondizionatore uguale all'identità), dopo aver calcolato il condizionamento diA. A riguardo, conviene utilizzare il comandocondestper matrici sparse, molto meno costoso dal punto di vista computazionale rispetto al comandocond. Segnaliamo tuttavia checondestrestituisce il condizionamento in norma 1, ovvero K1(A) .%CondizionamentodiA: condA = condest(A) condA = 439.8246 %Sistemanonprecondizionato: % vettoredeiresidui:resvec% maxiterazioni:2000% tolleranzasulresiduo(normalizzato):1e12%PCG [x, flag, relres, iter, resvec1] = pcg(A, b, 1e 12, 2000);iter1 = iter iter1 = 136 Il vettore dei residui memorizzato in resvec1potrà poi essere visualizzato consemilogy(si veda la Figura 2). b)Per quanto riguarda l'uso del precondizionatorePdato dalla parte tridiagonale diA(che è invece pentadiagonale), si ottiene :%Precondizionatore=partetridiagonale: P = diag(diag(A)) + diag(diag(A, 1),1) + diag(diag(A,1),1);%PesimmetricapercheloeA;inoltreedefinita %positiva.Infattiilminoredegliautovalori(chesicalcola %comesegue)epositivo: eigs(P, 1,SA) ans = 2.0041 %Condizionamentomatriceprecondizionata condPreA = condest(P\A) condPreA = 221.0121 %PCG,precondizionatoreP [x, flag, relres, iter, resvec2] = pcg(A, b, 1e 12, 2000, P);iter2 = iter iter2 = 95 2 Osserviamo che il condizionamento del sistema precondizionato si riduce rispetto al caso precedente ; coerentemente, anche il numero di iterazioni necessarie per convergere alla soluzione è minore. c)Inne, utilizziamo la fattorizzazione di Cholesky incompleta per precondizionare il metodo del gra-diente coniugato. La fattorizzazione di Cholesky incompleta, produce un fattore di Cholesky inesatto^ Hche però mantiene le stesse caratteristiche di sparsità della matriceA(Figura 1 ; per il problema dell l-inassociato ai metodi diretti si vedano le soluzioni alle precedenti Esercitazioni).%CalcolodelfattorediCholeskyIncompleto H = ichol(A); spy(H); %Calcolodellamatriceprecondizionatacentralmente Astar = (H)\A/H; %edelvettorebstar bstar = H *b;%Condizionamentomatriceprecondizionata condAstar = condest(Astar) condAstar = 64.2629 %PCGsulsistemaprecondizionatocentralmente: [xstar, flag, relres, iter, resvec3] = pcg(Astar, bstar, 1e 12, 2000);%Calcolodixapartiredaxstar x = H\xstar; iter3 = iter iter3 = 50 %Equivalentemente: %(implementazioneefficentedelmetodoprecondizionato %centralmente) [x, flag, relres, iter, resvec3] = pcg(A, bstar, 1e 12, 2000, H, H);Abbiamo qui sfruttato il fatto che, quando si danno in input a pcgdue precondizionatori diversi,Figure 1  La matriceA(a sinistra) e il suo fattore di Cholesky incompleto^ H(a destra)30 500 1000 1500 0 200 400 600 800 1000 1200 1400 16001800 nz = 8710 0 500 1000 1500 0 200 400 600 800 1000 1200 1400 16001800 nz = 5262 Figure 2  Abbattimento del residuo normalizzato rispetto al numero di iterazioni, metodo del gradiente coniugato precondizionato. Le curve corrispondono a : metodo non precondizionato (cerchi), precondiziona- torePdato dalla parte tridiagonale diA(croci) e precondizionatore simmetrico ottenuto a partire da^ H (quadrati) il comando considera il loro prodotto come precondizionatore. Osserviamo che il numero di condi- zionamento in norma 1 diA , ovveroK 1(A ), è notevolmente più piccolo di quello diA. Dunque, il precondizionatore centrato basato su^ Hè estremamente ecace. Il confronto fra i tre casi visti (Figura 2) mostra come un buon precondizionatore può migliorare drasticamente la velocità di conver- genza del metodo iterativo utilizzato. Tuttavia, precondizionatori di complessità crescente implicano un aumento del tempo di calcolo associato ad ogni singola iterazione (e, nel caso del precondizionatore di Cholesky incompleto, al calcolo della fattorizzazione incompleta). Esercizio 2 Vogliamo confrontare le prestazioni del metodo del gradiente coniugato precondizionato applicato a un sis- tema lineare derivante dalla discretizzazione a elementi niti di un problema, al variare della matrice di precondizionamento. La matrice e il vettore del sistema lineare vengono generati con i comandiA = gallery(poisson, 40); b = ones(size(A,1),1); Una buona scelta per il precondizionatore Pè quella di considerare una matrice simile adAma facilmente invertibile. La sceltaP=RT R, doveRè la matrice triangolare superiore ottenuta mediante la fattorizzazione di Cholesky diA, è ottimale ma decisamente inecace a causa del fenomeno del ll-in. Una possibile alterna- tiva consiste nel considerare la fattorizzazione incompleta di Cholesky. L'idea di base è quella di determinare una matrice~ RRe porreP=~ RT ~ R. Due possibili alternative sono le seguenti : a)Cholesky no-l l: si pone~ Rij= 0 seA ij= 0 ; b)Cholesky droptol: si pone~ Rij= 0 sejR ijj < droptol, dovedroptolindica una tolleranza da passare al metodo. Nel primo caso, si usano i seguenti comandiMatlabper generare il precondizionatore :40 20 40 60 80 100 120 140 10 −14 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 10 2 opts.type =nofill; P = ichol(A, opts); %precondizionatorementre nel secondo si usano i seguenti comandi : opts.type =ict; opts.droptol = 1e 2;P = ichol(A, opts); %precondizionatorePer risolvere il sistema lineare, si usa in entrambi i casi la seguente sintassi del comando pcg:[X,FLAG,RELRES,ITER,RESVEC] = pcg(A,b,tol,maxIter,P, P); Confrontare, in termini di ecienza (tempo di calcolo, numero di iterazioni), accuratezza (andamento dei residui) e pattern della matrice di precondizionamento le seguenti scelte : 1.metodo non precondizionato ; 2.precondizionatore dato dalla diagonale diA; 3.precondizionatore Cholesky no-ll ; 4.precondizionatore Cholesky droptol =10 2 ; 5.precondizionatore Cholesky droptol =10 4 . Soluzione 2 Il codice richiesto è il seguente :A = gallery(poisson, 40); b = ones(size(A,1),1); figure spy(A) maxIter = 800; tol = 1e 8;%%NoPrecon timePrecon(1) = 0; timeR = tic; [X,FLAG,RELRES,ITER(1),RESVEC1] = pcg(A,b,tol,maxIter, speye(size(A,1),size(A,2))); timeSolve(1) = toc(timeR) %%DiagPrecon timeP = tic; M = diag(diag(A)); timePrecon(2) = toc(timeP) timeR = tic; [X,FLAG,RELRES,ITER(2),RESVEC2] = pcg(A,b,tol,maxIter,M); timeSolve(2) = toc(timeR) %%ICHOLPrecon timeP = tic; opts.type =nofill; opts.droptol = 1e 2;M1 = ichol(A, opts); timePrecon(3) = toc(timeP) figure spy(M1) 5 timeR = tic; [X,FLAG,RELRES,ITER(3),RESVEC3] = pcg(A,b,tol,maxIter,M1, M1); timeSolve(3) = toc(timeR) %%ICHOLPrecon timeP = tic; opts.type =ict; opts.droptol = 1e 2;M2 = ichol(A, opts); timePrecon(4) = toc(timeP) figure spy(M2) timeR = tic; [X,FLAG,RELRES,ITER(4),RESVEC4] = pcg(A,b,tol,maxIter,M2, M2); timeSolve(4) = toc(timeR) %%ICHOLPrecon timeP = tic; opts.type =ict; opts.droptol = 1e 4;M3 = ichol(A, opts); timePrecon(5) = toc(timeP) figure spy(M3) timeR = tic; [X,FLAG,RELRES,ITER(5),RESVEC5] = pcg(A,b,tol,maxIter,M3, M3); timeSolve(5) = toc(timeR) [M4,p, s] = chol(A); Risulta poi possibile costruire una tabella che riporti le informazioni richieste : T = table(timePrecon,timeSolve,timePrecon+timeSolve,ITER,RowNames, ...{NoPrecon;Diagonal;ICno fill;ICdroptol1e2;ICdroptol1e4},...VariableNames,{timePrecon;timeSolve;TotalTime;IterationNumber}) 6 Figure 3  Abbattimento del residuo normalizzato rispetto al numero di iterazioni, metodo del gradiente coniugato precondizionato. Le curve corrispondono a : metodo non precondizionato, precondizionatore dato dalla parte diagonale diA, precondizionatori ottenuti con fattorizzazioni di Cholesky incompletetimePrecon timeSolve TotalTime IterationNumber __________ _________ _________ _______________ No Precon 0 0.053385 0.053385 74 Diagonal 0.00024075 0.015674 0.015915 74 IC no fill 0.00099568 0.013949 0.014945 35IC droptol 1e 2 0.00029248 0.0044056 0.0046981 20IC droptol 1e 4 0.0016928 0.010005 0.011697 6Inne, il seguente graco riporta l'andamento dei residui nei vari casi. figure semilogy(RESVEC1(2: end)./RESVEC1(2),ob,LineWidth,2);hold on semilogy(RESVEC2(2: end)./RESVEC2(2),k,LineWidth,2);semilogy(RESVEC3(2: end)./RESVEC3(2),r,LineWidth,2);semilogy(RESVEC4(2: end)./RESVEC4(2),g,LineWidth,2);semilogy(RESVEC5(2: end)./RESVEC5(2),m,LineWidth,2);grid on legend(NoPrecon,Diagonal,IncompleteCholno fill,...IncompleteCholdroptol1e 2,IncompleteCholdroptol1e4)70 50 100 150 200 250 300 350 10−10 10−8 10−6 10−4 10−2 100    No Precon Diagonal Incomplete Chol no−fill Incomplete Chol droptol 1e−2 Incomplete Chol droptol 1e−4 Esercizio 3 (P)A partire dai modelli di main le e dataleEser10_DiffReactMixedBC_template.me Eser10_DiffReactMixedBC_data_template.m, implementare inPoliFEMla soluzione del seguente problema di diusione e reazione, con condizioni miste di Dirichlet-Neumann ai due estremi dell'intervallo(0;1), usando elementi niti quadratici :( (u0 )0 + u=f(x);0< x