Mathematical Engineering - Matematica Numerica

Exercise - 07 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 7 - Soluzione Sistemi sovradeterminati e Metodi iterativi per sistemi lineari (1) Esercizio 1 Sia data la seguente matrice che possiede rango pieno per ogni >0 A =2 41 1 0 03 5; e i vettorix ex= [1 ;1]T eb= Ax ex. Vogliamo calcolare x, la soluzione del seguente problema ai minimi quadratikAxbk2 2= min y2R2k Aybk2 2; benché sappiamo che la soluzione esatta di tale problema è propriox=x ex. a)Per= 10 i ; i= 1;2;3;4calcolarex=x 1risolvendo il sistema delle equazioni normali ovvero AT Ax 1= AT b: Per ogni valore di, calcolare l'errorekx 1 x exk 2e visualizzarlo in funzione di su un graco in scala logaritmica, insieme all'andamento del numero di condizionamentoK 2(AT A). b)Per gli stessi valori didel punto precedente, calcolarex=x 2sfruttando la fattorizzazione~ Q~ Rridotta diA, ovvero applicando la relazione, x2=~ R 1 ~ QT b: Le matrici~ Qed~ Rtali cheA =~ Q~ Rsi ricavano in Matlab attraverso il comando[Q,R]=qr(A,0). Calcolare l'errorekx 2 x exk 2e visualizzarlo in scala logaritmica sul graco precedente. Che cosa si osserva ? Soluzione 1 Attraverso i seguenti comandi Matlab, risolviamo il problema ai minimi quadrati sia attraverso il metodo delle equazioni normali, sia attraverso la fattorizzazione~ Q~ R. Inne, tracciamo il graco degli errorikx 1 x exk 2e kx 2 x exk 2in funzione di .for i = 1:4a(i) = 10^(-i); A = [1 1; a(i) 0; 0 a(i)]; AA = A *A;xex = [1 1]; b = A *xex;%soluzioneattraversoleequazioninormali K(i) = cond(AA); bb = A *b;x1 = AA(A) *b;err1(i) = norm(xex-x1); %soluzioneattraversolafattorizzazioneQR [Q,R] = qr(A,0); x2 = R \ (Q *b);1 err2(i)=norm(xex-x2); end figure; loglog(a,K,b,a,err1,r,a,err2,g,LineWidth,2); grid; xlabel(\epsilon) legend(K_2(A^TA),errconeq.normali,errconQR,Location,East) Otteniamo il risultato della Figura 1, dove si osserva che l'algoritmo di calcolo basato sulla fattorizzazione ~ Q~ Rè notevolmente più stabile rispetto alla propagazione degli errori di arrotondamento. In altre parole, il metodo basato sulle equazioni normali risulta essere instabile poiché la matriceAT Atende ad essere molto mal condizionata.Figure 1  Numero di condizionamentoK 2(AT A)ed errorikxx exk 2nel caso di risoluzione mediante eq. normali e mediante fattorizzazione QR, in funzione di. Esercizio 2 Si consideri il sistema lineareAx=bdove la matriceAè denita da A =2 42 1 0 1 21 01 23 5eb=2 41 0 13 5: Si consideri ora il seguente metodo iterativo :Lx( k+1) = Lx( k) +(bAx( k) ); k0;(1) dove >0è un parametro reale e positivo eLè la parte triangolare inferiore diA: L =2 42 0 0 1 2 0 01 23 5: a)Vericare se il metodo (1) è consistente.210-410-310-210-110-2010-1510-1010-51001051010K2(ATA)err con eq. normalierr con QR Figure 2  Funzionij1j(linea continua) ej1=2j(linea tratteggiata) b)Riscrivere il metodo (1) nella formax( k+1) = Bx( k) +z , k0, esplicitando la matriceB . Stabilire per quali valori del parametro reale >0il metodo (1) converge. c)Determinare il parametro di convergenza ottimale opt. Soluzione 2a)Osserviamo immediatamente che, dettaxla soluzione esatta del sistema lineare, si ha Lx= Lx+(bAx): Quindi il metodo (1) è consistente per ogni valore del parametro. b)Si ha Lx( k+1) = Lx( k) +(bAx( k) ) = (LA)x( k) +b; e osservando cheLè invertibile poiché det(L)6 = 0, si ha x( k+1) = L 1 (LA)x( k) +L 1 b = (IL 1 A)x( k) +L 1 b: Poiché L 1 =2 41 =2 0 0 1=4 1=2 0 1=8 1=4 1=23 5; si ottiene la matrice di iterazione seguente (dipendente da) : B=2 41  =2 0 0 13=4=2 0=8 13=43 5: Gli autovalori della matriceB sono le radici del polinomio (1)" 2 2 134  + 134  2  216 # = 0 ovvero : 1= 1 , 2= 1 e 3= 1 =2. Allora(B ) = max fj1j;j1=2jg, >0, si veda il graco riportato in Figura 2.30 1 2 0 1 2 3 4 Si vede che l'ordinata del punto di intersezione tra le due curve j1j,j1=2jper >0, fornisce il valore minimo del raggio spettrale (quest'ultimo si ottiene a sua volta come il massimo fra le due curve). Tale punto ha coordinate(4=3;1=3), da cui si ottiene : 0<