Mathematical Engineering - Matematica Numerica

Exercise - 06 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 6 - Soluzione Metodi diretti per sistemi lineari (2) Esercizio 1 Si consideri la matrice A = " b 1c conb; c2R, e0< "1, ovvero molto piccolo. a)Si applichi il metodo di eliminazione di Gauss (senza pivoting) per determinare la fattorizzazione LUdiA. b)Si applichi il metodo di eliminazione di Gauss con pivotazione per righe per determinare la fattoriz-zazione LU diA. c)Consideriamo la matrice non singolare A =2 41 1 + 0 :5
10 15 3 2 2 20 3 6 43 5; si esegua con Matlab la fattorizzazione LU mediante il metodo di eliminazione di Gauss senza pivoting per righe (usando la funzionelugauss.m) e con pivoting per righe (usando il comando Matlablu). Quale dierenza si nota ? Sapreste motivarla ? Soluzione 1a)Applicando il MEG (senza pivoting) alla matriceAsi hanno L = 1 0 1="1 ;U = " b 0cb=" con elemento pivotalea(1) 11= ". Il calcolo della fattorizzazione LU con MEG al calcolatore genera quindi errori di arrotondamento rilevanti per"molto piccolo dovuti al calcolo diu 22= cb=", essendo1=" molto grande. In pratica, se il moltiplicatore è grande, può esserci perdita di cifre signicative nella sottrazione. b)Applicando il MEG con pivotazione per righe diA, ovvero usando P = 0 1 1 0 si ottengono : la matrice permutata~ A = PA = 1c " b ; ~ L = 1 0 "1 ;~ U = 1c 0b"c con elemento pivotale~ a(1) 11= 1 . In tal caso, il calcolo della fattorizzazione LU tramite MEG con pivoting per righe al calcolatore consente di contenere gli errori di arrotondamento ; infatti, per" molto piccolo, si ha~ u 22= b"c'b.1 c)Con i seguenti comandi, si ottiene : A = [1 1+0.5e-15 3; 2 2 20; 3 6 4] [L,U] = lugauss(A) L =1.0e+15 * 0.0000 0 0 0.0000 0.0000 0 0.0000 -3.3777 0.0000 U =1.0e+16 * 0.0000 0.0000 0.00000 -0.0000 0.0000 0 0 4.7288 Durante il calcolo della fattorizzazione con la funzionelugauss.mnon si generano elementi pivotali nulli. Tuttavia, i fattoriLeUcalcolati sono assai inaccurati, come si verica calcolandoALU: A-L*U ans =0 0 0 0 0 0 0 0 4 Se indicassimo con"= 0:5
10 15 , e procedessimo con carta e penna, otterremmo ~ L =2 41 0 0 2 1 0 332 (1 1=") 13 5;~ U =2 41 1 + "3 02"14 0 0 2126="3 5 e dunque lo stesso eetto riscontrato nel punto a). L'elemento pivotaleu 22= a(2) 22= 2"è molto piccolo, e genera errori di arrotondamento grandi. Eettuando invece il pivoting, implementato nella funzioneludi Matlab, si ottiene [L1,U1,P1] = lu(A) L1 =1.0000 0 0 0.6667 1.0000 0 0.3333 0.5000 1.0000 U1 =3.0000 6.0000 4.00000 -2.0000 17.3333 0 0 -7.0000 P1 =0 0 1 0 1 0 1 0 02 e in questo caso i fattori calcolati risultano accurati. Un'interpretazione quantitativa del ruolo del pivoting per righe è illustrato di seguito. Accuratezza della fattorizzazione LU con pivoting per righe Per comprendere l'importanza del pivoting parziale per righe, necessario anche in presenza di elementi pivotali non nulli, enunciamo il seguente risultato.Teorema [Wilkinson, 1961] Data la matriceA2Rn n e il sistema lineareAx=b, sia^ xla soluzione ottenuta al calcolatore tramite il metodo della fattorizzazione LU, ovvero applicando il MEG. Allora, si hanno : (A+A)^ x=bekAk 1 C n3 n" Mk Ak 1; doven=max i;j;k=1;:::;nj a( k) ijjmax i;j=1;:::;nj a ijj è il fattore di crescita,C >0una costante positiva e" Ml'epsilon macchina ; invece, n=max i;j;k=1;:::;nj ~ a( k) ijjmax i;j=1;:::;nj a ijj se il MEG è applicato con una tecnica di pivoting per righe o totale, dove~ Aè la matrice ottenuta permutando le righe (ed eventualmente le colonne) diA.In generale,  n 1, ma può essere arbitrariamentegrandese il MEG viene applicato senza alcuna tecnica di pivoting. Al contrario, n 2n 1 se si utilizza il MEG con pivoting per righe, mentre n n1 =2 n(log n)=4 per il MEG con pivoting totale. Tuttavia, tali stime sono pessimistiche quando una tecnica di pivoting viene applicata ; in tali casi infatti, si ha generalmente n 50. Si deduce che lo scopo delle tecniche di pivoting consiste nel contenere e limitare il fattore di crescita ne, in ultima istanza, mantenere stabile il metodo numerico. Quando tecniche di pivoting sono applicate contestualmente al MEG per determinare la fattorizzazione LU di A, l'assunzioneA0(fatta per dimostrare alcune stime nell'analisi a priori e a posteriori) è giusticabile. Tornando al caso dei punti a) e b), e prendendo per comoditàb=c= 1, applicando il MEG senza pivoting, si ottiene n= j11="j 1=", da cui, usando il risultato enunciato nel Teorema, la perturbazione (si ha kAk 1= 2 ,n= 2) kAk 1 16C" M" può risultare non trascurabile, ovverogrande, quando"èpiccoloo comparabile al valore dell'epsilon macchina "M. Al contrario, quando viene applicato il MEG con pivoting per riga, si ottiene che  n= 1 , da cui si deduce che la perturbazionekAk 1 16C " M è del tutto trascurabile.3 Esercizio 2 Utilizzare opportunamente i comandi Matlabonesediagper costruire una matriceA2Rn n conn= 20 tale che,a 1;i= 1 ea i;1= 1 peri= 1; : : : ; n,a ii= 4 peri= 2; : : : ; ned innea i;i+1= 1ea i+1;i= 1per i= 2; : : : ; n1. Tutti gli altri coecienti diAsono nulli. Tale matrice si dicesparsapoiché il numero degli elementi non-nulli diAèO(n)n2 , dunque molto elevato1 . Per visualizzare la posizione degli elementi non nulli diAsi utilizzi il comandospy. Calcolare inoltre la fattorizzazioneLUdiA. Utilizzare il comandospyper determinare se le matriciLed Usiano o meno sparse. A vostro parere, le caratteristiche di sparsità dei fattoriLeUsono vantaggiose o svantaggiose ? Soluzione 2 Se la matrice del sistema lineare è sparsa (e non strutturata),A2Rn n ha un numero di elementi non nulli dell'ordine din. In tal caso durante il processo di fattorizzazione si possono generare un gran numero di elementi non nulli in corrispondenza di elementi nulli della matrice di partenza. Questo fenomeno, noto con il nome di riempimento ol l-in, è particolarmente gravoso da un punto di vista computazionale perché non consente di memorizzare la fattorizzazione di una matrice sparsa nella stessa area di memoria necessaria per memorizzare la matrice stessa. La tecnica di pivotazione totale può essere convenientemente applicata per prevenire e/o contenere il fenomeno del ll-in. La matriceAsi denisce attraverso i seguenti comandi n = 20; A = 4*diag(ones(n,1)) - diag(ones(n-1,1),-1)- diag(ones(n-1,1),1); A(1,:) = ones(1,n); A(:,1) = ones(n,1); e la posizione dei suoi elementi non nulli si visualizza ad esempio come segue (Figura 1) figure; spy(A); title('A'); Calcoliamo ora la fattorizzazioneLUdiAe visualizziamo le matriciLeU [L,U,P] = lu(A); figure; spy(L); title('L'); figure; spy(U); title('U'); Otteniamo così la Figura 1. Si osserva che le matriciLeUnon sono sparse, poiché il loro numero di elementi non nulli è circan2 =2. Il fenomeno per cui una matrice sparsa dà luogo a fattori (LU o Cholesky) generalmente non sparsi è chiamatol l-in(riempimento delle matrici con elementi non nulli). Dunque, in generale l'algoritmo di fattorizzazione sostituisce ad elementi nulli inAdegli elementi non nulli inLeU.Figure 1  Rappresentazione degli elementi non nulli diA,LeU1. In tal caso è vantaggioso memorizzare solo gli elementi non-nulli di A. In Matlab ciò viene fatto attraverso il comando sparse, ad esempioB=sparse(A).40 5 10 15 20 nz = 76 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 nz = 210 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 nz = 210 0 2 4 6 8 10 12 14 16 18 20 Il l l-inè svantaggioso in termini di gestione della memoria. Nelle applicazioni al calcolo scientico, ovvero conngrande, si ha quasi sempre a che fare con matrici sparse per le quali viene riservata memoria solo per gli elementi non nulli, che tipicamente sonoO(n) O(n2 ). Ad esempio, denendo le matrici del presente laboratorio pern= 100, si osserva che l'occupazione di memoria nel formato sparso di Matlab consente un notevole risparmio per la matriceA, ma non per i fattori LU : n = 100; A = 4*diag(ones(n,1)) - diag(ones(n-1,1),-1)- diag(ones(n-1,1),1); A(1,:) = ones(1,n); A(:,1) = ones(n,1); [L,U,P] = lu(A); As = sparse(A); Ls = sparse(L); Us = sparse(U); % Esaminiamo il Workspace e i Bytes di memoria allocati whos Name Size Bytes Class Attributes A 100x100 80000 double As 100x100 6332 double sparse L 100x100 80000 double Ls 100x100 61004 double sparse P 100x100 80000 double U 100x100 80000 double Us 100x100 61004 double sparse n 1x1 8 double Si noti come la memoria allocata perAsè meno di un decimo di quella allocata perA; lo stesso guadagno non si ottiene invece perLseUs, in quanto la loro struttura non è sparsa. Quando la fattorizzazione LU comporta una generazione di elementi non nulli dell'ordine diO(n2 )(l l-in), l'allocazione ulteriore di memoria è proibitiva perngrande e rende impossibile l'uso della fattorizzazione stessa. Per questo, possibili rimedi consistono in un riordino opportuno degli elementi diAin maniera tale da minimizzare ill l-indei fattori L e U. La tecnica di pivotazione totale può essere convenientemente applicata per prevenire e/o contenere il fenomeno del ll-in. La pivotazione totale applicata alla matriceA2Rn n consiste nell'introdurre, oltre alla matrice di permutazione per righeP2Rn n da pre-moltiplicare adA, una matrice di permutazione per colonne Q2Rn n da post-moltiplicare adAcomeAQ. La matrice di permutazione per colonneQè ortogonale, cioèQT =Q 1 (per cuiQT Q=I) ; seQ=I, allora non vi sono permutazioni di colonne della matriceA. Applicando la tecnica di pivotazione totale per la determinazione della fattorizzazione LU della matrice non singolareA, si ha P AQ= LU: Osserviamo che la permutazione totale comporta costi computazionali aggiuntivi rispetto alla pivotazione per righe ; infatti, per ognik= 1; : : : ; n1, la ricerca dell'elemento pivotalemiglioreavviene su in insieme di(n+ 1k)2 elementi  della sottomatriceA( k) 2R( n+1k)(n+1k)  che è molto più grande di quello corrispondente alla rigakdella sottomatriceA( k) , ovvero gli(n+ 1k)elementif(A( k) )ikgn i=k. Per tale motivo, in Matlab la pivotazione totale viene operata solo su matrici sparse. Nel caso in esame, con i comandi As = sparse(A); Ls = sparse(L); Us = sparse(U); [Ls,Us,Ps,Qs] = lu(As) è possibile ottenere fattoriLeUsparsi. Nella gura 2 sono riportati i pattern diLseUs(ovvero dei fattori LeUin formato sparso), oltre che della matriceP AQ. Risulta evidente l'eetto del pivoting totale sulla sparsità dei fattoriLeUottenuti in presenza di pivoting totale.5 Figure 2  Rappresentazione degli elementi non nulli diL,Ue diP AQ. Concludiamo osservando un'importante proprietà della fattorizzazione LU. Una matriceA2Rn n è dettaa banda, conampiezza di bandap, sea ij= 0 perjijj p. La fattorizzazione LU diAconserva la struttura a banda, ossia, postoL = (l ij) eU = (u ij) ) aij= 0 perjijj p=)( lij= 0 perijp; uij= 0 perjip: Per questo motivo, l'applicazione del metodo di eliminazione di Gauss alla matrice A (che ha banda di ampiezzan) comporta il riempimento della banda di molti elementi diversi da zero (ll-in della banda) ; è dunque coerente trovare dei fattoriLeUavendo entrambi banda di ampiezzan. Esercizio 3 Per limitare il fenomeno del riempimento (l l-in) è possibile eettuare operazioni preliminari di permutazione tra righe, che prendono il nome direordering2 . L'algoritmo di Cuthill-McKee (e la sua variante inversa) permette di permutare una matrice sparsa simmetrica in una matrice a banda più compatta. La variante inversa di tale algoritmo è implementata in Matlab nel comandosymrcm; in particolare, r = symrcm(S) restituisce una permutazionertale cheS(r,r)tende ad avere tutti gli elementi non nulli vicini alla diagonale. a)Generare con il comandoB = buckyuna matriceB2R60 60 dallagallerydi Matlab. Vericare che tale matrice ha una banda piuttosto larga. b)Eettuare la fattorizzazione LU della matriceB, e mostrare il pattern dei fattori L ed U. c)Eettuare il reordering della matriceBcon il comandosymrcme, indicando conBrla matrice riordinata, determinare la fattorizzazione LU di tale matrice. Mostrare il pattern dei fattori L ed U, confrontando i risultati con quelli ottenuti al punto precedente. Soluzione 3a)Con i seguenti comandi si possono denire la matriceBe quella ottenuta con l'algoritmo di reordering.2. Da non confondersi con il pivoting60 2 4 6 8 10 12 14 16 18 20 nz = 56 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 nz = 40 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 nz = 76 0 2 4 6 8 10 12 14 16 18 20 B = bucky; figure(); subplot(1,2,1), spy(B), title('B') p = symrcm(B); Br = B(p,p); subplot(1,2,2), spy(Br), title('B(p,p)') La matriceBrappresenta la matrice di adiacenza degli elementi di una cupola geodetica, detta di Buckminster Fuller dal nome del suo inventore. La struttura del grafo si può visualizzare con i comandi G = graph(bucky); p = plot(G) ed è riportata nella Figura 3.Figure 3  Grafo di adiacenza degli elementi della struttura codicata nella matricebuckyFigure 4  Pattern della matriceB(a sinistra) e della matrice riordinataB r(a destra)7 0 10 20 30 40 50 60 nz = 180 0 10 20 30 40 50 60 B 0 10 20 30 40 50 60 nz = 180 0 10 20 30 40 50 60 B(p,p) L'ampiezza della banda può essere calcolata con i comandi [i,j] = find(B); bw = max(abs(i-j)) +1; b)Eseguiamo ora la fattorizzazione LU della matriceBe della corrispondente matrice riordinataB r: [L,U,P] = lu(B) [Lr,Ur,Pr] = lu(Br) figure() subplot(1,2,1), spy(L), title('L') subplot(1,2,2), spy(U), title('U') figure() subplot(1,2,1), spy(Lr), title('Lr') subplot(1,2,2), spy(Ur), title('Ur') Si ottengono i graci riportati in Figura 5.Figure 5  Pattern dei fattoriLedU(P B=LU)Figure 6  Pattern dei fattoriL red U r( P rB r= L rU r) c)Si può dunque vericare come l'ampiezza della banda sia pari a 36 perB, mentre si riduce a 12 per la matrice riordinataB r. I fattori ottenuti con il comando LU hanno presentano la stessa ampiezza di banda delle matrici di cui sono la fattorizzazione, ossiaP BoppureP B r.80 10 20 30 40 50 60 nz = 451 0 10 20 30 40 50 60 L 0 10 20 30 40 50 60 nz = 419 0 10 20 30 40 50 60 U 0 10 20 30 40 50 60 nz = 448 0 10 20 30 40 50 60 Lr 0 10 20 30 40 50 60 nz = 447 0 10 20 30 40 50 60 Ur Esercizio 4 Si consideri una matriceAtridiagonale (caso particolare di matrice a banda) e invertibile : A =2 6 6 6 6 4a 1c 10 e2a 2: :: :::c n1 0e na n3 7 7 7 7 5; dovefa ign i=1, fc ign 1 i=1, fe ign i=22 R. Si vuole fattorizzare tale matrice nel prodottoA = LU(caso particolare difattorizzazione LU) doveLeUsono due matrici bidiagonali, della forma L =2 6 6 6 41 0 21 :::: :: 0 n13 7 7 7 5; U =2 6 6 6 6 4 1 10 2: :: ::: n1 0 n3 7 7 7 7 5; allo scopo di risolvere il sistema lineareAx=bin modo eciente. Essendo infattiLeUcasi particolari di matrici rispettivamente triangolari inferiori e triangolari superiori, tale riscrittura permette di semplicare la risoluzione del sistema, adottando gli algoritmi disostituzione in avantiesostituzione al l'indietro. Questi algoritmi prevedono la risoluzione sequenziale di ciascuna riga di un sistema triangolarenn, andando rispettivamente, dalla prima all'ultima, o dall'ultima alla prima. a)Fornire l'espressione dei coecienti i,  ie iin funzione dei coecienti di A. b)Applicare le formule trovate in a) e costruire un algoritmo per risolvere il sistema lineareAx=b, con b= [b 1; : : : ; b n]T 2Rn , mediante sostituzioni in avanti e all'indietro (algoritmo di Thomas). c)Quante operazionioating-pointsono richieste dall'algoritmo di Thomas ? d)Implementare in Matlab, mediante unafunction, l'algoritmo di Thomas per risolvere un generico sistema lineare tridiagonale. L'intestazione dellafunctiondovrà essere : function [L,U,x] = thomas(A,b) Soluzione 4a)I coecienti i,  ie isi trovano imponendo che risulti A = LU. Ne derivano le relazioni seguenti : 1= a 1;  i=e i i1;  i= a i  ic i1; i= c i: (1) b)La soluzione del sistemaAx=bsi riconduce a risolvere due sistemi bidiagonali,Ly=beUx=y, per i quali abbiamo le formule seguenti : (Ly=b) :y 1= b 1; y i= b i  iy i1; i = 2; : : : ; n; (Ux=y) :x n=y n n; x i=y i c ix i+1 i; i =n1; : : : ;1:(2) c)L'algoritmo di Thomas richiede8n7ops:3(n1)per la fattorizzazione (1) e5n4per le sostituzioni (2). d)La funzione che implementa il metodo di Thomas è riportata nel seguente algoritmo function [L,U,x] = thomas(A,b) % function [L,U,x] = thomas(A,b) % % Risoluzione di un sistema tridiagonale Ax = b con % l'algoritmo di Thomas9 % % INPUT: % A: matrice tridiagonale % b: termine noto % OUTPUT: % L: matrice triangolare inferiore bidiagonale % U: matrice triangolare superiore bidiagonale % x: soluzione N=length(b); if ((size(A,1)~=N)||(size(A,2)~=N))error('Dimensioni di A e f incompatibili') end if (sum(sum(abs(tril(A,-2)+triu(A,2))))~=0)error('ERRORE: matrice A non tridiagonale') end a = diag(A); e = [0; diag(A,-1)]; c = diag(A,1); alpha = zeros(N,1); beta = zeros(N,1); gamma = c; alpha(1) = a(1); for i=2:Nbeta(i) = e(i)/alpha(i-1); alpha(i) = a(i)-beta(i)*c(i-1); end L = diag(ones(N,1), 0) + diag(beta(2:end), -1); U = diag(alpha,0) + diag(gamma,1); % soluzione di Ly = b con sostituzioni in avanti y=zeros(N,1); y(1)=b(1); for i=2:Ny(i) = b(i) - beta(i)*y(i-1); end % soluzione di Ux = y con sostituzioni all'indietro x=zeros(N,1); x(N)=y(N)/alpha(N); for i=N-1:-1:1 x(i)=(y(i)-c(i)*x(i+1))/alpha(i); end10 Esercizio 5 Si consideri la seguente matrice A =2 41 0 2 1 0 1 33 5: a)Fornire un intervallo di valori ditale che la matriceAammetta fattorizzazione di Cholesky. b)Calcolare analiticamente la fattorizzazione di Cholesky di A in funzione dinell'insieme dei valori ammissibili trovato al punto precendente. Vericare il calcolo nel caso= 1utilizzando il comando choldi Matlab. c)Sempre nel caso= 1, risolvere in Matlab i tre sistemi lineari seguenti, utilizzando opportunamente la fattorizzazione di Cholesky trovata al punto precedente e le funzionifwsubebksub: Ax i= e i; i = 1;2;3; dovee 1= [1 ;0;0]T ,e 2= [0 ;1;0]T ee 3= [0 ;0;1]T . Utilizzare questi risultati per calcolareA 1 . Si spieghi perché la fattorizzazione di Cholesky è conveniente dal punto di vista computazionale per il calcolo dell'inversa diA. Soluzione 5a)Per potere applicare la fattorizzazione di Cholesky, è necessario cheAsia simmetrica e denita positiva. Essendo la simmetria automaticamente vericata, per vericare la denita positività diAapplichiamo il criterio di Sylvester (in maniera equivalente si può imporre che gli autovalori siano tutti positivi). Precisamente, imponiamo che tutti i minori principali diAsiano positivi. Ciò fornisce il seguente sistema di disuguaglianze :8 > < > :1 >0 22 >0 532 >0: La seconda disuguaglianza è soddisfatta perp2 <