1.CONSORZIO NETTUNO POLITECNICO DI TORINO DIPLOMI UNIVERSITARI TELEDIDATTICI Eugenio BRUSA Cristiana DELPRETE Paolo GAY Tutorato di CALCOLO NUMERICO Settembre 2001 1. Questaraccoltadiesercizienote,prodottaadusointerno,vieneutilizzataperitutorati dellinsegnamento diCalcolo Numerico dei Corsi di Diploma Universitario Teledidattico in Ingegneria Automatica e Informatica, Ingegneria Elettrica, Ingegneria Elettronica e Ingegneria Meccanica. Nel corso dellastesura si fatto riferimento al testo G. MONEGATO,100 Pagine di ... Elementi di Calcolo Numerico, Levrotto & Bella, Torino, 1996. Il materiale proposto non vuole in alcun modo sostituire il libro di testo, a cui si fa riferimento per quanto riguarda lebasi teoriche e limpostazione formale dellinsegnamento, ma viene fornito esclusivamente a titolo di supporto didattico volto a facilitare lapproccio alla materia. Eugenio Brusa Cristiana Delprete Paolo Gay Torino, settembre 2001 Vietate la riproduzione e qualsiasi forma di commercializzazione. 1.TUTORATO DI CALCOLO NUMERICO Contenuti Fondamenti del calcolo scientifico. Sistema di numerazione e sistemi aritmetici. Errore assoluto ed errorerelativo.Cifredecimalicorretteecifresignificativecorrette.Erroridiarrotondamento. Troncamentosempliceearrotondamentosimmetrico.Cancellazionenumerica.Esercizisvolti. Esercizi proposti. Richiamisulcalcolomatriciale.Matricieoperazionitramatrici.Matricetrasposta.Matrice simmetrica. Determinante. Matrice simmetrica definita o semidefinita positiva. Matrice non singolare. Rango. Matrice inversa. Norma. Esercizi svolti. Esercizi proposti. Sistemi di equazioni lineari. Generalit sui metodi di soluzione: diretti e iterativi. Soluzione di sistemi triangolare superiori e inferiori. Algoritmo Backsost. Metodo di Gauss: procedura, decomposizione GA=U, fattorizzazione PA=LU, algoritmo Factor e algoritmo Solve. Metodo di Choleski: procedura ealgoritmo.Esercizisvolti.Eserciziproposti.Cennisulcondizionamento dei sistemi lineari. Metodi iterativi:proceduraecennisullaconvergenza.MetododiJacobi.MetododiGauss-Seidele algoritmo Gseidel. Cenni sul metodo SOR. Autovalorieautovettori.Autovaloridimatrici.Metododellepotenze:autovaloredimodulo massimoeminimo.AlgoritmoPow.Metododellepotenzeinverse:autovaloredinota approssimazione. Algoritmo Invpow. Cenni sulle trasformazioni di similitudine. Cenni sul metodo QR per il calcolo di tutti gli autovalori. Approssimazione di funzionie dati sperimentali. Classi delle funzioni interpolanti. Criteri di scelta dellafunzioneinterpolante.Interpolazionepolinomiale:metododiLagrangeemetododiNewton. AlgoritmoDifdivealgoritmoInterp.Interpolazionepolinomialeatratti.Splineesplinecubica. Algoritmo Spline e algoritmo Valspl. Approssimazione di dati. Minimi quadrati. Regressione lineare. Esercizi svolti. Esercizi proposti. Equazioninonlineari. Metodo di bisezione e algoritmo Bisez. Regula falsi. Metodo di Newton o delletangenti.MetododellesecantiealgoritmoSecant.Confrontotraidiversimetodi.Criteridi arresto. Esercizi svolti. Esercizi proposti. Quadratura.Generalit.Formulediquadraturadibase(Newton-Cotes): rettangolo, punto medio, trapezio, di Simpson. Formule di quadratura Gaussiane. Formule di quadratura composte. Tecniche di quadratura automatica. Esercizi svolti. Esercizi proposti. Equazionidifferenzialiordinarie.Generalit.CondizionidiLipschitz.Metodiesplicitiapasso singolo e a passo multiplo. Metodo di Eulero. Metodi di Runge-Kutta. Temi desame. fornito il testo di tre prove di esonero. Testo di riferimento G. MONEGATO,100 Pagine di Elementi di Calcolo Numerico, Levrotto & Bella, Torino, 1996. 1.SOMMARIO SISTEMI NUMERICI Sistema di numerazione posizionale1 Sistema aritmetico intero1 Sistema aritmetico reale1 Errore assoluto ed errore relativo2 Cifre decimali corrette2 Cifre significative corrette3 Errori di round-off3 Troncamento semplice3 Arrotondamento simmetrico3 Precisione di macchina4 Cancellazione numerica4 Considerazioni finali5 Esercizi svolti5 Esercizi proposti9 OPERAZIONI TRA MATRICI Operazioni elementari10 Matrice trasposta11 Matrice simmetrica11 Determinante di matrice12 Matrice simmetrica definita (o semidefinita) positiva13 Criterio di Sylvester13 Matrice a diagonale dominante14 Vettori linearmente indipendenti14 Matrice non singolare14 Rango di matrice14 Matrice inversa14 Norma di vettori e matrici15 Norma di vettore15 Norma di matrice16 Compatibilit di norma16 Esercizi svolti17 Esercizi proposti18 SISTEMI DI EQUAZIONI LINEARI Introduzione20 Sistemi triangolari21 Metodo di Gauss22 Osservazioni sul Pivoting22 Decomposizione GA=U23 Fattorizzazione PA=LU23 Utilit della fattorizzazione24 Metodo di Choleski25 1.Condizionamento di sistemi lineari26 Metodi di soluzione iterativi27 Convergenza28 Metodo di Jacobi28 Metodo di Gauss-Seidel29 Metodo SOR29 Algoritmi30 Esercizi svolti32 Esercizi proposti36 AUTOVALORI DI MATRICI Introduzione38 Metodo delle Potenze38 Metodo delle Potenze Inverse41 Trasformazioni di similitudine42 Metodo QR43 APPROSSIMAZIONE DI DATI E FUNZIONI Introduzione44 Classi di funzioni approssimanti44 Criteri di scelta della funzione45 INTERPOLAZIONE DI DATI Criteri di interpolazione polinomiale di dati46 Metodo di Lagrange: polinomi fondamentali di Lagrange46 Esercizi svolti48 Metodo di Newton: differenze divise48 Esercizi svolti50 Algoritmi52 Esercizi proposti53 Interpolazione polinomiale a tratti54 Interpolazione lineare a tratti54 Interpolazione con funzioni splines54 Definizione di spline55 Numero di incognite per definire la spline55 Numero di equazioni disponibili55 Equazioni supplementari55 La spline cubica naturale56 Esercizi svolti57 Algoritmi59 Esercizi proposti61 APPROSSIMAZIONE DI DATI Approssimazione di dati62 Regressione lineare63 Minimi quadrati: calcolo dellapprossimazione63 1.Elaborazione del risultato63 Interpretazione geometrica del risultato64 Errore quadratico64 Esempio: trovare una retta approssimante65 Esercizi proposti66 EQUAZIONI NON LINEARI Introduzione67 Metodo di bisezione67 Generalit su Regula Falsi e metodi delle tangenti e secanti68 Regula Falsi69 Metodo delle tangenti69 Metodo delle secanti70 Confronto tra i diversi metodi71 Criteri di arresto71 Esercizi svolti72 Esercizi proposti73 CALCOLO DI INTEGRALI MEDIANTE FORMULE DI QUADRATURA Introduzione73 Formule di quadratura di base (Newton-Cotes)76 Formule di quadratura composte77 Polinomi ortogonali79 Formule di quadratura Gaussiane80 Formule di quadratura automatiche80 Esercizi proposti81 EQUAZIONI DIFFERENZIALI ORDINARIE (ODE) Introduzione82 Soluzione numerica di ODE83 Condizione di Lipschitz83 Metodi one-step e multi-step83 Metodo di Eulero84 Metodi di Runge-Kutta84 Convergenza di un metodo one-step e ordine di convergenza85 TEMI DESAME Esonero del 05 maggio 199487 Esonero del 11 luglio 199588 Esonero del 14 luglio 199789 E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO1 1.SISTEMI NUMERICI Riferimento al testo: Cap. I SISTEMA DI NUMERAZIONE POSIZIONALE Qualunque numero intero1 > Npu essere scelto come base del sistema di numerazione e, in sistema posizionale con base N, qualsiasi reale si pu scrivere come ( ) ( ) t + + + t 1 0 1 111 0 111a . a a a a = N a . + N+a a N +a N a = am mmmmm N dove le cifre ai sono numeri interi appartenenti allintervallo [0, 1 N ]: 1 0 N ai Esempi 1)( )2 1 1 21010 7 10 2 5 10 4 10 2 27 . 245 + + + + 2)( ) ( )106 5 4 3 2 1 12921875 . 2 2 1 2 1 2 0 2 1 2 1 2 1 0 2 1 111011 . 10 + + + + + + + Isistemiaritmeticidiuncalcolatoresonoaprecisione finita (spazio di memoria finito) ed effettuano operazioniconprecisionefinita(operazionidimacchina: , ...); lebasiutilizzatesonoN = 2 o N = 16. I numeri di macchina sono numeri rappresentabili esattamente nello spazio di memoria disponibile. SISTEMA ARITMETICO INTERO un sottoinsieme finito dellinsieme infinito degli interi Z:I Z ( ) { }max max1 0 1 , i ... ... i = l N I con1max lN i . Nello spazio di memoria si usano 1 bit per il segno e 7 bit per il numero (l=7). Si ha Overflow o Underflow sugli interi quandoa Z , maI , cio a>imax o a0 e quindi C definita positiva; se y=0 allora xTAx>0 e quindi A definita positiva. (cvd) ESERCIZI PROPOSTI 1.Verificare che il prodotto di 2 matrici triangolari inferiori una matrice triangolare inferiore. 2.SvilupparelalgoritmopercalcolareilprodottomatricetridiagonalevettoreAx=y conARn,n tridiagonale e x,yRn. Utilizzare esclusivamente vettori. [Tridvet] 3.SvilupparelalgoritmopercalcolareilprodottomatricematriceAB=CconA,B,CRn,n. Ottimizzare quindi lalgoritmo nel caso di matrici triangolari superiori. [Matmat] 4.Calcolare i risultati delle seguenti operazioni tra matrici. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO19 1.B A C B A + 111]1

111]1

3 - 2 12 - 2 31 1 4=5 1 2 -1 0 42 3 1 -=111]1

111]1

2 3 1 -1 - 2 73 4 3= C A B A k k 1]1

21=4 - 1 21 2 - 3=B=32

1]111

1]1111121122 AB C B A 111]1

111]1

1 1 -2 - 31 - 1=6 5 44 3 23 2 1=111]1

111]1

8 - 134 - 72 - 4= C 5.Calcolare il determinante delle seguenti matrici. 1]1

5 20 0A111]1

3 0 12 1 03 2 1B111]1

3 0 02 1 03 2 1C [det(A)=0 , det(B)=4 , det(C)=3] 6.Calcolare il rango delle seguenti matrici. 1]1

1 22 1A=111]1

1 1 03 2 11 1 1= B 111]1

1 2 - 1 -2 - 4 21 - 2 1= C [rango(A)=2 , rango(B)=3 , rango(C)=1] 7.Data la matrice diagonale 1]1

baA=00 calcolare la matrice inversa A1. 1 00 111]1

1]1

/b/a= A- 8.Calcolare la distanza tra i vettoriu=(3, 5, 4) ev=(6, 2,1) e tra i vettorix=(5, 3,2, 4, 1) e y=(2, 1, 0, 7, 2). [ 83 =2v u ,47 =2y x ] 9.Calcolare le norme , 1 e 2 del vettore v=(6, 5). [ 6 v,111 v,612 v ] 10.Calcolare le norme e 1 della matrice 111]1

