Mathematical Engineering - Matematica Numerica

Exercise - 04 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 4 - Soluzione Propagazione dell'errore in sistemi lineari Esercizio 1 Si consideri la matrice di HilbertA2Rn n A =2 6 6 6 6 6 41 12 13 14 : : : 12 13 14 : : : 13 : : : : : :3 7 7 7 7 7 5; a ij=1i +j1: a)Calcolare pern= 1;2; : : : ;10la matriceA(scrivendo un'opportunafunctionMatlab oppure uti- lizzando il comandohilb) e il suo numero di condizionamentoK 2(A) mediante il comandocond. Salvare i numeri di condizionamento calcolati in un vettore di 10 componenti. b)Visualizzare su un graco in scala semilogaritmica (lineare nelle ascisse e logaritmica nelle ordinate)l'andamento diK 2(A) in funzione din, utilizzando il comandosemilogy. Dedurre cheK 2(A) si comporta comeen . c)Si costruisca un vettore colonnax exunitario di ncomponenti, utilizzando il comandoonese si calcoli il vettoreb= Ax ex. Per ogni n= 1;2; : : : ;10, si risolva, quindi, il sistema lineareAx=bmediante il comandone si calcoli l'errore relativo" r=k xx exkk x exk. Si visualizzi tale errore su un graco in scala semilogaritmica e si deduca che esso ha lo stesso comportamento del numero di condizionamento diA. d)Si mostri numericamente che il rapporto tra" re K 2(A) , per ognin= 1;2; : : : ;10, è dell'ordine dell'unità diround-o. e)Denito ilresiduor=bAx, vale la seguente stimaa posterioridi"r : kxx exkk x exk K(A)k rkk bk: (1) Per ognin= 1;2; : : : ;10, si calcoli la norma del residuokrke del residuo relativokrk=kbk. Si noti che il residuornonè esattamente nullo, a causa degli errori diround-o. f )Si verichi che la norma del residuo ottenuta è dell'ordine dell'unità diround-oper ognin. g)Si verichi la stima a posteriori (1). h)Si giustichino i risultati ottenuti al punto (d) sulla base dei punti precedenti. Soluzione 1a)Calcoliamo le quantità richieste dall'esercizio mediante le seguenti istruzioni Matlab, all'interno di uncicloforpern= 1;2; : : : ;10: help hilbHILB Hilbert matrix. HILB(N) is the N by N matrix with elements 1/(i+j-1), which is a famous example of a badly conditioned matrix.1 KA = zeros(1,10); err = zeros(1,10); nn = [1:1:10]; for n = 1:10A = hilb(n); KA(n) = cond(A); xex = ones(n,1); b = A*xex; x = A\b; err(n) = norm(x-xex)/norm(xex); end b-c)Il condizionamento della matriceAè salvato nel vettoreKA, e possiamo visualizzare il suo andamento rispetto adntramite le seguenti istruzioni (Figura 1) % Grafico semilogaritmico: semilogy(nn, KA, 'r'); % Creiamo una griglia e dei nomi per gli assi coordinati: grid; xlabel('n'); ylabel('KA(n)');Figure 1  Andamento del numero di condizionamentoK 2= K 2(A) della matrice di Hilbert e dell'errore relativo" r=k xx exkk x exkin funzione di n; scala lineare in ascissa e logaritmica in ordinata Visualizziamo ora l'andamento del numero di condizionamento e dell'errore relativo" rsullo stesso graco (Figura 1) : semilogy(nn, KA, 'r--o', nn, err, 'b-s'); grid; xlabel('n'); legend('KA(n)', 'err(n)'); Si nota che il numero di condizionamento ha un andamento lineare nel graco ; questo signica dipen- denzaesponenzialediK 2(A) dan. Si nota inoltre che l'errore relativo ha un comportamento analogo al numero di condizionamento, anche se molto più piccolo. Questo conferma il fatto che, più il condi- zionamento della matrice è grande, più gli errori diround-ovengono amplicati nella risoluzione del sistema lineare portando a un errore relativo sempre maggiore, in accordo con la stima kxkk xk K (A)1 K(A)kAk=kAk kbkk bk+ k Akk Ak :(2) d)Calcolando il rapporto" r=K 2( A), si trova21 2 3 4 5 6 7 8 9 10 10−20 10−15 10−10 10−5 100 105 1010 1015    K 2 "r n max(err./KA) ans =0.3283e-16 e si vede come esso sia eettivamente dell'ordine dell'unità diround-o: eps ans =2.2204e-16 e)Aggiungiamo le righe per il calcolo del residuo normalizzato e del termine noto della stima (1).err=zeros(1,10); res=zeros(1,10); stima=zeros(1,10); for n=1:10A=hilb(n); KA(n)=cond(A); xex=ones(n,1); b=A*xex; x=A\b; err(n)=norm(x-xex)/norm(xex); res(n)=norm(A*x-b)/norm(b); stima(n) = KA(n)*res(n); end f-g-h)Procediamo alla verica dei punti richiesti, aggiungendo al graco precedente l'andamento al variaredindel residuo normalizzato e della maggiorazione (1). Si ottiene il graco in gura 2. figure; semilogy(nn,KA,'-*r','DisplayName','cond'); hold on; semilogy(nn,err,'-ob','DisplayName','err-rel'); semilogy(nn,res,'-^k','DisplayName','res-norm'); semilogy(nn,stima,'-sg','DisplayName','stima'); grid; xlabel('n'); legend show Osserviamo dai graci ottenuti che eettivamentekrk=kbksi mantiene moderatamente costante, in- torno a10 16 , ovvero l'ordine di grandezza dell'unità diround-o. Abbiamo inoltre visualizzato la stima superioreK(A)krk=kbk: vediamo come essa segua l'errore relativo" re rimanga sempre mag- giore ad esso, il che conferma la stima a posteriori (1). Tutto ciò conferma quanto osservato al punto (d) : infatti, essendokrk=kbkcirca costante, la stima a posteriori prevede che, nel caso peggiore, l'errore" raumenti come K(A); questo è quello che eettivamente accade. Un discorso a parte merita il cason= 2. In tal caso, le operazioni di arrotondamento e di aritmetica oating point fanno sì chex ex6 =x, come ci si aspetta, e che quindi l'errore relativo non sia nullo. Tuttavia, il residuor= Axbè nullo ; si tratta di un caso sfortunato in cui gli errori compiuti nelle operazionib= Ax ex, x= Anber= Axbsi elidono. Si dovrebbe quindi ricorrere ad un'analisi dell'errore più generale, in cui si considerino anche gli errori di arrotondamento e cancellazione intro- dotti dal MEG (che possono essere descritti ad esempio tramite una perturbazioneAsulla matrice del sistema lineare). Si noti inoltre che, utilizzando dierenti versioni di Matlab, è possibile che il residuo si annulli anche per altri valorin. Ciò è dovuto a possibili dierenze nell'implementazione delle funzioni Matlab necessarie a svolgere l'esercizio.3 Figure 2  Andamento diK(A)(asterischi), dell'errore relativo" r=k xx exkk x exk(cerchi), del residuo relativo krkk bk(triangoli) e della stima a posteriori K(A)k rkk bk(quadrati) in funzione di n; scala lineare in ascissa e logaritmica in ordinata. Esercizio 2 SiaA2Rn n una matrice simmetrica denita positiva, indichiamo con ie v i, i= 1; : : : ; ni suoi autovalori e autovettori. Si verichi che sexè la soluzione del sistema lineareAx=b, allora x=n X i=1c i iv i; dovec iè l' i-esima componente del vettorebnel suo sviluppo rispetto alla base ortonormale degli autovettori diA. Si consideri poi il seguente sistema lineareAx=b 1001 1000 1000 1001  x1 x2 = b1 b2 : Si osservi che la matriceAè mal condizionata avendo autovalori 1= 1 , 2= 2001 . Datob= [2001 2001]T , vericare che una piccola perturbazioneb= [1 0]T produce una grande variazione nella soluzione (si utilizzi il risultato dell'esercizio precedente). Analogamente, si verichi che una piccola perturbazionex= [0:001 0]T nella soluzione produce una piccola variazione nel termine notob. Soluzione 2 PoichéAè una matrice simmetrica, grazie al teorema spettrale sappiamo che i suoi autovalori isono tutti reali e i corrispondenti autovettoriv iformano una base ortogonale di Rn (e possono essere normalizzati in modo tale chevT jv l=  j l, dove  j lè il simbolo di Kronecker). Sia dunque b=P n j=1c jv jlo sviluppo del termine notobrispetto alla base degli autovettori diA. Esprimendo anche la soluzionexnella base degli autovettori diA, otteniamo Ax=An X i=1x iv i=n X i=1c iv i e dunque, ricordando cheAv i=  iv iper 1; : : : ; n, n X i=1( x i i c i) v i= 0:(3)41 2 3 4 5 6 7 8 9 10 10-20 10-15 10 -10 10-5 100 105 10 10 1015 cond err-rel res-norm stima Poiché gli autovettori v isono linearmente indipendenti, avremo che x i i= c i, i= 1; : : : ; n, da cui x=n X i=1c i iv i: Gli autovettori della matriceAproposta nel presente esercizio sonov 1= [1 1]T (corrispondente all'auto- valore 1= 1 ) ev 2= [1 1]T (corrispondente a 2= 2001 ). Sianob= [2001 2001]T eb= [1 0]T . Allora, b+b= 2001 2001 + 1 0 = 2001v 2+12 ( v 1+ v 2) =12 v 1+40032 v 2: Esprimendo la soluzione perturbatax+xcome combinazione lineare degli autovettori otteniamo x+x=12 1 v 1+40032 2001 v 2=12 v 1+40034002 v 2' 1:5 0:5 : La perturbazionebproduce quindi un errorex'[0:50:5]T sulla soluzione esattax= [1 1]T del problema non perturbato. Consideriamo ora il vettoreb= [1;1]T come termine noto del sistema lineare, la soluzione esatta èx= [11]T . Esprimiamo la soluzione perturbatax+xrispetto alla base degli autovettori : x+x= 1:001 1 =2 :0012 v 1+0 :0012 v 2; da cui ricaviamoc 1= 2 :001=2ec 2= 2 :001=2. Dunque b+b= 2:001 0 eb= [1:001 1]T . Una possibile interpretazione delle conclusioni raggiunte è la seguente. Il sistema lineare considerato in questo esercizio può essere ottenuto dall'analisi di un'asta rigida collegata nella parte centrale ad una molla di rigidezzaK 3= 4000 e connessa alle estremità a due molle di rigidezzaK= 1(si veda la Figura 3) ; alle estremità dell'asta sono applicate due forzeb 1e b 2. Indicando con x 1e x 2lo spostamento verticale alleFigure 3  Asta rigida collegata ad un sistema di tre molle estremità dell'asta e imponendo l'equilibrio delle forze (occorre notare che lo spostamento nel centro dell'asta è pari a(x 1+ x 2) =2), otteniamo le seguenti equazioni : 8 > > < > > :K x 1+ K 3x 1+ x 24 = b 1; K x2+ K 3x 1+ x 24 = b 2; ovvero il sistema lineare considerato. Se le forzeb 1et b 2sono equilibrate (ad esempio b 1= 2001 ,b 2= 2001 ), allora delle piccole perturbazionibcausano degli spostamenti signicativi della barra (xgrande). Viceversa, se le forze non sono equilibrate (ad esempiob 1= 1 ,b 2= 1), anche imponendo delle perturbazionibelevate otteniamo degli spostamentixmolto contenuti.5x1x 2 b1 b2