Mechanical Engineering - Metodi Analitici e Numerici per L'Ingegneria

Fondamenti di Calcolo Numerico

Divided by topic

METODI ANALITICI E NUMERICI PER L'INGEGNERIA Dipartimento di Ingegneria Industriale Corso di Laurea in Ingegneria Meccanica Samuele D'Elia Matricola: 10663808 A.A. 2021/221 Contents I Fondamenti di Calcolo Numerico 4 1 Introduzione 41.1 Insieme Numerico al Calcolatore . . . . . . . . . . . . . . . . . . . . . . .4 1.2 Norma (matematica) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.2.1 Norma dei Vettori inRn . . . . . . . . . . . . . . . . . . . . . . . .5 1.2.2 Norma delle Matrici QuadrateRn n . . . . . . . . . . . . . . . . .5 1.2.3 Norma delle FunzioniR2 !R. . . . . . . . . . . . . . . . . . . . .5 2 Risoluzione di sistemi lineari 62.1 Metodi Diretti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .72.1.1 Risoluzione di Matrici Diagonali . . . . . . . . . . . . . . . . . . .7 2.1.2 Risoluzione di sistemi triangolari inferiori . . . . . . . . . . . . . .7 2.1.3 Risoluzione di sistemi triangolari superiori . . . . . . . . . . . . . .8 2.1.4 Fattorizzazione LU . . . . . . . . . . . . . . . . . . . . . . . . . . .8 2.1.5 Condizioni di esistenza per la fattorizzazione . . . . . . . . . . . .10 2.1.6 Pivoting (per righe) . . . . . . . . . . . . . . . . . . . . . . . . . .11 2.1.7 Fattorizzazione di Cholesky . . . . . . . . . . . . . . . . . . . . . .11 2.1.8 Risolvere sistema tridiagonale . . . . . . . . . . . . . . . . . . . . .12 2.2 Controllo degli errori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132.2.1 Norma due vs raggio spettrale . . . . . . . . . . . . . . . . . . . .14 2.2.2 Numero di Condizionamento . . . . . . . . . . . . . . . . . . . . .16 2.2.3 Stimatore dell'errore . . . . . . . . . . . . . . . . . . . . . . . . . .16 2.2.4 Residuo piccolo implica Errore piccolo? . . . . . . . . . . . . . . .17 2.3 Metodi Iterativi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .182.3.1 Convergenza, Consistenza, Stabilita . . . . . . . . . . . . . . . . .18 2.3.2 Metodi di Richardson per decomposizione Additiva . . . . . . . . .20 2.3.3 Metodi Di Richardson per de nizione di residuo . . . . . . . . . .21 2.3.4 Come scelgo la matrice di precondizionamento P? . . . . . . . . .22 2.3.5 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . .23 2.3.6 Metodo di Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . .24 2.3.7 Valore ottimale di. . . . . . . . . . . . . . . . . . . . . . . . . .25 2.4 Metodo del Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . .28 2.5 Metodo del gradiente coniugato . . . . . . . . . . . . . . . . . . . . . . . .30 2 List of Algorithms 1 Sostituzione in avanti . . . . . . . . . . . . . . . . . . . . . . . . . . . . .7 2 Sostituzione all'indietro . . . . . . . . . . . . . . . . . . . . . . . . . . . .8 3 Metodo di Richardson in forma generale . . . . . . . . . . . . . . . . . . .22 4 Algoritmo di Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .24 5 Algoritmo di Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . .25 6 Gradiente Precondizionato (P6 =I) . . . . . . . . . . . . . . . . . . . . . . .29 3 Part I Fondamenti di Calcolo Numerico 1 Introduzione Al giorno d'oggi ci sono moltissimi problemi di interesse scienti co che vengono risolti tramite l'ausilio di calcolatori, perche troppo dicili da risolvere a mano oppure perche non esiste una soluzione analitica al quesito posto. In questo testo si incontreranno continuamente entita matematiche elementari che dovreb- bero far parte del bagaglio culturale del lettore, ma il cui ricordo e talvolta appannato. Appro ttiamo di questo capitolo introduttivo per rinfrescare quelle nozioni che saranno utili nella trattazione successiva. 1.1 Insieme Numerico al Calcolatore Il calcolatore e in grado di lavorare con un insieme nito di numeri discreti essendo impossibile rappresentare su una macchina (le cui risorse sono necessariamente nite) l'in nita dei numeri reali. Foating PointFSottoinsieme diRdi dimensione nita caratterizzato da proprieta diverse. In particolare vale che: x2F,x= (1)s M e t dovesvale 0 o 1, (un numero intero positivo maggiore od uguale a 2) e labase,me un intero dettomantissala cui lunghezza t e il numero massimo di cifre memorizzabili ede e un numero intero dettoesponente. Non sempre il risultato di un'operazione matematica appartiene all'insiemeF, per- tanto la macchina arrotonda automaticamente al numero rappresentabile precedente o successivo. Possiamo de nire un errore di arrotondamento associato: E rrore di arrotondamento=k xf l(x)kk xk< 12  m dove:{x2R {f l(x) e il oating-point piu vicino adx { m( Epsilon Macchina) e il piu piccolo numero appartenente all'insiemeF(siano a; b2Ft.c.ab < )ab= 0) Minimo e massimo diFIl minimo coincide con l'Epsilon Macchinamentre il massimo e un numero dell'ordine di 10308 4 1.2 Norma (matematica) Una norma su uno spazio vettoriale reale o complessoXe una funzione: k
k:X![0;1) che veri ca le seguenti proprieta:1.kxk 0;8x2X(Positivita) 2.kxk= 0,x= 0 3.kxk=jjkxk2R(Distributiva) 4.kx 1+ x 2k  k x 1k +kx 2k (Disuguaglianza Triangolare) 1.2.1 Norma dei Vettori inRn Possiamo de nire diverse norme:ˆNorma due o euclidea: kxk =px
x= pxT xˆ Norma uno: kxk 1= jx 1j +jx 2j +