9 8 76 5 43 2 1A[ 24 A,181 A ] E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO20 1.SISTEMI DI EQUAZIONI LINEARI ( Riferimento al testo: Cap. III ) INTRODUZIONE Unsistema di n equazioni lineari in n incognite xi (i = 1, n) pu essere scritto in forma matriciale come: b Ax =... .......... ............. .......... .................3 3 2 2 1 13 3 2 2 1 11 1 3 13 2 12 1 11' + + + + + + + + + + + +n n nn n n ni n in i i in nb x a ....... x a x a x ab x a ....... x a x a x ab x a x a x a x a dove n nR, A lamatricedeicoefficienti(dimensione(n,n)), nR x ilvettoredelleincognite (dimensione (n, 1)), nR bil vettore dei termini noti (dimensione (n, 1)). La soluzione del sistema di equazioni lineari esiste ed unica se e solose la matrice A non singolare, ciolesuecolonne(righe)sonovettorilinearmenteindipendenti.Inquestocasodet(A)0 (ovvero rango(A) = n) e quindi esiste la matrice inversa A1 tale che x=A1b. I metodi di soluzione si dividono in: METODI DIRETTIGauss (decomposizione GA=U e fattorizzazione LU) Choleski si applicano a matrici A dense e di ordine non elevato (102, 103) effettuanounnumerofinitodioperazionisuA e subchetrasformanoilsistemainizialeinun sistema equivalente pi semplice con matrice dei coefficienti triangolare occorre memorizzare tutti gli elementi della matrice A sono affetti soltanto da errori di round-off METODI ITERATIVIJacobi Gauss-Seidel SOR si applicano a matrici A sparse e di ordine elevato (104, 106) operanoesclusivamenteconlamatriceA iniziale e generano una successione (infinita) di vettori convergente alla soluzione x sufficiente memorizzare soltanto gli elementi non nulli della matrice A sono affetti sia da errori di round-off sia da errori di troncamento analitico (discretizzazione) Nel seguito si affronteranno per primi i metodi di soluzione diretti e successivamente quelli iterativi. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO21 1.SISTEMI TRIANGOLARI Se il sistema di equazioni lineari ha matrice dei coefficienti triangolare superiore, a partire dallultima equazione si risolve il sistema a ritroso con tecnica diback-sostitution: dallultima equazione si ricava il valore dellincognita xn e, per sostituzione progressiva nelle equazioni via via precedenti, si ricavano tutte le altre xi. Le equazioni quindi si disaccoppiano via via e le soluzionixi si ricavano mediante la formula ricorsiva riportata nel seguito. Dato il seguente sistema lineare con matrice dei coefficienti triangolare superiore: ' + + + + + + +n n nnn nn nb x a

b x a ...... x a x ab x a ...... x a x a x a.... ...... .......... ..........2 2 3 23 2 221 1 3 13 2 12 1 11con aii 0, cio A non singolare la sua soluzione si ricava con la formula ricorsiva di back-sostitution (crf. Algoritmo: Backsost): ...n k= na) x a (bxabxkknk jj kj kknnnn' + 1 , , 2 , 1 con1 Esempi 1) 2 1 1 2 00 1 1529414291029101 2 3 42 3 43 44x x x xx x xx xx + + + '2

da cui si ricavano: rrr2: 2r4291029101352941410 1 1 11 2 1 1 2 0 14 43 32 21 1:::x xx xx xx x + + + In modo analogo se il sistema di equazioni lineari hamatrice dei coefficienti triangolare inferiore, a partiredallaprimaequazionesirisolveilsistemacontecnicadiforward-sostitution:dallaprima equazionesiricavailvaloredellincognitax1e,persostituzioneprogressivanelleequazioniviavia successive, si ricavano tutte le altre xi. Il costo per risolvere un sistema triangolare (inferiore o superiore) pari a circa n2/2 operazioni. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO22 1.METODO DI GAUSS Lideadibaseconsisteneltrasformare,conunnumerofinitodicombinazionilineariedeventuali permutazioni di equazioni (cio scambi di righe), il sistema iniziale denso ' + + + + + + + + + + + +n n nn n n ni n in i i in nb x a ....... x a x a x ab x a ....... x a x a x ab x a x a x a x a3 3 2 2 1 13 3 2 2 1 11 1 3 13 2 12 1 11... .......... ............. .......... .................: = b Ax inunsistemaequivalente ~= b Ux conmatricedeicoefficientiUtriangolaresuperiore,chesirisolve mediante tecnica di back-sostitution. Procedura 1) Si effettua leliminazione delle variabili in (n1) passi. Al generico passo k-esimo (k = 1, ... , n1): a)si individua la riga pivot (cio la riga il cui elemento sulla diagonale principale diA aii 0) e la colonna da annullare (operazione di pivoting) b)sicalcolailmoltiplicatoremikcomen...i=k /a a m(k)kk(k)ik ik, , 1 + ,daapplicareagli elementi della colonna da annullare che seguono quello in esame c)sitrasformanoglielementiaij(k) ebj(k)dellamatrice dei coefficienti e del vettore dei termini noti (con i e j > k, cio appartenenti alle righe sottostanti quella in esame) come '+ + + + ++nj=k b m b bni=k a m a a(k)k ik(k)i) (ki(k)kjik(k)ij) (kij, ... , 1, ... , 111Lelemento ) (kkka il pivot. Il costo delleliminazione pari a n3/3 operazioni. 2)Si risolve il sistema triangolare superiore equivalente mediante back-sostitution. Ilcostodisoluzione(n2/2operazionicirca)moltoinferiorealcostodelleliminazioneequindiil costo complessivo del metodo di Gauss circa pari an3/3 operazioni (conta soltanto il processo di eliminazione delle variabili). OSSERVAZIONI SUL PIVOTING 1)La scelta dellelemento pivot viene solitamente effettuata con tecnica di pivoting parziale. Siscegliecomepivotlelementodimodulomassimo della colonna da annullare in modo da avere tutti i moltiplicatori minori di 1 in modulo. Si sceglie cior pari al pi piccolo intero maggiore o uguale dik e tale che ) ( ) (maxkikn i kkrka a e si scambia la riga k-esima con la riga r-esima. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO23 1.2)Se ilpivot ) (kkka nullo, il metodo di Gauss si blocca perch impossibile calcolare i moltiplicatori mik(divisioneper0).Sidevequindisceglierelelementopivotnonnulloscambiandodipostodue equazioni,p.e.lak-esimaconunadellesuccessive.Questopossibileperchilsistemanon singolare e quindi se akkk ( )=0 necessariamente qualche altro elemento della colonna k-esima 0. 3)Si migliora la stabilit del metodo permutando le righe anche se lelemento pivot piccolo in modulo rispettoaglialtrielementi(unpivotpiccolopotrebbederivaredalladifferenzadiduenumeriquasi uguali e quindi essere affetto da cancellazione numerica). 4)Leliminazionedellevariabilinonnecessitadipermutazionidiequazioninelcasodimatricia diagonaledominanteosimmetrichedefinitepositive,purchglielementipivot ) (kkkasiano tutti non nulli. DECOMPOSIZIONE GA=U Il metodo di Gauss pu essere interpretato come una successione finita di trasformazioni diA e b, cio come moltiplicazione di A e b per un numero finito di matrici opportune. Questo tipo di interpretazione consente di riformulare il metodo in due parti distinte: laprima(crf.Algoritmo:Factor)determinaunamatricenonsingolareGtalecheGA=Uuna matrice triangolare superiore, la seconda (crf. Algoritmo: Solve) utilizza la matrice G e risolve il sistema Ax=b. Il sistema lineare iniziale Ax=b pu essere riscritto come: GAx=Gb Ux=~bcon: 1 1 2 2 1 n 1 nP M P M ...... P M G Scambidirighe-LoscambiodidueequazionidelsistemaAx=b, p.e. la rigai-esima con la rigaj-esima,equivaleamoltiplicare(dasinistra)entrambiimembridelsistemaperlamatricePcheuna matrice identit con le righe i e j scambiate. Combinazionilineari-Lasostituzionenelsistemadellarigai-esimaconlamedesimapilarigaj-esimamoltiplicatapermij,equivaleamoltiplicare(dasinistra)ilsistemaperlamatriceMcheha diagonale principale unitaria e in posizione (i, j) il moltiplicatore mij. FATTORIZZAZIONE PA=LU La decomposizione GA=U pu essere riformulata ipotizzando di effettuare prima tutti gli scambi di righe e successivamente tutte le combinazioni lineari. Le matrici che contengono i moltiplicatorimij saranno ovviamente ordinate in modo diverso perch gli scambi di righe avranno agito anche su di esse; quindi invece di avere: ( ) U A P M P M ...... P M GA1 1 2 2 1 n 1 n E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO24 1. ora si ha: ( )( ) U A P P ...... P M M ...... M GA1 2 1 n 1 2 1 n U PA M GA dove M la matrice triangolare inferiore, a diagonale unitaria, dei moltiplicatori con le righe riordinate dal vettore pivot. In particolare si pu scrivere: U M PA1 cioPA=LU dove 1M L unamatricetriangolareinferiorecondiagonale unitaria ed elementi fuori diagonale pari ai moltiplicatori mij cambiati di segno e riordinati come detto dal vettore pivot. Lalgoritmo Factor fornisce direttamente le matrici L e U se si effettuano le seguenti modifiche ai passi: 6.akj ai jj=1, ..., n (scambio di righe anche nella parte inferiore di A che via via contiene i moltiplicatori) 9.aik= mik=aaikkk (si memorizza il moltiplicatore cambiato di segno) 10.aij=aijaikakj ,j=k+1, ..., n (cambia il segno della combinazione lineare) Nota la fattorizzazione PA=LU, per ricavare la soluzione del sistema iniziale Ax=b sufficiente risolvere in cascata i sistemi triangolari: ' b Pb Lyy Ux InfattidaAx=bsiricavaPAx=Pb(moltiplicandoamboimembriperP),cioLUx=Pb(perch PA=LU) da cuiy Ux e b Ly . UTILIT DELLA FATTORIZZAZIONE La fattorizzazione PA=LU, come anche la decomposizione GA=U, possono essere utilizzate per: 1.Calcolare la matrice inversa A1 PA=LU (PA)1 = (LU) 1 = U1L1 A1P1 = U1L1 Dato che P1 = P (ci sono soltanto 0 e 1) si ha A1 = (U1L1)P A1 = BP E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO25 1.FareilprodottoPAsignificaeseguiresuAgliscambidirigheriportatinelvettorepivot; analogamente il prodotto BP corrisponde a eseguire gli stessi scambi sulle colonne. 2.Risolvere p sistemi lineari multipli Ax bAx bAx b1 12 2p p'....... .... lacuisoluzionemediantedecomposizioneGA=Uesoluzionedeisistemitriangolaricosterebbe ( ) 2 32 3n n p +operazioni. La fattorizzazionePA=LU (costo33n ) viene calcolata una sola volta perch le matriciL, U eP sono le stesse per tutti i p sistemi lineari. Si ottengono quindi i seguenti sistemi triangolari: Ly bUx yLy bUx y1 1*1 1p p*p p'''costo costo nn2222...........il cui costo di soluzione ( )2 2 22 2 pn n n p + Ilcostocomplessivo(fattorizzazione+soluzionedeisistemitriangolari)quindiparian pn3 23 +ed decisamente inferiore al costo da sostenere con la decomposizione GA=U. METODO DI CHOLESKI SelamatriceAadiagonaledominanteosimmetricadefinitapositiva,ilmetododiGaussprocede senzanecessitdieffettuarescambidirighe.LaconseguentefattorizzazioneA=LUpuessere interpretata come prodotto di particolari matrici triangolari, una inferiore e laltra superiore, A=L1L1T. La fattorizzazione A=LU pu infatti essere riscritta come: A=LDU1 dove D=diag(U), L e U1 sono triangolari (superiore e inferiore) con diagonale unitaria. Poich la matrice A simmetrica allora U1=LT e quindi A=LDLT. LamatriceAanchedefinitapositivaequindiglielementidiagonale(D)iisonotuttipositivi;sipu introdurreDscrivendo A=LD1/2 D1/2LT Poich la matrice D diagonale allora DT=D e quindi A=(LD1/2) ((D1/2)TLT)=L1L1T con L1=LD1/2. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO26 1.111111111]1

nn n n nii j jl l l ll l ll l ll ll3 2 12 133 32 3122 21110 00 00 00 0 0=1L Gli elementi lij si calcolano con la seguente formula ricorsiva: '

,_

1 , ... , 1 22 11121111 11i j= ,...,ni= l a lll l al a =l/i-k=ik ii iijjjk=jk ik ijij Il costo pari a63n , cio la met della fattorizzazione LU vera e propria. CONDIZIONAMENTO DEI SISTEMI LINEARI LasensitivitdellasoluzionediunsistemalineareallevariazionideicoefficientidellamatriceA e dei termini noti del vettore b viene esaminata attraverso lo studio del condizionamento. Sia nel caso di perturbazione del solo vettoredei termini noti b sia nel caso di perturbazione del vettore dei termini noti b e della matrice dei coefficienti A, si definisce indice di condizionamento del problema il numero: 1) ( A A A K cherappresentailfattorediamplificazionedelleperturbazionirelative(ciodeglierrorirelativi) introdotte nel vettore b e nella matrice A. Se K(A)>>1 il sistema lineare Ax=b mal condizionato; K(A) tanto maggiore quanto pi la matrice dei coefficienti A tende a essere singolare. Tipici esempi di mal condizionamento sono la matrice di Hilbert Hn e la matrice di Vandermonde Vn: E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO27 1.H Vn n + + +

1]11111

1]11111 1 1 2 11 2 1 3 1 11 3 1 4 1 21 1 1 1 2 11 1 11 21222 21121 1/ // / / ( )/ / / ( )/ / ( ) / ( )nnnn n nx x xx x xx x xnnn nnn METODI DI SOLUZIONE ITERATIVI Questi metodi si applicano a matrici dei coefficienti A sparse e di ordine elevato (104, 106). OperanoesclusivamenteconlamatriceinizialeAegeneranounasuccessioneinfinitadivettoriche converge alla soluzione (vettore delle incognite x). Procedura La matrice dei coefficienti A pu essere pensata come somma di due matrici reali, quadrate di ordine n: A=C+D In questo modo si ha: Ax=b (C+D)x=b Cx +Dx=b Dx=bCx Dx=d con b, x Rn e d funzione di x. Il procedimento iterativo costituito dai seguenti punti: 1.si assume quale vettore soluzione iniziale un vettore qualsiasi x(0) (p.e. x(0)=0 vettore nullo) 2.si costruisce il vettore d(0)=bCx(0) 3.si risolve il sistema lineare Dx(1)=d(0) ricavando il vettore x(1) 4.si costruisce il vettore d(1)=bCx(1) 5.si risolve il sistema lineare Dx(2)=d(1) ricavando il vettore x(2) 6.. e cos via 7.siconsideraraggiuntalaconvergenzaquandoladifferenzatraduevettorisoluzionesuccessivi minore di un prefissato valore di soglia Assunto il vettore iniziale x(0), il procedimento iterativo viene quindi espresso dalle formule: ') ( ) 1 + () ( ) ( iterazioni di numero , ... 2, , 1 , 0 = con=k kk kkd x DCx b d Per calcolare il vettore soluzione ) 1 + (kxsi deve quindi risolvere il sistema lineare ) ( ) 1 + ( k kd x D . La matrice dei coefficientiD deve perci essere non singolare ed conveniente se ha struttura idonea da consentire calcoli rapidi e agevoli. Einoltrenecessariochelasuccessionedeivettorisoluzioneconvergaallasoluzionex:lamatriceD deve essere scelte in modo da garantire tale convergenza (la matrice C=AD viene di conseguenza). E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO28 1.CONVERGENZA Unmetodoiterativoconvergenteselerroreassolutoe x x( ) ( ) k k tende a zero allaumentare del numero k di iterazioni. LasceltadellamatriceD(equindidella matriceC) che garantisce la convergenza viene fatta proprio ragionando sullerrore assoluto. Dalla scrittura ) ( ) ( ) 1 + ( k k kx C b d x D si ricava( )) ( 1 ) ( ) 1 + ( k k kx C b D d x e quindi: ( ) ( ) ( )) 0 ( 1 ) 0 ( ) 1 () ( ) ( 1 ) ( 1 ) ( 1 1 ) 1 ( ) 1 + (... e B e BB BBeBe Ce D x x C D x C b D x C b D x x e+ + k k kk k k k k k dovee(0) lerroreassolutoinizialeeB=D1Clamatricediiterazioneche,aogniiterazione, moltiplica lerrore. La matrice di iterazioneB=D1C rimane costante nel corso di tutta la soluzione perch dipende dalla matrice dei coefficienti iniziale A che a sua volta non viene modificata dai metodi iterativi. Il processo iterativo di soluzione certamente convergente se: ||ID1A||j. La matrice input A stataottenutadallalgoritmoFactor.Ilvettorepivotcontienegliscambidiriga effettuati da Factor. Al termine il vettore b contiene la soluzione x. Parametri.Input: n, A, pivot, b Output: b 1.Ciclo 1: k=1, ..., n1 2.j=pivot(k) 3.se jk, bj bk 4.Ciclo 2: i=k+1, ..., n 5.bi=bi+aik bk 6.Fine Ciclo 2 7.Fine Ciclo 1 8.bn=bn / ann 9.Ciclo 3: i=n1, ..., 1 E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO32 1.10.bi=b a b ai i l ll ini i

_,

+1/11.Fine Ciclo 3 12. Exit End Algoritmo: Gseidel (n, A, b, toll, kmax, x, ier) Commento.Lalgoritmo,utilizzandoilprocessoiterativodiGauss-Seidel,migliora lapprossimazione inizialex=x(0) della soluzione del sistema non singolare di ordine n Ax=b.Lesuccessiveapprossimazionivengonomemorizzatenelvettorex.Sela precisione richiesta toll raggiunta con un numero di iterazioni kmax, si pone ier=0, altrimenti ier=1. Parametri.Input: A, , b, toll, kmax, x Output: x, ier 1.Ciclo 1: k=1, ...,kmax 2.y=x1 3.x1=1121 1a x a bnjj j

,_

4.ermax = yx1 5.Ciclo 2: i=2, ..., n 6.y=xi 7.xi=iini jj ijijj j ia x a x a b

,_

+ 1111 8.er= y-xi 9.se ermax < er , ermax=er 10.Fine Ciclo 2 11.se ermax < toll x, ier=0; Exit 12.Fine Ciclo 1 13. ier=1 14. Exit End ESERCIZI SVOLTI 1.Risolvereilsistemalineareproposto mediante la decomposizioneGA=U e con tecnica di pivoting parziale. 1111]1

1111]1

1111]1

40331 0 1 21 1 2 01 0 1 10 1 0 24321xxxx IndicareinoltrelematriciUeM(matricetriangolareematricedeimoltiplicatori)eivettori ~b(termine noto trasformato) e pivot (contenente gli scambi di righe) e calcolare det(A). E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO33 1. passo 1: a11 0edipigrandedeglialtrielementidellaprimacolonnaquindila1^riganonsiscambia (rimane al suo posto) e pivot(1)=1 2 0 1 01 1 0 10 2 1 12 1 0 13304

1]11111 si azzera la colonna1: 2^ 2^121124^ 4^ 1 1 12141 '^^

mm ottenendo 2 0 1 00 1 1 2 10 2 1 10 1 1 133 201

1]11111/ / passo 2: a22 0mapipiccolo(inmodulo)dia32:sieffettualoscambio2^riga3^rigaequindi pivot(2)=3 2 0 1 00 2 1 10 1 1 2 10 1 1 1303 21

1]11111/ / si azzera la colonna2: 3 3122124 4122123242^ ^ ^m^ ^ ^ m + + ' ottenendo 2 0 1 00 2 1 10 0 0 3 20 0 1 2 3 2303 21

1]11111// // passo 3: a33 = 0: si effettua lo scambio 3^ riga 4^ riga e quindi pivot(3)=4 2 0 1 00 2 1 10 0 1 2 3 20 0 0 3 23013 2

1]11111/ // / la colonna3 gi azzerata: m43 = 0 Alla fine del processo di eliminazione delle variabili si giunti a: 1111]1

1111]1

1111]1

2 / 31032 / 3 0 0 02 / 3 2 / 1 0 01 1 2 00 1 0 2

~=4321xxxxb x Upivot=(1,3,4)T E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO34 1.1111]1

1 0 2 1 10 1 2 1 00 0 1 2 10 0 0 1///=M det(A)= 32321) 2 ( 2 ) 1 (2 ,_

,_

Risolvendo Ux=~bcon back-sostitution: '+ =k= nux u bx n= /u b xkkjnk jkj kknn n n1 , 2 , 3 1 , ,... 141 ' + + + + + +2323123210 1 1 23 0 1 0 244 34 3 24 3 2 1xx xx x xx x x x x4=1, x3=1, x2=1, x1=1 x=(1, 1, 1, 1)T 2.Ricavare la fattorizzazione LU della matrice A proposta nellesercizio, noto il vettore Pivot. A

1]11112 0 1 01 1 0 10 2 1 12 1 0 1 Pivot=13412 33 4

_,

r kr r kr r knon scambia (passo= 1) (passo= 2) (passo= 3) Rispetto alla decomposizione GA=U sufficiente applicare le seguenti modifiche: -cambiare di segno i moltiplicatori mij; -effettuare gli scambi memorizzati nel pivot anche sulla parte inferiore di A. passo 1: Pivot(1)=1 quindi la riga r1 non si scambia. Si azzera la colonna 1 e si memorizza m21=1/2em41=1 nelle rispettive posizioni di A 1111]1

1 1 1 11 1 2 01 2 1 1 2 10 1 0 2/ / passo 2: Pivot(2)=3 quindi si scambia 2^ riga 3^ riga E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO35 1.2 0 1 00 2 1 11 2 1 1 2 11 1 1 1

1]1111/ / Si azzera la colonna 2 e si memorizza m32=1/2em42=1/2 nelle rispettive posizioni di A 2 0 1 00 2 1 11 2 1 2 0 3 21 1 2 1 2 3 2

1]1111/ / // / / passo 3: Pivot(3)=4 quindi si scambia 3^ riga 4^ riga 2 0 1 00 2 1 11 1 2 1 2 3 21 2 1 2 0 3 2

1]1111/ / // / / La colonna 3 gi azzerata; si memorizza m43=0 nella rispettiva posizione di A Alla fine del processo di eliminazione delle variabili si giunti a: 1111]1

1111]1

2 3 0 0 02 3 2 1 0 01 1 2 00 1 0 21 0 2 1 2 10 1 2 1 10 0 1 00 0 0 1// /=/ //=U L Si pu verificare che PA=LU. 3.Sviluppare lalgoritmo di soluzione del sistema lineareAx=b, conA tridiagonale (si pu eliminare il pivoting) di ordine n, utilizzando esclusivamente vettori. Algoritmo: Tridigsolu (n, c, d, f, b, l, u, x) Commento.CalcolalasoluzionediAx=bconAtridiagonale.Utilizzasoltantoivettorid=diag., c=codiag.inf.ef=codiag.sup.diA.RicavalafattorizzazioneLUinterminidivettori l=codiag. inf. di L, u=diag. di U e f=codiag. sup. di U (coincide con la codiag. sup. di A). La soluzione effettuata come cascata di sistemi triangolari. Parametri.Input: n, c, d, f, b Output: x 1.u1=d1 2.Ciclo 1: i = 2,..., n (fattorizzazione LU di A tridiagonale) 3.li= ci/ui1 4.ui= dilifi1 E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO36 1.5.Fine Ciclo 1 6.y1=b1 7.Ciclo 2: i = 2,..., n (soluzione sistema triangolare Ly= b) 8.yi=biliyi1 9.Fine Ciclo 2 10.xn=yn/un 11.Ciclo 3: i = n1,..., 1 (soluzione sistema triangolare Ux = y) 12. xi=yifixi+1/ui 13.Fine Ciclo 3 14.Exit End ESERCIZI PROPOSTI 1.Sviluppare lalgoritmo per la soluzione di un sistema diagonale Dx=b di ordine n. [Diagsolu] 2.Sviluppare lalgoritmo per la soluzione di un sistema triangolare inferiore Lx=b di ordine n. [Forwsost] 3.Risolvere i seguenti sistemi lineari triangolari superiori con tecnica di Forward-sostitution. 1111]1

1111]1

1111]1

40331 0 1 20 1 2 20 0 1 10 0 0 24321xxxx [x=(1.5, 1.5, 6, 0.5)T] 1111]1

1111]1

1111]1

66331 2 1 20 3 5 40 0 2 30 0 0 14321xxxx [x=(3, 3, 3, 3)T] 4.Ricavare la fattorizzazione LU della seguente matrice A, noto il relativo vettore Pivot. A

1]11114 7 9 25 4 6 13 3 4 15 7 2 8 Pivot

_,

224 1111]1

1111]1

1111]1

139 1 0 0 019 115 19 139 0 05 6 5 21 5 19 01 6 4 51 139 5 19 3 5 30 1 19 15 10 0 1 5 40 0 0 1// // / / U = / / ///L = E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO37 1. 5.Ricavare la fattorizzazione LU della seguente matrice A, noto il relativo vettore Pivot. A

1]11111 2 8 22 4 3 25 1 2 15 3 4 2 Pivot

_,

323 1111]1

1111]1

1111]1

117 1 0 0 01 2 13 0 05 8 5 11 5 18 01 2 1 51 117 14 9 5 10 1 2 1 5 10 0 1 5 20 0 0 1/// / / U = / // //L = 6.Ricavare la fattorizzazione LU della seguente matrice A, noto il relativo vettore Pivot. A

1]11116 5 3 82 2 10 37 2 8 21 4 9 7 Pivot

_,

343 1111]1

1111]1

1111]1

281 1 0 0 026 9 26 281 0 07 47 7 55 7 26 02 8 2 71 281 122 13 5 7 20 1 26 23 7 60 0 1 7 10 0 0 1// // / / U = / / // //L = E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 38 AUTOVALORI DI MATRICI Riferimento al testo: Cap. IV INTRODUZIONE Moltiproblemidellingegneriasonodescrittidasistemidiequazionilineariconfunzionichedipendono, oltre che dalle incognite, anche da un parametro .Ilsistemaammettesoluzionediversadallasoluzionenullaperparticolarivaloriidettiautovaloriele corrispondenti soluzioni xi sono dette autovettori ffnn n1 1100( , ..., ; )...........................( , ..., ; ) ' in forma matriciale si haAx=x cio(A-I) x=0 RautovaloredellamatriceARn,nseesoloselamatrice[A-I]singolare,cioseesolose det(A-I)=| A-I |= 0 (equazione caratteristica). Se ARn,n simmetrica allora R e i corrispondenti xRn sono un sistema di vettori ortogonali. Se A simmetrica definita positiva allora i suoi R e sono tutti positivi. Se autovalore di A allora 1/ autovalore di A-1. A e AT hanno gli stessi autovalori. Lequazione caratteristica det(A- I)=0 un polinomio di grado n det(A- I)=0cio0 ..... ) 1 (111 + + + + n nn n n Lasoluzionedirettadellequazionecaratteristicanonvieneeffettuatainquantositrattaspessodi problema mal condizionato. Gli autovaloriisicalcolanoutilizzandodiversimetodinumericiesuccessivamentesirisolveciascun sistema lineare omogeneo associato (A-iI)xi=0 per calcolare i corrispondenti autovettori xi. METODO DELLE POTENZE Serve per calcolare lautovalore di modulo massimo o quello di modulo minimo della matriceARn,n. Autovalore di modulo massimo Data una matriceA si vuole calcolare lautovalore di modulo massimo (e il corrispondente autovettore x). Vengono fatte due ipotesi. Hp1: esiste un solo 1Rn di modulo massimo( ) 1 2 3> ....n, cio i 11 < . Hp2: A diagonalizzabile (cioX-1AX= con=diag(1, ...,n)) e quindi tutti gli autovettorixi sono linearmente indipendenti tra loro. Grazie allHp2, il generico vettore v0 pu essere scritto come combinazione lineare degli autovettori. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 39 v x x x v x0 1 2 n i + + + 1 2..... ,nnR Nel metododelle potenze si sceglie un vettore inizialev0 qualsiasi (p.e.v0 =(1, 1, 1,..., 1)T ) e quindi si genera la seguente successione di vettori: v A v Ax Ax x x x x xv A v A Av x x x x1 0 1 n 1 n 1 2 n2 1 0 1 n 1 2 + + + + +

_, + +

_,

_,

+ + +

_, + + 1 1 1 1 121211 12 212121221... ... ...... ...n n nnnn nn

_,

_,

+ + +

_, + +

_,

_21 1 1 12121 nmn nm mmnmnxv A v A Av x x x x xnm m 1 0 1 n 1 2 n.................................................... ... ...,

Al passom il vettorevm, per effetto della moltiplicazione progressiva per la matriceA, tende a disporsi parallelamente allautovettore x1. Inoltre, grazie all Hp1 ( i 11 < ) si ha: lim limmimmm

_, 1 1101quindi v xm 1. Lautovaloredimodulomassimocercatosiottienecomerapportotralegenerichecomponentik0 dei vettori vm e vm+1: ( )( )( )( )vvxxm 1m11+

_, + + kkkk0011 01 0100 terminitermini Algoritmo.Pow(n, A, toll, mmax, , y) Commento.Determina lautovalore di modulo massimo della matrice A e lautovettore corrispondente. Se la precisione relativa richiesta toll viene raggiunta con un numero di iterazioni mmax la variabile ier assume il valore 0, altrimenti ier=1. Parametri.Input: n, A, toll, mmax Output: , y 1.y0=(1, 1, ..., 1)T 2.Ciclo 1: m=0, ..., mmax 3. wm+1=A ym 4. 11 00( )( )( )mkk++wym 1m (approssimazione dellautovalore di modulo massimo cercato) 5. ywwm 1m 1m 1+++(normalizzazione dellautovettore per evitare overflow o underflow.) 6.er=( ) ( ) 1 11 m m+ 7se ertoll( )11 m+| oppure ertoll allora poni ier=0 e vai al punto 10. 8.Fine Ciclo 1 E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 40 9.ier=1 10. y=ym+1 11. ( )11 m+ 12. Exit End Osservazioni Se1=0e|2|>|3|,inteoriailmetododellepotenzedovrebbeconvergeresullautovalore2. Per, causa degli inevitabili errori di round-off , dopo pochi passi si ha 10 e il metodo ricade su 1. Se 1 2 laconvergenzapuessereeccessivamente lenta allora si utilizza il metodo delle potenze per avere una stima iniziale p di 1 e poi si raffina con il metodo delle potenze inverse. Autovalore di modulo minimo Data una matriceA si vuole calcolare lautovalore di modulo minimo (e il corrispondente autovettore x). Vengono nuovamente fatte 2 ipotesi. Hp1: esiste un solo 1Rn di modulo minimo( ) 1 2 3< ....n. Hp2: Tutti gli autovettori xi sono linearmente indipendenti tra loro. Poich se autovalore diA, lautovalore diA-1 pari a 1/ (crf. Pag. 7-1), per calcolare lautovalore dimodulominimodellamatriceAsufficiente calcolarelautovaloredimodulomassimo dellamatrice inversa A-1 e quindi calcolarne il reciproco. Lo schema di calcolo il seguente A x x Bx xB A-1 =metodo delle potenze per calcolareautovalore di modulo diquindi si calcola =11111 maxmin OvviamenteanzichcalcolareesplicitamentelamatriceinversaA-1(coston3)convieneeffettuarela fattorizzazione LU (costo n3/3) e risolvere i due sistemi triangolari risultanti (costo n2/2 ciascuno). Il metodo delle potenze viene portatoavanti finch 111 11 ( ) ( ) ( ) m m m + + ; a questo punto si assume 11 ( ) m+ quale migliore approssimazione cercata e se ne calcola il reciproco. 1.y0=(1, 1, ..., 1)T 2.Fattorizzazione PA=LU (algoritmo Factor) 3.Ciclo 1: m=0, ..., mmax 4.LUwm+1=Pym '++ +Lz PyUw z11 1m mm m E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 41 5.1110 0( )( ) ( )mm k m kw y++(approssimazione di 1 autovalore di modulo massimo di A-1) 6.y w wm m m + + +1 1 1 (normalizzazione autovettore) 7.Fine Ciclo 1 8.n 11 (autovalore di modulo minimo di A cercato) METODO DELLE POTENZE INVERSE E una generalizzazione del metodo delle potenze e serve per calcolare un particolare autovalore di cui si conosca una stima p. Utilizzando la stima p, il sistema Ax=x pu essere riscritto come (ApI)x=Axpx=xpx=( p)x cio( p) autovalore della matrice (ApI) e lautovettore corrispondente sempre x. Se ( p) autovalore della (ApI) allora 1/(p) autovalore della matrice inversa (ApI)-1; inoltre se punabuonastimadi(pvicinoa)allora1/( p)lautovaloredimodulomassimo della matrice (ApI)-1 che pu essere calcolato con il metodo delle potenze. Lo schema di calcolo il seguente ( )max ( )A I x x Bx xB A I pp ppp-= con metodo delle potenze per calcolareautovalore di modulo diquindi si calcola = + 11111 11 E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 42 Algoritmo.Invpow (n, A, toll, mmax, p, x, ier) Commento.UtilizzailmetododellepotenzeinverseperdeterminarelautovaloredellamatriceA, di ordine n, pi vicino al numerop stimato e il corrispondente autovettore x. Se la precisione relativarichiestatollvieneraggiuntaconunnumerodiiterazionimmaxsiponeier=0, altrimenti ier=2. Se la matrice (ApI) singolare, allora ier=1. Parametri.Input: n, A, toll, mmax, p Output: p, x, ier 1.( ) ( ) , , ...., A Aii iip i n 12.richiama lalgoritmo Factor (n, A, pivot, det, ier) 3.se ier=1 Exit 4.( ) yo 11 1 , ,...T 5.p(0)=p 6.Ciclo 1: m=0, ..., mmax 7.Um+1=Gym m+1 8.a= m 1 +; sia k0 la posizione della prima componente di modulo massimo di m+1 9. ym 1m 1++ a (normalizzazione del vettore m+1) 10.p(m+1)=p+( )( ) ymmkkoo+1 11.er=( ) ( ) pmpm+ 1 12.se ertoll( )mp+1| oppure ertoll poni ier=0 e vai al punto 15 13.Fine Ciclo 1 14.ier=2 15.x=ym+1 16.p=p(m+1) 17.Exit End TRASFORMAZIONI DI SIMILITUDINE Sono necessarie alcune definizioni preliminari. 1. Due vettori sono ortogonali se e solo se il loro prodotto scalare nullo (xTy= x yiini=1 =0). 2. Un sistema di vettori ortonormale se i vettori a1, ..., an sono ortogonali due a due. 3. Unamatriceortogonaleseesoloselesuerighe(colonne)formanounsistemadivettori ortonormale. 4. Se A ortogonale si ha ATA=AAT=I, A-1=AT; se A ortogonale e simmetrica allora A-1=AT=A. UnriflettoreelementareUunamatricediordinennonsingolareeortogonale(U-1=UT).Una trasformazionedisimilitudineunatrasformazionecheassociaaunamatrice A la matrice simileU-1AU=UTAU. Esseresimili significa che la matrice iniziale e quella trasformata hanno gli stessi autovalori e autovettori semplicemente trasformati. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 43 METODO QR Serve per calcolare tutti gli autovalori di una matrice A. Poich il calcolo di tutti gli autovalori risulta estremamente costoso se applicato a matrici dense, il metodo trasformalamatriceinizialeAinunamatricesimile,diformapisemplice(tridiagonaleselaA simmetrica), di cui calcola autovalori e autovettori. DatalamatriceAsicostruisceunamatriceRtriangolaresuperiorecomeprodottodi(n1) riflettori elementari Ui Un-1Un-2...U1A=R che si scrive anche come QTA=R dove QT=Q-1, ortogonale in quanto prodotto di matrici ortogonali. Quindi si pu scrivere A=QR dove la matriceR, simile alla matrice inizialeA, triangolare superiore e ha gli autovalori posizionati sulla sua diagonale principale. GliautovaloridellamatriceinizialeAsonoglielementisulladiagonaleprincipaledellamatrice simile A ottenuta dalla seguente successione di trasformazioni 1.A1=A 2.Ciclo 1: i=1, ..., imax 3.Ai=QiRi (fattorizzazione QR) 4.Ai+1=RiQi 5.Fine Ciclo 1 6.=diag(Aimax) IlprocessoiterativovienefermatoquandoglielementifuoridiagonaledellaAi,chedovrebbero convergereazerodopoinfinitipassi,sonosufficientementepiccoli.Peraccelerarelaconvergenzasi utilizzano parametri di accelerazione ti opportunamente scelti e modificare il metodo nei passi 3.(AitiI)=QiRi (fattorizzazione QR) 4.Ai+1=RiQi+tiI Il costo della fattorizzazione QR doppio rispetto alla fattorizzazione LU (2n3/3 contro n3/3 operazioni ), per applicabile anche se la matrice singolare (al contrario della fattorizzazione LU) o rettangolare (si utilizza nel metodo dei minimi quadrati). E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 44 APPROSSIMAZIONE DI DATI E DI FUNZIONI Riferimento al testo: Cap. VINTRODUZIONE Approssimazione di funzioni: laf(x) nota analiticamente, ma risulta difficile o impossibile eseguire operazioni (p.e. integrazione) con strumenti dellanalisi matematica; si approssima la f(x) con una fn(x) pi semplice. Approssimazionedidati: laf(x) non nota analiticamente, ma si ha a disposizione una raccolta di dati di cui si conoscono i valori{ } yi (ordinate) corrispondenti ai nodi{ } xi (ascisse). Sicostruisceunmodellomatematicofn(x) (cio una funzione) che approssima in modo attendibile il valore della y in punti diversi dai nodi. Per effettuare lapprossimazione di dati e funzioni necessario: 1.individuare la classe delle funzioni approssimanti Fn={fn(x)} in base alle caratteristiche del fenomeno 2.scegliere il particolare elemento fn(x) applicando un criterio di scelta. CLASSI FNDI FUNZIONI APPROSSIMANTI Polinomi algebrici di grado n Pn={ } f x a a x a x a xn nn( ) ... + + + +0 1 22. Si utilizzano per funzioni continue su intervalli chiusi e limitati. Richiedono la determinazione di (n+1) parametri a0, ..., an. Polinomi trigonometrici di grado n e pulsazione Tn= f x a a k x b k xn kknk( ) ( cos( ) sin( ) + +'; 01 Si utilizzano per funzioni periodiche. FunzionirazionaliRn,d=PPnd= f xa a x a x a xb b x b x xnnnd( )......+ + + ++ + + +';0 1 220 1 22.Siutilizzanoperfunzioni aperiodiche su intervalli illimitati. Richiedono la determinazione di (n+d+1) parametri (il denominatore viene normalizzato in modo che il coefficiente del termine xd sia unitario). PolinomialiatrattiFn,dassociateasottointervallidellintervallodiinteresse.Ognitrattoun polinomio di ordine basso (primo o secondo). Spline Sn,d di gradod. E un particolare sottoinsieme delle funzioni polinomiali a tratti che garantisce la continuit della funzione e di tutte le sue derivate di ordine (d1) in tutto lintervallo dinteresse. Sommeesponenzialidiordinen En= f x a en kb xknk( ) ';1.Richiedonoladeterminazionedi2n parametri ak, bk. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 45 CRITERI DI SCELTA DELLA FUNZIONE fN(x) FN Il criterio di scelta garantisce lindividuazione della funzione approssimante fn(x) in modo univoco. I criteri pi comunemente utilizzati per lindividuazione della funzione fn(x) sono: 1.Interpolazione di dati - si applica nel caso di dati (punti) privi di errore e serve per selezionare la funzione che passa per i punti in esame 2.Approssimazione (o Smoothing) di dati - si applica nel caso di dati (punti) affetti da errore o dispersi e serve per selezionare la funzione che passa (al meglio possibile) tra i punti in esame. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 46 INTERPOLAZIONE DI DATI Riferimento al testo: Cap. V CRITERIO DI INTERPOLAZIONE POLINOMIALE DI DATI Dati gli (n+1) punti noti {xi, yi}, con i=0, ..., n, si determina il polinomio interpolatore Pn(x) di grado n fn(x)= Pn(x) = a0 + a1x + a2x2 + ... + anxn imponendone il passaggio per i punti, cio imponendo le seguenti (n+1) condizioni sulle ordinate fn(xi)=yiperi=0, ..., nconfn(x)Pn(x). Esiste unteorema (Weierstrass) che garantisce lesistenza e unicit del polinomio interpolatore Pn(x) tale che ||f(x)-Pn(x)||0) per tutti gli x [ ] a b ,intervallo di interesse. Se si impongono le condizioni indicate (passaggio della funzione per i punti) e si ricerca il polinomio come risoluzione del sistema di equazioni lineare in (n+1) incognite a a x a x ya a x a x ynnn n nnn0 1 0 0 00 1+ + + + + + '........... .................. 1110 01 10101x xx xx xaaayyynnn nnn n......... ... .... .... ........... ...

_,

_,

_,

= VnTa=y si ricava una matrice dei coefficienti che la Vn (matrice di Vandermonde). Se gli (n+1) nodixi sono distinti, cioxixj, det(Vn) 0 la Vn non singolare e quindi esiste un unico polinomio interpolante Pn(x). Il sistema per mal condizionato e quindi non possibile risolverlo direttamente per determinare i coefficienti del polinomio. Ladeterminazionedellunicopolinomiodiinterpolazionetraipuntidatisieffettuamediantedue tecniche numeriche: il metodo di Lagrange e il metodo di Newton. METODO DI LAGRANGE: POLINOMI FONDAMENTALI DI LAGRANGE Il polinomio di interpolazione Pn(x) si ricava dalla combinazione lineare di polinomi fondamentali Pn(x)= y l xj jjn( )0 Glilj(x) sono detti polinomi fondamentali di Lagrange associati ai nodi, costituiscono una base dello spazio dei polinomi e si rappresentano come E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 47 l xx xx xjii i jnj ii i jn( )( )( ),, 00tali che l xi=ji jj i ij( ) '10 se se cio i polinomi fondamentali di Lagrange sono nulli in tutti i nodi tranne che nel nodo j-esimo. Esempio 1) Costruire il polinomio fondamentale di Lagrange l0(x). l0(x) nullo in tutti i nodi tranne che nel nodo 0, cio l0(x0)=1 e l0(xj)=0 per j=1, ...,n. l0(x) deve annullarsi in x1, x2, ..., xn, quindi necessaria una forma del tipo ( )( ) ( )x x x x x xn 1 2.... che deve contemporaneamente essere unitaria in x0, cio ( )( )( ) ( )( )( ) ( )( )( )l xx x x x x xx x x x x xx xx xnniiniin01 20 1 0 2 0101 ........ Si spiega quindi la forma generale dei polinomi fondamentali di Lagrange lj(x). Il calcolo di ciascuno deglin+1 polinomi fondamentali costa 2n2 operazioni, quindi il costo totale paria(n+1)(2n2)=2n22cherisultamoltomenoonerosodellasoluzionedelsistemadiequazioni lineari (coston3/3). Inoltre il metodo di Lagrange preferibile visto il mal condizionamento del sistema di equazioni per la presenza della matrice di Vandermonde. Esempio 2)Datiiseguentin+1=3punti(0,0),(1,1),(2,0)scrivereilpolinomiodiinterpolazioneP2(x) utilizzando il metodo di Lagrange. Il polinomio si presenta nella forma P2(x) = yx xx xyx xx xyx xx xoiiiiii iii iii iii i( )( )( )( )( )( ),,,,++ 1201210 1210 1220 2220 22 ( )( )( )( )( )( )( ) + + + 0 10 21 0 1 20 1212 22x x x xx x x x Poich una sola ordinata diversa da zero, lapplicazione del metodo particolarmente rapida. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 48 Nelcasogeneralediordinatetuttenonnullesufficienteeffettuarelacombinazionelineare dei polinomi fondamentali lj(x) associati ai j nodi in esame per j=0, ..., n. ESERCIZI SVOLTI 1.Calcolareilpolinomiointerpolatoreperi3punti(0,0),(1,1),(2,3)utilizzandoilmetododi Lagrange. Sono dati n+1=3 punti il polinomio cercato di grado n=2. P20 0200 0210 1210 1220 2220 220 120 120 220 220 11321( )( )( )( )( )( )( )( )( )( )( )(,,,,,,,,,,x yx xx xyx xx xyx xx xx xxx xxx xoii iii iii iii iii iii iii iii iii iii i++ ++ 0 20 20 10 12 221 132 210 21 0 1 230 12 0 2 11213122321 232321212)( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( ) ( )x xx xx x x xx xx x x xx x x xx x x x x x x xx x + + + + + + + METODO DI NEWTON: DIFFERENZE DIVISE Il polinomio di interpolazione Pn(x) viene costruito esplicitamente nella forma Pn(x) =a a x x a x x x x a x x x x x xn n 0 1 0 2 0 1 0 1 1+ + + + ( ) ( ) ( ) ..... ( ) ( ) ...( ) dove i coefficientiai sono dati dalle cosiddettedifferenze divise, rapporti incrementali (numeri) che si definiscono sugli n+1 nodi( , , , ..., ) x x x xn 0 1 2 [ ] [ ] [ ]a f x xf x if x x f x xx xii i i ii >' 001 0 1000, ...,( ), ..., , ..., = . E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 49 Nel caso di 2 nodi( , ) x x0 1, si definisce la differenza divisa di ordine 1 [ ]( ) ( )( )[ ] f x xf x f xx xf x x0 11 01 01 0, , Nel caso di 3 nodi( , , ) x x x0 1 2, si definisce la differenza divisa di ordine 2 [ ][ ] [ ]( )[ ] [ ] f x x xf x x f x xx xf x x x f x x x0 1 21 2 0 12 01 0 2 2 1 0, ,, ,, , , , Nel caso generale di n+1 nodi( , , , ..., ) x x x xn 0 1 2, si definisce la differenza divisa di ordine n [ ][ ] [ ]( )[ ] f x x x xf x x f x xx xf x x x xnn nnn 0 1 21 0 101 2 0, , , ...,, ..., , ...,... , , , ..., Ladifferenzadivisaunnumeroinvarianteallepermutazioni,ciodipendesoltantodainodienon dallordine in cui si trovano. Le differenze divise associate aglin+1 nodixi (coni=0, ...,n) si costruiscono mediante una tabella a partire dai punti noti (xi, yi=f(xi)) ordine 0ordine 1ordine 2ordine 3 ordine 4ordine n xiyi=f(xi) x0f(x0) x1f(x1)[ ] f x x0 1,x2f(x2)[ ] f x x1 2, [ ] f x x x0 1 2, ,x3f(x3)[ ] f x x2 3, [ ] f x x x1 2 3, , [ ] f x x x x0 1 2 3, , ,x4f(x4)[ ] f x x3 4, [ ] f x x x2 3 4, , [ ] f x x x x1 2 3 4, , , [ ] f x x x x x0 1 2 3 4, , , , . .... . . .... . . .... . . .... . xnf(xn)[ ] f x xn n 1, [ ]f x x xn n n 2 1, , ........................................................... [ ] f x x xn 0 1, , ..., Icoefficientidelpolinomiointerpolatorecoincidonoconglielementisulladiagonaleprincipaledella tabella cos costruita Pn(x) =a a x x a x x x x a x x x x x xn n 0 1 0 2 0 1 0 1 1+ + + + ( ) ( ) ( ) ..... ( ) ( ) ...( ) con E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 50 ( ) a f x0 0[ ] a f x x1 0 1 ,[ ] a f x x x2 0 1 2 , ,[ ] a f x x x x3 0 1 2 3 , , ,[ ] a f x x x x x4 0 1 2 3 4 , , , ,[ ] a f x x xn n0 1, , ................., Il metodo di interpolazione di Newton preferibile a quello di Lagrange perch: 1.richiede meno operazioni (O(n) a fronte di O(n2)) - minore complessit 2.ha una maggiore stabilit numerica (conseguenza di 1.) 3.bisogna memorizzare meno dati per l'interpolazione (gli elementi diagonale della tabella) 4.esistenza di un algoritmo semplice per il calcolo dei coefficienti ai del polinomio (Difdiv) 5.permigliorareilpolinomioecostruireilPn+1(x) non bisogna ripartire dallinizio, ma sufficiente aggiungere un punto noto, calcolare una nuova riga della tabella delle differenze divise, calcolare il nuovo elemento sulla diagonale e ... usarlo. Esempio 1)Datii3punti(0,0),(1,1),(2,3)costruirelatabelladelledifferenzediviseequindiilpolinomio interpolatore utilizzando il metodo di Newton. Sono dati n+1=3 punti il polinomio interpolante di grado n=2 e si scrive nella forma: P2(x) =a a x x a x x x x0 1 0 2 0 1+ + ( ) ( ) ( ) dove i coefficienti ai sono le differenze divise che si calcolano sulla diagonale della tabella: ord.0 ord.1ord.2 xi yi=f(xi) 00 11(10)/(10)= 1 23(31)/(21)= 2(21)/(20)= 1/2 Il polinomio cercato : ( ) ( )( )( )P220 1 0120 1212 2( ) x x x xxxxx x + + + + ESERCIZI SVOLTI 1.Calcolare il polinomio interpolatore per i 4 punti (0, 0) , (1, 2) , (2, -1) , (3, 0) utilizzando il metodo di Newton il metodo di Lagrange. E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 51 Sono dati n+1=4 punti il polinomio cercato di grado n=3. Ovviamente il polinomio sar lo stesso (unicit del polinomio interpolatore). Metodo di Newton Il polinomio si presenta nella forma: ( ) ( ) [ ] [ ] [ ] ( ) ( ) ( )P3 0 0 0 1 1 0 1 2 2 0 1 2 3( ) , ( ) , , ( ) , , , x f x x x f x x x x f x x x x x f x x x x + + + Si costruisce la tabella delle differenze divise a partire dai punti noti: ord.0ord.1ord.2ord.3 xiyi=f(xi) 00 12(20)/(10)= 2 21(12)/(21)= 3(32)/(20)= 5/2 30(0(1))(32)= 1(1(3))/(31)= 2(2(2.5))/(30)= 3/2 Si costruisce il polinomio utilizzando i coefficienti che compaiono sulla diagonale della tabella. ( ) P323 20 0 2 15223223211232112327152( ) ( ) ( ) x x x xx x x xx x x + + +

_,

_,

+ +

_, + Metodo di Lagrange Il polinomio si presenta nella forma: ( ) P303030303( )( )( ),,x y l x yx xx xj jjjii i jj ii i jj Sostituendo i punti noti si ottiene: E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 52 P3 00 0300 0310 1310 1320 2320 2330 3330 3301 2 30 1 0 2 0 320( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )(,,,,,,,,x yx xx xyx xx xyx xx xyx xx xx x x x xii iii iii iii iii iii iii iii i+++ + ( )( )2 31 0 1 2 1 310 1 32 0 2 1 2 300 1 23 0 3 1 3 20 22 31 1 211 32 1 1022 321 323271523 2)( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( )( ) ( )( )( )x x x x x x xx x x x x xx x x x x xx x x + + + + + + + ALGORITMI Algoritmo.Difdiv (n, x, f) Commento.I vettorix ef (dimensionen+1), inizialmente contengono i datixi eyi=f(xi) (cio glin+1 puntinotidainterpolare)coni=0,...,n.Lecolonnedellatabellaalledifferenzedivise vengonosuccessivamentedeterminateememorizzatenelvettoref.Allafineilvettoref contiene gli elementi diagonale della tabella. Parametri.Input : n, x, f Output: f 1. fi=yi per i=0, ..., n (inizializzazione fi) 2.Ciclo 1: i=1, ..., n 3.Ciclo 2: j=n, ..., i 4. ( )( )i j ji j jjx xf ff5.Fine Ciclo 2 6.Fine Ciclo 1 7.Exit End Algoritmo.Interp (n, x, f, t, p) Commento.I vettorix e f (dimensione n+1), inizialmente contengono i dati {xi} e gli elementi diagonale dellatabelladelledifferenzedivise(OutputdallalgoritmoDifdiv).Assegnatounvaloret, lalgoritmovalutailvalorecheilpolinomiointerpolatorediNewtonassumeintelo memorizza in p. Parametri.Input : n, x, f, t Output: p 1.p=fn 2.Ciclo 1: i=n1, ..., 0 3.( ) p f t x pi i + E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 53 4.Fine Ciclo 1 5.Exit End ESERCIZI PROPOSTI 1.Calcolareilpolinomiointerpolatoreperi3punti(-2,0),(0,0),(2,1),utilizzandoimetodidi Lagrange e di Newton. 18142x x +

1]1 2.Calcolareilpolinomiointerpolatoreperi3punti(-2,0),(0,6),(2,1),utilizzandoimetodidi Lagrange e di Newton. + +

1]11181462x x3.Calcolare il polinomio interpolatore per i 5 puntikk, sen4

_,

_,

con k=0, ..., 4, utilizzando i metodi di Lagrange e di Newton. [ ]00143 01144 0 0360 0 77124 3 2. . . . x x x x + +4.Calcolare il polinomio interpolatore per i 5 puntikk, exp2

_,

_,

con k=0, ..., 4, utilizzando i metodi di Lagrange e di Newton. [ ]0007379 00012248 015509 04850 14 3 2. . . . x x x x + + +E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 54 INTERPOLAZIONE POLINOMIALE A TRATTI Inpresenzadiunnumeroelevatodipuntidainterpolarepossonomanifestarsiproblemiconilmetodo della interpolazione ordinaria, che utilizza un solo polinomio interpolatore su tutto lintervallo. Ilpolinomiointerpolantedigradoprogressivamentesemprepielevatoallaumentaredelnumerodi punti e tende ad avere un profilo molto irregolare (a carattere oscillatorio), in molti casi indesiderato. Tale problema pu essere eliminato adottando una tecnica di interpolazione polinomiale a tratti, ovvero suddividendo lintervallo in sottointervalli allinterno dei quali effettuata una interpolazione con polinomi di grado inferiore. INTERPOLAZIONE LINEARE A TRATTI Nelcasopisempliceipolinomiutilizzatineisottointervallisonodiprimogrado.Taleinterpolazione lineare a tratti corrisponde a unire in successione con tratti rettilinei i punti da interpolare. xyx0= a xn= b xi= a+(i-1)hbreakpointsn+1 nodin sottointervallin-1breakpointsxyx0= a xn= b xi= a+(i-1)hbreakpointsn+1 nodin sottointervallin-1breakpoints Notouninsiemedi(n+1)punti{(xi,yi),i=0,...,n}dainterpolare,conx0=a,exn=b, si definiscono n sottointervalli dellintervallo [a,b]: [xi1 , xi],i = 1,... , n il generico sottointervallo le (n+1) ascisse xi per i = 0,... , n sono, come sempre, dette nodi; le (n1) ascisse xi per i = 1, ... , n1, di frontiera tra i sottointervalli, sono dette breakpoints. INTERPOLAZIONE CON FUNZIONI SPLINES Questotipodiinterpolazioneriproduceinformamatematicaquellochenellapraticaildisegno effettuato con il tracciacurve (flessibile di materiale plastico e modellabile a piacere). E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 55 DEFINIZIONE DI FUNZIONE SPLINE La funzione S(x) la spline di grado m associata agli (n+1) nodi xi in [a,b] con i = 0,... , n SE si verificano le seguenti CONDIZIONI: 1.S(x) un polinomio di grado m in ogni sottointervallo [xi1 , xi] per i = 1,... , n: ( )mi i i i i ix m coeff x d x c x b a x S ) 1 ( ...3 2+ + + + + + 2.ognipolinomioSi(x)elesueprimem1derivatesonocontinueintuttolintervallo[a,b]ein particolare negli (n1) breakpoints xi per i = 1,... , n1: Si (k) (xi) = Si+1(k) (xi) per i = 1 , ... , n1 per k = 0 ,... , m1 3.la S(x) interpola gli (n+1) punti: S(xi)=yi per i = 0 , ... , n Si noti che la funzione S(x) definita a tratti: S(x) = Si(x) in ogni sottointervallo[xi1 , xi] per i = 1,... , n NUMERO DI INCOGNITE PER DEFINIRE LA SPLINE Dato che ciascuna dellen splines (condizione 1.) un polinomio di grado m e quindi ha m+1 coefficienti incogniti, la funzione interpolante sullintero intervallo [a,b] univocamente definita se si determinano gli n(m+1) coefficienti incogniti totali. NUMERO DI EQUAZIONI DISPONIBILI Leequazioniutilialladeterminazionedeicoefficientiincognitisiricavanodallecondizioniimpostealla funzione S(x) per essere una spline: 2)m(n1) equazioni di continuit della funzione e delle sue derivate nei breakpoints: Si (k) (xi) = Si+1(k) (xi) per i = 1 , ... , n1 per k = 0 ,... , m1 3) (n+1) equazioni di interpolazione:S(xi)=yi QUINDI SI HANNO: n(m+1)(m1) equazioni lineari con n(m+1) coefficienti incogniti. EQUAZIONI SUPPLEMENTARI Condizionenecessariaaffinchlaspline sia univocamente determinata che il numero di equazioni sia pari al numero delle incognite, ovvero i coefficienti dei polinomi Si(x); occorrepertantoaggiungerealleprecedentiequazionilemancanti(m1)supplementari,in corrispondenza degli estremia e b dellintervallo di definizione della spline; E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 56 queste nuove equazioni sono le cosiddette (m1) condizioni al contorno, utilizzando la terminologia generalmente applicata alle equazioni differenziali. LA SPLINE CUBICA NATURALE la funzione spline pi utilizzata; dettasplinecubicaperchsibasasupolinomidigradom=3e richiedem1=2 condizioni al contorno; se le condizioni al contorno implicano lannullamento della derivata seconda agli estremi dellintervallo [a,b] si parla convenzionalmente di spline cubica naturale. CONDIZIONI che si devono verificare: 1.Si(x) un polinomio di grado 3 in ogni sottointervallo [xi1 , xi] per i = 1,... , n: ( )3 2x d x c x b a x Si i i i i+ + + 2.ogni polinomio Si(x), la sua derivata prima e la sua derivata seconda sono continue in tutto lintervallo [a,b] e in particolare negli (n1) breakpoints xi per i = 1,... , n1: Si (k) (xi) = Si+1(k) (xi) per i = 1 , ... , n1 per k = 0, 1, 2 3.la spline S(x) interpola gli (n+1) punti: S(xi)=yi per i = 0 , ... , n La spline cubica ha 4n coefficienti incogniti totali (4 in ogni sottointervallo, n sottointervalli) essendosplinecubicanaturale,le2equazionisupplementarinecessarieimplicanolannullamento della derivata seconda negli estremi a e b dellintervallo: S1 (2) (a) = 0Sn (2) (b) = 0 Il problema si riconduce alla soluzione di un sistema di equazioni lineari di ordine 4n. DIVERSA IMPOSTAZIONE possibileutilizzareunimpostazionedifferentecheportaaunsistemadiequazionilinearidiordine decisamente pi basso, pari a (n+1), e di forma pi semplice: sistema simmetrico tridiagonale a diagonale dominante. Comeincognite,anzichi4ncoefficientidellesplinecubica(4inognisottointervallo,n sottointervalli),siscelgonogli(n+1)valoriMichelederivatesecondedellasplinecubica assumono nei nodi: Si (2) = Mi per i = 0 , ... , n E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 57 Dato che la spline in esame cubica (polinomio di grado 3) la sua generica derivata seconda un polinomio di grado 1, cio una retta. In ogni sottointervallo si scrive lequazione della retta che passa per i due punti (inizio e fine) del sottointervallo [xi1 , xi] stesso: 11) 2 (11i ii ii iiM MM Sx xx x cio 11) 2 (1 i ii iiiM MM Shx x con hi ampiezza del sottointervallo. Lequazione della derivata seconda incognita risulta quindi: ( )( ) ( )ii i i iihM x x M x xx S1 1 ) 2 ( + derivata seconda Integrando lequazione della derivata seconda si ricavano le espressioni della derivata prima e della funzione: ( )( ) ( )iii i i iiChM x x M x xx S + + 221 12) 1 ( derivata prima ( )( ) ( )( )i i iii i i iiD x x ChM x x M x xx S + + + 131 136 funzione LecostantidintegrazioneCieDisiricavanoimponendoilrispettodellacondizione dinterpolazione 3.: S(xi1)=yi1 e S(xi)=yi da cui: ( )1211 166 iii ii i iii iiMhy DM M hhy yC Perchlafunzione( ) x Sisiaunasplinebisognaancoraimporreilrispettodellacondizionedi continuit 1. sulla derivata prima nei breakpoints: Si (1) (xi) = Si+1(1) (xi) per i = 1 , ... , n1 In questo modo si ricava il sistema di equazioni lineari di (n1) equazioni in (n+1) incognite (le Mi per i = 0 , ... , n): E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 58 ( ) 1 , 1 6 6 21111 1 1 1 + + ++++ + + n ihy yhy yM h M h h M hii iii ii i i i i i i(*) a cui vanno aggiunte le 2 equazioni al contorno dovute alla naturalit della spline: ( )( ) 00) 2 (10 0) 2 (1 n nM x SM x S per giungere al sistema di equazioni lineari di ordine (n+1) finale: ( )' + + ++++ + + 01 , 1 6 6 201111 1 1 10nii iii ii i i i i i iMn ihy yhy yM h M h h M hM Si nota che il sistema di equazioni (*), di ordine (n1), tridiagonale a diagonale dominante e simmetrico con matrice dei coefficienti pari a: 1111]1

+++) ( 2 ... 0 0... ... ... ...0 ... ) ( 20 ... ) ( 213 2 22 2 1n nh hh h hh h hA e pertanto risolvibile con il metodo di Gauss senza pivoting. ESERCIZI SVOLTI 1.Calcolare la spline cubica naturale interpolante a tratti i tre punti: (0, 0) , (1, 1) , (2, 0). Sono datin+1=3 puntiisottointervallisonoquindi2ein ciascun sottointervallo la spline cubica si scrive come: ( )( ) ( )( )i i iii i i iiD x x ChM x x M x xx S + + + 131 136 con: ( )1211 166 iii ii i iii iiMhy DM M hhy yC [ ]x x x i h x xi i i i i 1 11 2 , , E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 59 Trattandosi di spline naturale le condizioni al contorno sono: M M0 20 Il sistema per determinare i coefficienti della spline cubica naturale si riduce allequazione: ( ) 2 61 2 12 121 01h h My yhy yh+

_,

hh121 0 12 1 1 ;; da cuiM13 Esaminando ciascun sottointervallo: nel primo sottointervallo [(0, 0) , (1, 1)] ( )( )( ) ( )( )( ) ( ) ( )( ) [ ]Cy yhh M MD yhMS xx x M x x MhC x x Dx xx x x x11 011 1 01 02103130 03111 0 13 336326061 0 0 36 1320 012320 1 + + + + + + + ;;, nel secondo sottointervallo [(1, 1) , (2, 0)] ( )( )( ) ( )( )( )( ) ( ) [ ]Cy yhh M MD yhMS xx x M x x MhC x x Dxx x x x22 122 2 12 12213231 13222 1 23363263262 3 06 132132122323 12 + + + + + + ;;, ALGORITMI Algoritmo : Spline (n, x, y, z) E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 60 Commento. Lalgoritmo determina i coefficienti M1, M2,...,Mn-1 necessari per rappresentare la spline cubica naturale passanteperipunti(xi,yi),i=0,1,...n.Talinumerivengonomemorizzatinelvettorez.Ivettorid ecvengono introdotti per memorizzare rispettivamente la diagonale e la codiagonale del sistema tridiagonale simmetrico la cui sola soluzione fornisce i valori {Mi}. Il termine noto del sistema viene memo rizzato nel vettore b. Parametri. Input: n, x, y Output : z Ciclo 1 : i= 1,..., n-2 di = 2(xi+1-xi-1) ci = xi+1 - xi bi =6 y yx xy yx xi ii ii ii i++

_,

1111 Fine ciclo 1dn-1 = 2(xn-xn-2) bn-1 =6y yx xy yx xn nn nn nn n

_,

111 21 2 Processo di eliminazione di Gauss per il sistema tridiagonale Ciclo 2 : i = 2,..., n-1 di =di-c2i-1/ di-1 bi = bi-(ci-1 /di-1)bi-1 Fine ciclo 2 Soluzione sistema bidiagonale zn-1 =b n-1 / d n-1 Ciclo 3 : i=2,..., n-1 z n-i =(b n-1 - c n-i z n+1-i) / d n-i Fine ciclo 3 Exit End Algoritmo : Valspl (n, x, y, z, t, s, ier) Commento . Noti i coefficienti M1,M2, ...Mn-1,contenuti nel vettore z ,della spline cubica naturale passante per i punti (xi, yi), i = 0, 1,..., n, Valpl valuta il valore s che tale spline assume nel punto t ,x0 t xn. Se t[x0, xn] la variabile ier = 0, altrimenti ier=1. Parametri. Input :n, x, y, z, t Output :s, ier Ciclo 1 : i=1,.., n se xi-1 t xi poni ier = 0; vai al punto 6. Fine ciclo 1 ier = 1 Exit h= xi -xi-1 s=( ) ( )( ) ( )x t z t x zhy yhhz z t x yhzi i i ii ii i i i i + +

1]1 + 31 1311 1 1216 6 6 Exit End E. Brusa, C. Delprete, P. Gay Politecnico di Torino CALCOLO NUMERICO 61 ESERCIZI PROPOSTI 1. Calcolare la spline cubica naturale per l'insieme di dati {(0,0) (1,2), (2,0)} [3x-x3 ,-2+9x-6x2+x3] 2. Calcolare la spline cubica naturale per l'insieme di dati {(0,0), (1,2), (2,-2), (3,0)} [4x-2x3 , -6+22x-18x2+4x3 , 42-50x+18x2-2x3] 3. Determinare la spline interpolante per linsieme di datik sinkk,204

_,

_,

3212192312192312264526123 2 3 2 3 2 3x x x x x x x x x x x + + + + +

1]1, , , Politecnico di Torino CALCOLO NUMERICO 62 APPROSSIMAZIONE DI DATI Riferimento al testo: Cap. V CRITERIO DI APPROSSIMAZIONE DI DATI Ilproblemadellapprossimazioneconsistenel,assegnatim+1dati(x0,y0),...,(xm,ym),individuare unafunzionef(x),appartenenteaduncertoinsiemeoclasse,chemeglioliapprossimi,secondoun qualche criterio che andremo ad indicare. Cichedistinguelapprossimazionedeidatidallinterpolazionechelafunzioneapprossimante nonhailvincolodipassareperipuntiassegnati(dati),ovveronontenutaadassumerein corrispondenza delle ascisse dei dati {xi} il valore delle corrispondenti ordinate {yi}. Sipreferisceall'interpolazionequandoidatisonodisponibiliinnumeroelevato,eventualmente affettidaerrore(rumore).Inquestocaso,ilprocedimentodiapprossimazionemiraaridurreinparte l'effetto degli errori. SOLUZIONE DEL PROBLEMA Loscostamentotradatiefunzioneapprossimanteingeneredefinitodallanormadelvettoredegli erroriei = f ( xi ) - yi, i = 0, ..., m. Lasceltadeltipodinormaconduceaproblemidiapprossimazionediversi:nelcasodellanormal2 abbiamoilproblemadeiminimiquadratimentrenelcasodinormalsihaapprossimazione minimax. Criterio dei Minimi Quadrati (caso discreto) Notim+1punti{xi,yi} coni=0,,m, laf(x)vienesceltainmododaminimizzarelanormadue delvettoredegliscarti.Questoequivaleaminimizzarelasommadeiquadratidegliscartitrala funzione approssimante valutata nei nodi {xi} ed i corrispondenti dati { yi} mii ifny x f x f02) ) ( ( min arg ) ( Seesistonodatiacuisiattribuiscemaggioreimportanza(es.derivantidamisurazionipi precise) allora utile introdurre dei pesi wi mii i ifny x f w x f02) ) ( ( min arg ) (Criterio Minimax (caso discreto) Lafn(x)vienesceltainmododaminimizzarelanormainfinitodelvettoredegliscarti.Questo corrisponde a minimizzare lo scarto massimo i im ifny x f f ) ( max min0K

Politecnico di Torino CALCOLO NUMERICO 63 REGRESSIONE LINEARE Abbiamounproblemadiregressionelineareogniqualvoltalafunzioneapprossimantepossaessere espressa e ricercata come combinazione lineare di un insieme di funzioni di base. Piindettaglio,assegnaten+1funzionielementarikx ( ),k=0,...,n,lafunzioneapprossimante pu essere espressa come combinazione lineare queste f x c xk kkn( ) ( ) 0 MINIMI QUADRATI: CALCOLO DELL'APPROSSIMAZIONE Il metodo si applica a problemi di regressione lineare. L'errore quadratico tra i dati e la funzione approssimante in questo caso esprimibile come: ( ) R c y c xi k k iknim

1]1 20 02( )Lobiettivodelcriterioquellodideterminareicoefficientiincognitickdellafunzione approssimantecheminimizzinolerrore2.Laricercaditaleminimopuesserecondotta annullando le derivate parziali della funzione R(c) rispetto ai coefficienti ck. Derivando R(c) rispetto ai coefficienti ck si ottengono le n+1 equazioni: y c x x k ni j j ijnimk i

1]1 ( ) ( ) ,0 00 0necessarie per il calcolo dei coefficienti che definiscono la funzione f(x).Tali equazioni possono essere trascritte sotto forma di sistema di equazioni lineari 020 1 01 0 1210 1201( ) ( ) ( ) ... ( ) ( )( ) ( ) ( ) ... ( ) ( )........( ) ( ) ( ) ( ) ... ( )..x x x x xx x x x xx x x x xcciii iii n iii iiiii n iin i iin i iin ii

1]1111111 cy xy xy xni iii iii n ii

1]1111111

1]111111101( )( )..( ) che, raccolto in forma vettoriale, assume la forma Bc = d . Le equazioni costituenti tale sistema vengono definite come equazioni normali.La determinazione dei coefficienti ck , ossia lindividuazione della funzione f(x), ricondotta alla risoluzione del precedente sistema lineare. ELABORAZIONE DEL RISULTATO Ilsistemadelleequazioninormalipuessererielaboratonellaformaperagevolarelapplicazione delle tecniche di risoluzione gi note dalle precedenti esercitazioni. InparticolareilsistemaBc=dpuesserescrittocomeATAc=ATydove c il vettore colonna dei coefficientiincogniti,yilvettorecolonnadelleordinateyieAlaseguentematricedidimensioni (m+1)(n+1) : Politecnico di Torino CALCOLO NUMERICO 64 A

1]11111 0 0 1 0 00 1 1 1 10 1( ) ( ) ... ( )( ) ( ) ... ( )........( ) ( ) ... ( )x x xx x xx x xnnm m n m Il vettore c pu quindi essere espresso come c = ( AT A )-1 ATy dove la matrice (AT A )-1 AT prende il nome di pseudo-inversa di A.Questacoincideovviamenteconl'inversadiA,ogniqualvoltaArisultiesserequadrata(siricordaa taleriguardocheinognicasononpossibiledefinirelamatriceinversadiunamatrice rettangolare).SiosservicheB=ATAunamatricesimmetricadefinitapositivapercuilasoluzionedelsistema diequazionilinearipufarericorsoallafattorizzazionediCholesky,chepiefficientedella fattorizzazione LU mediante algoritmo di Gauss. INTERPRETAZIONE GEOMETRICA DEL RISULTATO EinteressanterilevarechepoichAT(y-Ac)=0(dallaATAc=ATy),ilvettorey-Acrisulta essereortogonaleallospaziogeneratodellecolonnediA(compostedaivaloriassuntidaciascuna kx ( )incorrispondenzadelleascissedeidati).L'approssimazionesecondoilmetododeiminimi quadratiequivaleadunaproiezionegeometricadelvettoredelleordinateysullospaziogenerato dalle colonne di A : la misura della distanza euclidea di y da questo spazio l'errore 2. DEGENERAZIONE DELL'APPROSSIMAZIONE IN INTERPOLAZIONE Qualorasiverifichicheilnumerodellefunzionielementarisiaugualealnumerodidatida approssimare,l'approssimazionedegeneraininterpolazione,Ainquestocasounamatrice quadratael'errorequadraticosiannulla(perchlinterpolazionegarantiscecheilvaloreassunto dalla funzione in corrispondenza dellascissa di un dato sia lordinata di quel dato). ERRORE QUADRATICO Trascrivendoilproblemadell'approssimazioneaiminimiquadratiutilizzandolamatriceA,anche l'errore quadratico pu essere riscritto in forma matriciale: ( )21

1]1y I A A A A yT T T (A) Essendo la definizione di errore quadratico medio: ( ) R c y c xi k k iknim

1]1 20 02( )(B) l'obiettivo quello di dimostrare la corrispondenza tra le due formulazioni (A) e (B). 1) Dapprima si sviluppi il prodotto in (A): Politecnico di Torino CALCOLO NUMERICO 65 ( )[ ]( ) 2 y I A A A A y y y y A A A A yT T1T T T T1T 2) si consideri che : c = ( ATA )-1ATy;inoltre cheBc = d ATAc = ATy da cui cTAT= yT 3) si riscriva l'errore quadratico medio nella forma: 2 y y y Ac y y c A AcT T T T T Quest'ultimatrascrizionecorrisponde,informamatriciale,alladefinizione(B)espressainformadi serie. ESEMPIO: TROVARE UNA RETTA APPROSSIMANTE Uncasospecificodiapplicazionedelcriteriodeiminimiquadratinelcampodellingegneriae delleconomiaconsistenelricercarelafunzionerettachemeglioapprossimiuncertoinsiemedato di punti ( xi , yi ). Questoequivaleaporrelafunzionef(x)=a+bx,conn+1=2funzionidibase, rispettivamente01 ( ) x ,1( ) x =x. Assegnati m+1 dati (xi, yi), le equazioni normali divengono: 111]1

1]1

111]1

+ iiiiiiiiiiiy xyccx xx m1021 111]1

1]1

111]1

+ iiiiiiiiiiiy xybax xx m21 Leequazioninormalisoprascritteseguonodallapplicazionedelledefinizioniedelmetodoespresse nel caso generale per descrivere la matrice B e i vettori c e d. Esempio 1)Determinareicoefficientidellarettadiregressionelineareel'errorequadratico2periseguenti dati: (0, 1), (2, 3), (3, 4), (4, 6). Le equazioni normali per la regressione lineare sono: 111]1

1]1

111]1

+ iiiiiiiiiiiy xybax xx m21 dove (m+1)=4( ) ( ) 0 11 x x x xi + + + 0 2 3 4 9; x i20 4 9 16 29 + + + ; yi + + + 1 3 4 6 14; Politecnico di Torino CALCOLO NUMERICO 66 x yi i + + + 0 1 2 3 3 4 4 6 42; 1]1

1]1

1]1

421429 99 4ba a = 0.8 b = 1.2 f(x) = 0.8 + 1.2 x L'errore quadratico medio dato da: ( ) [ ]202+ miibx a y [ ] [ ] [ ] [ ]1 08 1 12 0 3 08 1 12 2 4 0 8 1 12 3 6 0 8 1 12 4 042 2 2 2 + + + + . . . . . . . . . ESERCIZI PROPOSTI 1.Determinareicoefficientidellarettadiregressionelineareel'errorequadratico2periseguenti dati: (0, 5), (1, 3), (2, 0), (3, -1). [f(x) = 4.9 -2.1x , 2 = 0.7] 2.Determinarel'approssimazioneaiminimiquadratiel'errore2deidati:(0,2),(1,0),(2,2),(3, 6) con01 ( ) x , 1( ) x x , 22( ) x x . [f(x) = 1.9 -3.1x + 1.5x2 , 2 = 0.2] 3.Determinarel'approssimazioneaiminimiquadratiel'errore2deidati:(0,2),(1,0),(2,2),(3, 6), 01 ( ) x , 1( ) x x , 23( ) x x . [f(x) = ( 120 -97x + 22x3) /69 ,2 = 18 / 23] 4.Determinarel'approssimazioneaiminimiquadratiel'errore2deidati:(0,2),(1,0),(2,2),(3, 6) con01 ( ) x , 12 ( ) ( / ) x x sin , 22 ( ) cos( / ) x x . [f(x) = 2.5 - 3sin ( / ) x 2, 2 = 1] Politecnico di Torino CALCOLO NUMERICO 67 EQUAZIONI NON LINEARI Riferimento al testo: Cap. VI INTRODUZIONE Datalequazionef(x)=0conf(x)funzionenonlinearenelsuoargomentox,senevogliono calcolare le soluzioni o radici (valori della x che verificano luguaglianza a 0). Leradicidellequazionenonlinearef(x)=0possonononessere,ingenerale,esprimibiliinforma chiusaequindivengonocalcolatepervianumericamediantemetodiiterativiche,apartireda unaopiapprossimazioniiniziali,produconounasuccessionex0,x1,x2,...,xn,...convergente (sotto certe ipotesi) alla radice cercata. Verrannopresentatiiseguentimetodiiterativiperilcalcolodelleradicirealidiequazioninon lineari:metododibisezione,regulafalsi,metododelletangentiodiNewton-Raphsone metodo delle secanti. Sidicecheunasuccessionediapprossimazioni{xk}convergealvalorecercatoconordinedi convergenza p se esiste un numerop 1 tale che 0 lim lim1 1 + + ceex xx xpkkkpkkk, costante positiva. L'ordinediconvergenzarappresentaindicativamenteilfattoredicuiaumentanoidecimali corretti delle approssimazioni xk a partire da un certo punto in avanti. Imetodicheverrannopropostipresenterannodiversoordinediconvergenza:imetodidi bisezioneeregulafalsiconvergonoconordine1(convergenzalineare);ilmetododelle tangentiodiNewton-Raphsonconordine2(convergenzaquadratica);ilmetododellesecanti con ordine circa 1.6. METODO DI BISEZIONE Siapplicaafunzionicontinuedicuisiconosceunintervalloiniziale[a,b]chenecontieneuna radice.Sef(a)f(b) 0 nuovo intervallo dinteresse il semi-intervallo di sinistra [a, m1] se f(m1) f (b) < 0 nuovo intervallo dinteresse il semi-intervallo di destra [m1, b] se f(m1) f (b) = 0 m1 la radice cercata sef(m1)0siprendeilnuovosotto-intervalloesiricominciailprocedimento(k=k+1);aogni iterazione lampiezza dellintervallo che contiene la radice si dimezza. Lampiezza dellintervallo tende a 0 come 1/2k. Politecnico di Torino CALCOLO NUMERICO 68 Al passo iterativo k-esimo il punto medio mk=(ak-1+b k-1)/2 la stima della radice cercata. Algoritmo. Bisez (a0, b0, f, toll, x, ier) Commento. Determina,contolleranzarelativatoll,unaradicexdellequazionenonlinearef(x)=0, nellintervallo[a0,b0].Sef(a0)f(b0)0 poni ier=1; Exit 2.ier = 0, k = 0 3.k = k+1 4.mk = 12 (ak-1 + bk-1) 5.se|bk-1 ak-1| 0 nuovo intervallo dinteresse il semi-intervallo di sinistra [a, x0] se f(x0) f(b) 0, ier=1; Exit 2.Ciclo 1: n=1, ..., nmax 3.xn+1 = xnf(xn)( ) ( )x xf x f xn nn n11 4.se |xn xn+1|> xtoll|xn+1| oppure |xn xn+1|> xtoll poni ier=1 e vai al punto 10 5.se |f(xn+1)|> ftoll poni ier=2 e vai al punto 10 6.ier=0 7.x = xn+1 8.nmax = n 9.Exit 10.Fine Ciclo 1 11.x = xn+1 12.Exit End Politecnico di Torino CALCOLO NUMERICO 71 CONFRONTO TRA I DIVERSI METODI BisezioneeRegulaFalsisonometodiglobali;sonometodichiusi,ciosufficientelocalizzare la radice o meglio un intervallo iniziale che la contiene per garantire la convergenza. TangentieSecantisonometodilocali;sonometodiapertiinquantononesisteunintervallo prefissatodipartenzaerichiedonolaconoscenzadiunaodueapprossimazioniinizialineldominio dattrazionedellaradice.Ilmetododelletangentirichiedeanchelaconoscenzadellafunzione derivata f(x). Lasceltaoperativachevienesolitamenteeffettuataconsistenellutilizzodiunmetodoglobaleper avvicinarsi alla radice cercata e quindi lutilizzo di un metodo locale per il successivo raffinamento. APPLICABILITACONVERGENZA BISEZIONE ( ) f xcontinua [a, b] noto ( ) ( ) f a f b< 0 p 1 sempre REGULA FALSI ( ) f xcontinua [a, b] noto ( ) ( ) f a f b< 0 p 1 sempre TANGENTI ( ) f xcontinua e derivabile 1 approx. iniziale x0 nota p 2 soltanto se x0 vicina alla radice SECANTI ( ) f xcontinua 2 approx. iniziali x0 e x1 note p 1618 .soltanto se x0 e x1 vicine alla radice CRITERI DI ARRESTO Asecondadelvaloredelladerivatadif(x)nellaradicex,siadottaundiversocriteriodiarresto della successione di approssimazioni alla radice cercata. Sef x ' ( )>>1 conviene verificare la convergenza dif xn( ) . Il criterio di arresto del tipo f xn( ) < ftoll Ci possono essere problemi nel caso di funzioni particolarmente piatte o ripide. Sef x ' ( )

a=1.375 e b= 1.4375 ma b521375 143752140625 ++. .. EA3.110-2