+jx nj ˆNorma1: kxk 1= max ij x ij Le precedenti norme sono tutte convergenti, ossia xn ! xallora kxn
xk N! 08A dove N e la Norma di qualsiasi ordine (1;2; :::;1). 1.2.2 Norma delle Matrici QuadrateRn n Si usa la "norma operatoriale" per cui vale la seguente relazione:kAxk N k Ak N
k xk N dove N e la Norma di qualsiasi ordine (1;2; :::;1). Prendiamo il caso N=2 (Norma Spettrale). kAk 2=q max( AT A) La formula viene dimostrata piu avanti (Teo. 2.4). 1.2.3 Norma delle FunzioniR2 !R Siaf(x 1; x 2; :::; x n) !R, le norme della funzione sono: ˆkfk 2=qR b af 2 (x)dx(Scarto quadratico medio) ˆkfk 1= maxjf(x)j NOTA: In generale queste due normeNONsono equivalenti. 5 2 Risoluzione di sistemi lineari Il sistema lineare si compone di una serie di equazioni:8 > > < > > :a 11x 1+ a 12x 2+ :::+a 1nx n= b 1 . . .. . . an1x 1+ a n2x 2+ :::+a nnx n= b n Si puo scrivere nel seguente modo:Ax= bPer il Teorema di Crameril sistema lineare ha una sola soluzione(x) se e solo se laMatrice e non-singolare:9!x,det(A)6 = 0 In linea di principio, la soluzione del sistema potrebbe essere calcolata tramite il metodo di Cramer (x i=det (A i)det (A)), con l'inversa della matrice ( A 1 ), oppure utilizzando la MEG. Il problema sussiste quando calcoliamo il costo computazionale associato a questi metodi, che richiederebbe un tempo enorme anche per risolvere un sistema di piccole dimensioni. Per questo motivo vengono utilizzati altri metodi di risoluzione molto piu ecenti. Tali tecniche sono classi cabili in due gruppi: ˆMetodi Diretti, arrivano alla soluzione esatta del sistema in un numero nito di passi. ˆMetodi Iterativi, giungono alla soluzione esatta del sistema dopo un numero in nito di passi, l'algoritmo si blocca quando raggiunge il criterio di arresto. 6 2.1 Metodi Diretti Esistono diversi metodi diretti applicati a matrice particolari (triangolari sup/inf, Chom- sky, Thomas,...) e a matrici generiche (Fattorizzazione lu, pivoting). 2.1.1 Risoluzione di Matrici Diagonali Il sistema lineare da risolvere e:Dx= bdove: D=2 6 6 6 4d 110 : : :0 0d 22: : : 0 . . .. . .. . . 0 0: : : d nn3 7 7 7 5 Essendo che le matrici diagonali sono necessariamentenon-singolariper def.(d ii6 = 0), la soluzione e immediata e la si ottiene innpassaggi: xi=1d iib i8 i= 1;2; : : : ; n 2.1.2 Risoluzione di sistemi triangolari inferiori Il sistema lineare da risolvere e:Lx= bdove: L=2 6 6 6 4l 110 : : :0 l21l 22: : : 0 . . .. ... . . ln1l n2: : : l nn3 7 7 7 5 Hp:det(L)6 = 0,l ii6 = 0 L'agoritmo utilizzato e quello della sostituzione in avanti:Algorithm 1 Sostituzione in avantix 1=1l 11b 1 fori= 2 :ndo xi=1l ii biP i1 j=1l ijx j end forIl Costo Computazionale dell'algoritmo e: n X i=1(2 n1) =


n2 7 2.1.3 Risoluzione di sistemi triangolari superiori Il sistema lineare da risolvere e:U x= bdove: U=2 6 6 6 4u 11u 12: : : u 1n 0u 22: : : u 2n . . .. . .. . . 0 0: : : u nn3 7 7 7 5 Hp:det(U)6 = 0,u ii6 = 0 L'algoritmo necessario per la risoluzione del sistema prende il nome di sostituzione all'indietro:Algorithm 2 Sostituzione all'indietrox n=1u nnb n fori=n1 : 1do xi=1u ii biP n j=i+1u ijx j end for2.1.4 Fattorizzazione LU Hp:det(A)6 = 0,a ii6 = 0!!! Il metodo consiste nel costruire due matrice L e U (L; U2Rn n ) tali che: A=LU Segue che:Ax= bLU x= bSe chiamiamo y= U x, allora si possono risolvere due equazioni lineari (triangolare supe- riore e inferiore):( Ly= b, sostituzione in avanti U x= y, sostituzione all0 indietro Il problema si e spostato su come trovare le matrici L e U.Si puo dimostrare che le matrici L e U vengono calcolare con ilmetodo di eliminazione di Gauss (MEG): 1.ad ogni iterazione di riarragia A al ne di farla diventare triangolare superiore!U 2.gli elementi "scartatori" costituiscono la matrice tri. inferiore con diagonale princi-pale costituita da elementi unitari!L 8 esempio: A (0) =2 41 2 1 2 01 1 1 53 5 + l21=a (0) 21a (0) 11= 2; l 31=a (0) 31a (0) 11= 1 A(1) =2 41 2 1 043 0 3 63 5 + l32=a (1) 32a (1) 22= 34 ; A(2) =2 41 2 1 043 0 0154 3 5 Le matrici L e U sono: L=2 41 0 0 2 1 0 134 13 5;U=2 41 2 1 043 0 0154 3 5 La sistassi della fattorizzazione semplice inMatlabe:1 [ L , U ] = l u ( A ) ;Essa produce una matrice L che e uguale a P*M, dove M e triangolare inferiore e P e una matrice di permutazione generata dal pivoting per righe. Il comandoluattiva automati- camente il pivoting per righe. In particolare, quando la matrice e memorizzata in formato sparse, la pivotazione per righe e eettuata solo quando viene generato un elemento pivot nullo (o molto piccolo). Il Costo Computazionale della fattorizzazione LU e circa:23 n 3 + 2qn2 n3 dove q e il numero delle equazioni da risolvere. 9 2.1.5 Condizioni di esistenza per la fattorizzazione Teorema 2.1.(Condizione necessaria e suciente) Sia data una matriceA2Rn n Non-Singolare, esiste una sola fattorizzazione LU se e solo se le sottomatrici di Nord-OvestA isono tutte Non-Singolari: 9!A=LU,det(A i) 6 = 08 i= 1 ;2; : : : n Commenti: 1.9!, l'unicita si ha solo con la convenzionel 11= l 22=


=l nn= 1 2.Le sottomatrici nord-ovest sono tutte quelle sottomatrici quadrate il cui termine 11coincide con quello della matrice principale Condizione dicile da veri care perche calcolare il determinante e molto costoso dal punto di vista computazionale. Teorema 2.2.(Condizioni sucenti (e pratiche) per la fattorizzazione) Sia data una matriceA2Rn n Non-Singolare, esiste una sola fattorizzazione LU se UNAdel le seguenti aermazioni e veri cata: 1.A e simmetrica e de nita positiva 2.A e a dominanza diagonale stretta per le righe 3.A e a dominanza diagonale stretta per le colonne Commenti: Matrice a dominanza diagonale stretta per le righe si intende che ssata una riga il modulo dell'elemento appartenente alla diagonale principale e maggiore al modulo della somma degli altri elementi: ja iij >X j6 =ij a ijj 8 i Stesso discorso vale per le colonne. 10 2.1.6 Pivoting (per righe)  E utilizzato nel caso in cui, al ne di determinare la fattorizzazione LU di una matrice, sia presente un termine pivotiale NULLO. Il meccanismo alla base di questo metodo e quello di eettuare delle permutazioni su A, ossia uno scambio delle righe, per poi ritornare all'algoritmo di Gauss. Matrice di permutazione= matrice che si ottiene scambiando alcune righe della matrice identita. Esempio: Iperm=2 41 0 0 0 0 1 0 1 03 5Il pivoting puo essere utilizzato per due principali motivi: ˆnecessita, non esistono altri modi per risolvere il sistema lineare ˆaccuratezza, nel caso di a(0) ii6 = 0 ma di valore piccolo(es. 10 20 ). In questa casistica il pivoting tende a propagare minore incertezza rispetto alla fattorizzazione senza pivoting. La sintassi della fattorizzazione con pivoting inMatlabe:1 [ L , U , P ] = l u ( A ) ;Essa restituisce le matrici L, U e la matrice di permutazione (P). 2.1.7 Fattorizzazione di Cholesky E applicata a tutte le matriciA2Rn n simmetriche e de nite positive. Tale processo permette di determinare un'unica matrice R tale che: A=R
RT dove R e una matrice triangolare inferiore con elementi positivi sulla diagonale (a ii2 R; a ii> 0; a ii=pa iipa ii) 11 2.1.8 Risolvere sistema tridiagonale Matrice tridiagonale= gli unici elementi che possono essere non nulli appartengono alla diagonale principale ed alla prima sopra e sotto-diagonale. T=2 6 6 6 6 4t 11t 12: : : 0 t21t 22. . . . . .. .. 0: : : t nn3 7 7 7 7 5 Se la fattorizzazione LU di A esiste, i fattori L e U sono due matrici bidiagonali (inferiore e superiore, rispettivamente) della forma L=2 6 6 6 6 41 0 : : :0 21. . . . . .. ..0 0: : : n13 7 7 7 7 5; U=2 6 6 6 6 4 1 2: : : 0 0 2. . . . . .. .. n 0: : :0 n3 7 7 7 7 5 Si osserva immediatamente che le matrici L e U sono piu semplivi rispetto al caso tradizionale e di conseguenza si avranno gradi di sempli cazione sull'algoritmo. Il Costo Computazionale di questa procedura, conosciuta comeAlgoritmo di Thomas, e:8n1n NOTA: Per matrici Sparsee Strutturatenon si fa mai pivoting. ˆMatrice Sparsa, e una matrice i cui valori sono quasi tutti uguali a zero. ˆMatrice Strutturata, e una matrice in cui a ij6 = 0 sono messi in "punti strategici" (es. diagonale). 12 2.2 Controllo degli errori Il piu grande problema legato ai metodi diretti e sensa dubbio la propagazione dierrori macchina (), a causa degli errori di arrotondamento il prodotto LU non ritorna esat- tamente la matrice A (il pivoting in certi casi, non tutti, consente di contenere questo genere di errori). Si parla di accuratezza della soluzione in due fasi della risolzione del problema: 1.Durante la fattorizzazione si dice che essa e accurata se P ALU0 2.Durante il calcolo della soluzione ^xsi dice che essa e accurata se er=k x ^ xk k xk  0 1 Oimplica necessariamente2 O? Purtoppo NO. Esempio: ( Matrice di Hilbert) H=2 6 6 6 6 41 12 13 : : : 12 13 13 . .. . . .. ..3 7 7 7 7 5 Il sistema lineare da risolvere e il seguente:H x= bn doveb ne stato scelto in modo tale che che la soluzione esatta del sistema sia x= (1;1; : : : ;1)T . De niamo: ˆErrore di fattorizzazione, Rn= LUP A ˆErrore soluzione, En=k x ^ xk k xk 13 Figure 1: E nin linea continua e R nin linea trattegiata Prima di riportare i diversi stimatori dell'errore bisigna fare una piccola trattazione teorica sulla norma e raggio spettrale per andare a de nire ilNumero di Condizionamento, ossia uno strumento molto utile per paragonare l'errore relativo con il residuo relativo. 2.2.1 Norma due vs raggio spettrale De niamo la norma due e il raggio spettrale per poi confrontalrli:ˆNorma due= SiaA2Rn n , la norma e la piu piccola costante tale chekAxk  kAk 2k xk , ossia kAk 2= max x6 =0k Axk k xk = max kxk =1k Axk (1) ˆRaggio spettrale= e il massimo modulo degli autovalori diA. (A) = max i=1;:::;nj  i( A)j(2) In generale, possiamo aermare che: (A) kAk 2 Nel caso particolare diAsimmetrica e de nita positiva: (A) =kAk 2 Dimostrazione.Per de nizione di autovalore e autovettore della matriceAvale la seguente uguaglianza:v= Av+ jj
kvk =kvk =kAvk  k Ak 2
k vk + jj
kvk  k Ak 2
k vk 14 Teorema 2.3. (Quozioziente di Rayleigh) SiaBuna matrice simm. def>0, al lora vale che 0<  min( B)Bx
xk xk 2  max( B) =(B) dove:R(B; x) = Bx
xk xk 2 e dettoQuoziente di Rayleigh. Caso Particolare:Bv= v+ R(B; v) = Bv
vk vk 2=v
vk vk 2=  Teorema 2.4.SiaAuna matrice qualunque2Rn n , al lora kAk 2=q (AT A) Dimostrazione. kAk2 2= max x6 =0k Axk 2k xk 2= max x6 =0Ax
Axk xk 2=" Ax
y= x
AT y# = maxx6 =0A T Ax
xk xk 2=  max( AT A)Notare che AT A=Be una matrice simmetrica de nita positiva. 15 2.2.2 Numero di Condizionamento In norma 2: K2( A) =kAk 2k A 1 k2=s max( AT A) min( AT A) Spettrale:K(A) = max( A) min( A) NOTA: seAe simm. def>0 allora le due formule coincidono. 2.2.3 Stimatore dell'errore Quando si risolve un sistema lineare del tipoAx= bin realta si sta trovando la soluzione "esatta" di unsistema perturbatodella forma: (A+A)(x+ x) = b+ bdove A; bsono rispettivamente la matrice/vettore di perturbazione. Prima di continuare facciamo le seguenti osservazioni:ˆHp:A0 ˆLa soluzione calcolata ( ^x) NON COINCIDE con la soluzione esatta ( x) ma: ^ x= x+ xGli stimatori dell'errore sono: 1.Errore Relativo, quanto si discosta la soluzione calcolata (^ x) dalla soluzione reale (x), supponendola di conoscerla (Condizione quasi mai vera): erel=k x ^ xk k xk = k xk k xk 2.Residuo Relativo, stima dell'errore senza conoscere la vera soluzione ( x): rrel=k b A^ xk k bk = k bk k bk 16 2.2.4 Residuo piccolo implica Errore piccolo? NON NECESSARIAMENTE. Dimostrazione. C onto1:= ( (A+A)(x+ x) = b+ bAx= b+ kxk =kA 1 bk  kA 1 k2k bk +k xk  k A 1 k2k bk C onto2:= ( kbk =kAxk kAxk  k Ak 2k xk +1 k xk  k Ak 2k bk CombinandoC onto1con C onto2abbiamo che: kxk k xk  k A 1 k2k Ak 2k bk k bk +e rel K 2( A)
r relSe nel sistema lineare non ci sono uttuazioni dei coecenti della matrice ( A0) possiamo distnguere due casistiche: 1.K 2( A)"piccolo"(100 101 ), la matrice A e "Ben Condizionata" ed a piccoli errori sui dati corrispondono errori dello stesso ordine di grandezza sulla soluzione. 2.K 2( A)"grande"(109 ), la matrice e "Mal Condizionata" e potrebbe accadere che a piccole perturbazioni sui dati corrispondono grandi errori sulla soluzione. SeAe "Mal Condizionata" i Metodi Diretti non sono accurati. 17 2.3 Metodi Iterativi Un metodo iterativo per la risoluzione del sistema lineare conA2Rn n eb2 Rn consiste nel costruire una successione di vettorix( k) ; k0 diRn checonvergealla soluzione esatta x, ossia tale che lim k!1x( k) =x(3) per un qualunque vettore inizialex(0) 2Rn . Per realizzare questo processo una possibile strategia e quella di de nire ricorsivamente x( k+1) =Bx( k) +g; k 0;(4) essendo B una matrice opportuna (dipendente da A) egun vettore opportuno (dipendente da A e dab). 2.3.1 Convergenza, Consistenza, Stabilita1.Un metodo numerico si diceconvergenteperx(0) se lim k!1x( k) =x2.Un metodo e consistentese la soluzionexe un punto ssodell'algoritmo iterativo se x(0) =x= )x( k) =x; 8k Lacondizione di consistenzapuo essere scritta anche come x= Bx+ g(5) NOTA: generalmente i metodi sono consistenti per costruzione. 3.Un metodo estabilesex( k) si avvicina axad ogni iterazione. 18 Teorema 2.5. (Consistenza implica convergenza) Se l'algoritmo e consistente al lora converge8x(0) e (B)0al lora il metodo converge al piuniterazioni (in aritmetica esatta). 2.il metodo del gradiente coniugato converge axper ogni scelta di x(0) 3.L'errore e de nito come:ke( k) k ~ dK ke(0) k dove~ ddipende dapK (P 1 A) La sintassi del Gradiente coniugato inMatlabe:1 [ x , F L A G , R E S , I T E R ] = p c g ( A , b , T O L , M A X I T , P ) ;30