Biomedical Engineering - Calcolo Numerico

Dispensa Calcolo numerico

Complete course

pag. 1 CALCOLO NUMERICO Indice Cap. 1: Ricerca di radici di funzioni non lineari ................................ ................................ ........................ 5 Metodi i terativi locali ................................ ................................ ................................ ................................ .... 5 Newton ................................ ................................ ................................ ................................ ...................... 5 Metodo delle Corde ................................ ................................ ................................ ................................ ... 5 Metodo delle Secanti ................................ ................................ ................................ ................................ . 5 Iterazioni di punto fisso ................................ ................................ ................................ ................................ . 6 Teorema: ∃! punto fisso ................................ ................................ ................................ ............................. 6 Teorema: Convergenza locale ................................ ................................ ................................ ....................... 6 Applicazione del teorema ................................ ................................ ................................ .......................... 7 Crit eri di arresto ................................ ................................ ................................ ................................ ............ 7 Criterio sul residuo ................................ ................................ ................................ ................................ ..... 8 Criterio sull’incremento ................................ ................................ ................................ ............................. 8 Cap. 2a: Metodi diretti per sistemi lineari ................................ ................................ ............................... 9 Soluzione analitica ................................ ................................ ................................ ................................ ......... 9 Metodi numerici ................................ ................................ ................................ ................................ ............ 9 Fattorizzazione LU ................................ ................................ ................................ ................................ ......... 9 Teorema 1 ................................ ................................ ................................ ................................ .................. 9 Metodo dell’Eliminazione Gaussiana (MEG) ................................ ................................ ............................... 10 Applicabilità del MEG ................................ ................................ ................................ .............................. 10 Fattorizzazione di Choleski ................................ ................................ ................................ ...................... 10 Pivoting ................................ ................................ ................................ ................................ .................... 10 Errori di arrotondam ento ................................ ................................ ................................ ........................ 10 Propagazione degli errori di arrotondamento nel MEG ................................ ................................ .......... 10 Numero di condizionamento per matrice SDP ................................ ................................ ......................... 11 Stabilità della soluzione di un sistema lineare ................................ ................................ ......................... 11 Pivoting ( II) ................................ ................................ ................................ ................................ ............... 11 Fill -in ................................ ................................ ................................ ................................ ........................ 12 Cap. 2b: Metodi iterativi per sistemi lineari ................................ ................................ .......................... 13 Definizione metodo iterativo ................................ ................................ ................................ ...................... 13 Metodo di Jacobi ................................ ................................ ................................ ................................ ......... 13 Metodo di Ga uss -Seidel ................................ ................................ ................................ ............................... 13 Metodi iterativi lineari ................................ ................................ ................................ ................................ . 14 pag. 2 Convergenza dei metodi di Jacobi e GS ................................ ................................ ................................ ... 14 Metodo di Richardson stazionario ................................ ................................ ................................ .............. 14 Stabilità e convergenza ................................ ................................ ................................ ............................ 15 Scelta parametro ottimale ................................ ................................ ................................ ....................... 15 Massima velocità di convergenza ................................ ................................ ................................ ............ 15 Tecnica di precondizionamento ................................ ................................ ................................ .................. 15 Metodo di Richardson precond izionato ................................ ................................ ................................ ... 15 Scelta precondizionatore ................................ ................................ ................................ ......................... 16 Metodo del gradiente ................................ ................................ ................................ ................................ .. 16 Scelta del parametro dinamico ................................ ................................ ................................ ................ 17 Metodo del gradiente coniuga to ................................ ................................ ................................ ................. 17 Metodi del gradiente e del gradiente coniugato precondizionati ................................ .............................. 18 Criteri di arresto ................................ ................................ ................................ ................................ .......... 18 Cap 2c: Metodi numerici per sistemi non lineari ................................ ................................ .................... 19 Il metodo di Newton ................................ ................................ ................................ ................................ .... 19 Criteri di arresto ................................ ................................ ................................ ................................ .......... 20 Costi computazionali e varianti ................................ ................................ ................................ ................... 20 Aggiornamento Jacobiana ogni p iterazioni ................................ ................................ ............................ 20 Inexact Newton ................................ ................................ ................................ ................................ ........ 20 Cap 3: Approssimazione di dati e funzioni ................................ ................................ ............................. 22 Interpolazione Lagrangiana ................................ ................................ ................................ ......................... 22 Approssimante di tipo interpolatorio ................................ ................................ ................................ ...... 22 Polinomi caratterist ici di Lagrange ................................ ................................ ................................ .......... 22 Interpolatore Lagrangiano ................................ ................................ ................................ ...................... 22 Accuratezza dell’interpolatore Lagrangiano ................................ ................................ ........................... 23 Convergenza de ll’interpolatore Lagrangiano ................................ ................................ .......................... 24 Fenomeno di Runge ................................ ................................ ................................ ................................ . 24 Interpolazione Lagrangiana composita ................................ ................................ ................................ ....... 24 Accuratezza dell ’iterpolatore Lagrangiano composito ................................ ................................ ............ 24 Scelta di H e k ................................ ................................ ................................ ................................ ........... 25 Interpolazione con ubicazione specifica dei nodi ................................ ................................ ........................ 25 Nodi di Chebyshev ................................ ................................ ................................ ................................ .... 25 Interpolatore Lagrangiano sui nodi di Chebyshev ................................ ................................ ................... 25 Convergenza dell’interpolazione sui nodi di Chebyshev ................................ ................................ .......... 26 Nodi di Chebyshev su un intervallo generale ................................ ................................ ........................... 26 Interpolazione trigonometrica ................................ ................................ ................................ .................... 27 Interpolatore trigonometrico ................................ ................................ ................................ ................... 27 Spazio fisico e spazio delle frequenze ................................ ................................ ................................ ...... 28 pag. 3 Fast Fourier Transform ................................ ................................ ................................ ............................ 28 Espressione Lagrangiana dell’inter polatore trigonometrico ................................ ................................ ... 28 Accuratezza ................................ ................................ ................................ ................................ .............. 29 Il fenomeno dell’aliasing ................................ ................................ ................................ .......................... 29 Metodo dei minimi quadrati ................................ ................................ ................................ ....................... 29 Estrapolazione di dati ................................ ................................ ................................ .............................. 29 Relazione analitica semplice ................................ ................................ ................................ .................... 29 Minimi quadrati ................................ ................................ ................................ ................................ ....... 29 Retta di regressione ................................ ................................ ................................ ................................ . 30 Cap 4: Integrazione numerica ................................ ................................ ................................ ............... 32 Formule di quadratura (composita) ................................ ................................ ................................ ............ 32 Formula del punto medio composita ................................ ................................ ................................ ....... 32 Formula dei trapezi composita ................................ ................................ ................................ ................ 33 Formula di Simpson composita ................................ ................................ ................................ ................ 33 Cap 5: Approssimazione numerica di equazioni differenziali ordinarie ................................ ................... 35 Problema di Cauchy ................................ ................................ ................................ ................................ ..... 35 Problema di Cauchy scalare ................................ ................................ ................................ ..................... 35 Esistenza e unicità della solu zione continua ................................ ................................ ............................ 35 Approssimazione di derivate ................................ ................................ ................................ ....................... 35 Approssimazione numerica del problema di Cauchy ................................ ................................ ............... 36 Metodi di Eulero in avanti e all’indietro ................................ ................................ ................................ ...... 36 Metodo in avanti ................................ ................................ ................................ ................................ ..... 36 Metodo all’indietro ................................ ................................ ................................ ................................ .. 37 Soluzione dei metodi di Eulero ................................ ................................ ................................ ................. 37 Aspetti computazionali dei metodi di Eulero ................................ ................................ ........................... 37 Assoluta stabilità ................................ ................................ ................................ ................................ ......... 38 Utilizzabilità dei metodi di Eulero ................................ ................................ ................................ ............ 38 Assoluta stabilità ................................ ................................ ................................ ................................ ..... 38 Analisi dell’assoluta stabilità per il metodo di Eulero in Avanti ................................ ............................... 39 Assoluta stabilità del metodo di Eulero all’Indietro ................................ ................................ ................. 39 Prob lemi generali ................................ ................................ ................................ ................................ ..... 39 Metodo di Crank -Nicolson ................................ ................................ ................................ ........................... 40 Metodo di Heun ................................ ................................ ................................ ................................ ........... 40 Convergenza ................................ ................................ ................................ ................................ ................ 41 Accuratezza della soluzione numerica ................................ ................................ ................................ ..... 41 Ordine di convergenza ................................ ................................ ................................ ............................. 41 Proprietà di convergenza e stabilità dei metodi visti ................................ ................................ .............. 41 Consistenza ................................ ................................ ................................ ................................ .............. 41 pag. 4 Assoluta stabilità e convergenza ................................ ................................ ................................ ............. 42 Eulero in Avanti o Eulero all’Indietro? ................................ ................................ ................................ ..... 42 pag. 5 Cap. 1: Ricerca di radici di funzioni non lineari Problema: data f∶(a,b)→ R, trovare α ∈(a,b),f(α)= 0, dove α è detta zero di f o radice dell’equazione f(α)= 0. Metodi itera tivi locali Scelgo x(0); genero x(1),x(2),… ,x(k),… (successione); obiettivo limx→k x(k)= α (convergenza). Newton Considero la tangente a f in x(k): y(x)= f(xk)+ f′(xk)(x− x(k)) e scelgo x(k+1) tale che y(x(k+1))= 0, da cui si ottiene: x(k+1)= x(k)− f(x(k)) f′(x(k)) , con k≥ 0. Quindi considero l’approssimazione della tangente al punto x(k+1). Inoltre non posso usare Newton se f non è continua o se il punto x(k) è un punto stazionario. Pe r questo metodo dunque s erve la derivata prima . Si può quindi evitare usando il Metodo delle Corde o il Metodo delle Secanti. Metodo delle Corde Corde: x(k+1)= x(k)− f(x(k)) q , dove q è la pendenza della retta che unisce f(a) e f(b), quindi v ale q= f(b)−f(a) b−a ≠ 0. Quindi approssimo la derivata prima di f in x(k) con q e di conseguenza x(k) nuova è lo zero della retta e x(k+1) è la radice della retta parallela pas sante per f(x(k)) Metodo delle Secanti È simile al metodo di Ne wton ma sostituisce la derivata prima con il r apporto incrementale, ottenendo così: x(k+1)= x(k)− (x(k)− x(k−1))∙ f(x(k)) f(x(k))−f(x(k−1)) con k≥ 1. Questo metodo è più usato, rispetto ai precedenti, perché è più dinamico poiché consider a anche il valore precedente, quindi ha memoria, x(k−1). I 3 metodi precedentemente visti sono me todi “locali", cioè la convergenza è garantita solo se x(0) è preso sufficientemente vicino a d α, ossia se ∃�> 0∶|������− x(0)|< �. La precedente condiz ione è difficile da ap plicare a priori (δ è difficilmente stimabile) allora si possono usare consider azioni fisiche. Ma α e δ non li conosco perché α è il valore che sto approssimando e δ è l’intorno, quindi devo ipotizzare x(0) tale che ricada nell’int ervallo affinché converga , quindi x(0) deve essere dentro l’intorno. pag. 6 Iterazioni di punto fisso Detta Ф(x) la funzione di cui si cercano i punti fissi, cioè le soluzioni dell’equazione x= Ф(x). Il punto fisso è dato dall’intersezione tra la Ф(x) e la bisettrice del 1° -3° quadrante. Quindi il punto fisso si ha quando α= Ф(α). Posso riscrivere il prob lema come: Un metodo per la ricerca dei punti fissi di una funzione Ф ∶[a,b]→ R lo si può tro vare considerando il seguente metodo iterativo: Ф funz ione di iterazione di punto fisso. Se sto considerando f, cerco le radici, se sto considerando Ф, cerco i punti fissi. I metodi visti in precedenza possono essere visti come iterazioni di punto fisso : In generale non è vero che un’iterazione di punto fisso è sempre convergente, per esempio l’esponenziale semplice non ha intersezioni con la bis ettrice, e quindi le iterazioni di punto fisso non convergono. Teorema: ∃! punto fisso Sia Ф ∶[a,b]→ R continua e dato x(0) in [a,b]. Consideriamo le iterazioni di punto fisso: x(k+1)= Ф(x(k)) con k≥ 0. 1. Se ∀x ϵ [a,b] si ha Ф(x) ϵ [a,b] allora ∃ almeno un punto fisso α ϵ [a,b]. 2. Se inoltre ∃L< 1 tale che |Ф(x1)− Ф(x2)|≤ L|x1− x2| ∀x1,x2 ∈[a,b] allora: {come trovo un intervallo [a, b] tale che valga la proprietà 2 richiesta dal teorema?} Teorema: Convergenza locale Sia Ф ∶[a,b]→ R funzione di iterazione e supponiamo che sia α un suo punto fisso. Sia inol tre Ф ∈�1(������������), dove ������������ è un opportuno intorno di α. pag. 7 Si ha una convergenza più rapida in ordine 2 Applicazione del teorema Criteri di arresto Supponiamo che {x(k)}k sia una successione di approssimazione di uno zero di f(x) generata da uno dei metodi visti precedentemente. Supponiamo inoltre che f∈C1(Iα) ∀k e indichiamo con e(k)= α− x(k) l’errore al passo k. Si vuole dunque sapere quando arrestare il metodo iterativo. Esistono due criteri: Criter io sul residuo e Criterio sull’incremento. pag. 8 Criterio sul residuo Fissato ε> 0, si arresta il met odo se |f(x(k))|< �. Si può dimostrare che: Due possibili situazioni: |f′(α)|≪ 1 (a sinistra) in cui il criterio di non è valido perché si h a un errore troppo grande, |f′(α)|≫ 1 (a destra) per cui il criterio di arresto è valido ma è troppo restrittivo. Criterio sull’incremento Fissato ε> 0, si arresta il metodo se |x(k+1)− x(k)|< �. Si può dimostrare che: pag. 9 Cap. 2a: Metodi dire tti per sistemi lineari Considero il sistema A�= � lineare dove A∈Rn×n di componenti aij e �∈Rn sono assegnati e �∈Rn è il vettore delle incognite. Queste sono n equazioni lineari nelle incognite xj e l’i -esima equazione è: Soluzione analitica La soluzione del sistema esiste unica se det(A) ≠0, la sua espres sione è data da: (Cramer) con Aj= [a1… aj−1baj+1… an] e ai le colonne di A. Tuttavia essa è inutilizzabile poiché il calcolo richiede circa n! operazioni, le quali richiedono un tempo troppo elevato. Metodi numerici Servono quindi metodi numerici che si traducano in algoritmi efficienti, partendo da un caso semplice: la matrice triangolare inferiore (������� = �) lij= 0,per i< j, in generale quindi si ha: ed è detto Metodo di sostituzioni in avanti. Per tale metodo sono necessarie circa n2 operazioni. In maniera analoga si può risolvere un sistema con la ma trice triangolare superiore ( ������� = �) uij= 0,per i> j, in generale quindi si ha: ed è detto M etodo di sostituzioni all’indietro. Per tale metodo sono necessarie circa n2 operazioni. Fattorizzazione LU Definiamo Ap,i∈Ri×i,i= 1,… ,n le sottomatrici principali di A ottenute considerando le prime i -righe e i - colonne. Teorema 1 Data A∈Rn×n, se det( Ap,i)≠ 0 ∀i= 1,… ,n− 1 allora esistono un’unica matrice triangol are inferiore L e un’unica matrice triangolare superiore U tali che: A= LU (fattorizzazione LU). Il sistema lineare originale può essere scritto equivalentemente come cascata di due sistemi lineari triangolari: y è una nuova variabile ausiliaria. Il van taggio è quello di ottenere due matrici triangolari (superiore e inferiore) che posso risolvere coi metodi precedenti. pag. 10 Metodo dell’Eliminazione Gaussiana (MEG ) Si introducono m atrici A(k) e termini noti �(k) ad ogni passo siano A(1)= A,aij(k)= (������(k))ij. L’idea è quella di annullare gli elementi aij(k) con i≥ k e j< k. Per fare ciò si sottrae un multiplo della riga k alle successive in modo da annullare gli elementi desiderati. Al termin e del MEG si ottiene: A(n)= U lij→ L �(n)= � Questo metodo , MEG + sostituzioni indietro, ha un costo computazionale di circa 2 3n3 ⁄ . Applicabilità del MEG La verific a a priori delle ipotesi di applicabilità del Teo rema 1 è troppo costosa per essere utile in generale (richiede calcolo di n -1 determinanti). Esistono condizioni sufficienti facili da verificare: 1. Matrici a dominanza diagonale stretta per righe o per colonne (il modulo dell’elemento diagonale è maggiore della somma degli elementi della riga o della colonna); 2. Matrici simmetriche definite positive Se queste condizioni sono verificate allora posso applicare il teorema, se invece non lo sono allora non posso dire nulla. Fattorizzazione di Choleski Per matr ici simmetriche definite positive vale la seguente fattorizzazione di Choleski (che è unica ed è un caso particolare di MEG): A= RTR, con R triangolare superiore. Pivoting Se al passo k del MEG si ha akk(k)= 0 (cioè se det( Ap,k)= 0), scambio riga k con riga i>k tale che aik(k)≠ 0. Il pivoting può essere interpretato come una pre -moltiplicazione di A e b per una matrice di permutazione P e quindi risulta: Quindi ogni passo del MEG ha la sua matrice di permutazione P(k). In generale, se al passo k scambio riga k con riga r>k , sto di fatto pre -moltiplicando A(k) per una matrice ottenuta scambiando nella matrice identità la k -esima e la r -esima riga. Errori di arrotondamento Il MEG non è affetto da errore numerico, restituisce la soluzione esatta. Questo però è vero in arit metica esatta cioè svolgendo i conti “a mano”. Il calcolatore, nel memorizzare la v ariabile z, introduce un errore di arrotondamento: ẑ= z+ δz, dove ẑ è l’effettivo valore memorizzato (aritmetica finita). Questo avviene perché il calcolatore è dotato di memoria finita e non può quindi memorizzare esattamente i numeri. Propagazione degli e rrori di arrotondamento nel MEG In un algoritmo, d ove un numero enorme di operazioni elementari (addizione, sottrazione, moltiplicazione, divisione) viene eseguito, gli errori di arrotondamento (seppur piccolissimi) propagano, possibilmente amplificandosi, rendendo inaccurato l’output. N el nostro caso sti amo risolvendo il sistema lineare perturbato: (A+ δA)�̂= �+ δ� e introducendo un errore computazionale δ�= �̂− �. Si può dimostrare che sotto l’ipotesi ||δA||2||A−1||2< 1 vale il seguente risultato: pag. 11 Dove K2(A)= ||A||2||A−1||2 è il numero di condizionamento della matrice co n la norma euclidea e Ponendo si ottiene questo andamento: L’errore computazionale dovuto alla propagazione dell’errore di arrotondamento è quindi grande per matrici mal condizionate. Numero di condizionamento per matrice SDP Per una matrice sim metrica e definita positiva si ha: poiché vale (per una generica matrice reale) Si ottiene: valore del numero di condizionamento. Stabilità della soluzione di un sistema lineare In generale possia mo introdurre il residuo ������= ������− �������̂. Essendo una quantità calcolabile (a differenza dell’errore), ci chiediamo come sia legato all’errore: sottraendo l’espressione del residuo dal sistema lineare originario: avendo notate che la norma del prodotto è mino re o uguale del prodott o delle norme, notiamo che: attraverso dei passaggi matematici otteniamo: Da questa disequazione si deduce che il metodo del residuo è buono per K2(A) piccolo, perché se il residuo è piccolo, ma il condizionamento è grande, si ha un errore grande. Pivoting (II) L’operazione che maggiormente amplifica l’errore di arrotondamento è la sottrazione. Riprendendo il MEG, vorremmo che akk(k) sia il più grande possibile in modo da rendere più piccolo il sottraendo nella sottrazione. Ciò vale in aritme tica finita: Dove δaik(k) è l’errore di arrotondamento. Prendendo âkk(k) il più grande possibile, riduco anche l’effetto di δaik(k). per questo motivo, pivotin g si fa sempre , per limitare la propagazione dell’errore, (anche quando akk(k)≠ 0) scambiando la riga k con la riga r̅∶|ar̅k(k)|= maxi>k|aik(k)|. Nel caso più generale, il pivoting viene effettuato anche scambiando le colonne nella ricerca del massimo. Ciò equivale a introdurre un’altra matrice di permutazione Q: PAQ = LU Il pivoting, che in questo caso è totale, riduce al minimo l’errore computazionale, infatti per avere akk(k) molto grande si fa pivoting anche per le colonne e si riduce quindi ulteriormente l’errore. Niente però può fare per matrici mal co ndizionate perch é in generale anche K2(PAQ ) è grande. pag. 12 Fill -in Il pattern di una matrice mi permette di identificare facilmente gli elementi diversi da zero. Si ottiene ponendo un punto o una croce su ogni posizione con elementi non nulli. Una matrice è detta sparsa quand o il numero di elementi nulli è molto maggiore di quello degli elementi non nulli. Il fill -in consiste nel fatto che i fattori L e U di una matrice sparsa non sono in generale sparsi. L’implementazione del MEG è fatta in modo da memo rizzare L e U sovrasc rivendo lo spazio dedicato ad A. Infatti la fattorizzazione LU non è adatta per matrici sparse (dal punto di vista dell’occupazione di memoria). Per ridurre il fill -in può essere utile eseguire una tecnica di riordinamento, la quale co nsiste nel numerare diversamente le righe di A: Prima del riordino Dopo il riordino pag. 13 Cap. 2b: Metodi iterativi per sistemi lineari Riprendiamo il sistema lineare A�= � con A∈Rn×n , �∈Rn, �∈Rn e det(A) ≠0. La fatt orizzazione LU in generale n on è idonea per: - Sistemi di grandi dimensioni, visto il costo di circa n3 operazioni; - Sistemi con matrici sparse a causa del fenomeno del fill -in. Definizione metodo iterativo Si introduce una successione �(k) di vetto ri determinata da una legge ricorsiva che identifica il metodo. Al fine di innescare il processo iterativo, è necessario fornire un punto di partenza (vettore iniziale) �(0) sulla base, ad esempio, di aspettative fisico/ingegneristiche: �(0)→ �(1)→ ⋯ → �(k)→ ⋯ Affinché il metodo abbia senso, deve soddisfare la proprietà di convergenza: (la convergenza non deve dipendere dalla scelta di �(0)) Poiché la convergenza è garantita solo dopo infinit e iterazioni, dal punto di v ista pratico dovremo arrestare il processo iterativo dopo un numero finito di iterazioni, quando riterremo di essere arrivati “sufficientemente vicini” alla soluzione. Possiamo quindi concludere che anche in aritmetica esatta, un metodo iterativo (a diffe renza del MEG) sarà inevitabilmente affetto da un errore. Metodo di Jacobi Partiamo dalla i -esima riga del sistema lineare: formalmente la soluzione ������� è data da: ov viamente la formula precedente non è utilizzabile perché non conosciamo le altre soluzioni xj; potremmo pensare di introdurre un metodo iterativo che aggiorni al passo k+1 la xi utilizzando le altre xj ottenute al passo prec ed ente . Ogni iterazione costa n2 operazioni, quindi Jacobi è competitivo con il MEG se il numero delle iterazioni è inferiore ad n. Per matrici sparse, il costo è invece di solo circa n operazioni per iterazione. In questo caso Jacobi è competitivo con il MEG se il numero di iterazioni è lineare rispetto n. Si nota come le soluzioni xi(k+1) possono essere calcolate in parallelo, quindi con grandi vantaggi per sistemi di grandi dimensioni. Metodo di Gauss -Seidel Partendo dal metodo di Jacobi, pre ndiamo j>i. Al passo k+1 conosciamo già la xj(k+1) perché è già calcolata. Possiamo quindi pensare di usare le quantità al passo k+1 se ji: I costi computazionali sono paragonabili a quelli del metodo di Jacobi; però, a differenza di Jacobi, con GS non è possibile calcolare le soluzioni in parallelo, cioè GS si usa in serie, cioè calcola in sequenza xi, poi xi+1 e così via. pag. 14 Metodi iterativi lineari In generale consideriamo metodi iterat ivi lineari della seguente forma: �(k+1)= B�(k)+ � dove B∈Rn×n (matrice di iterazione) e �∈Rn identificano il metodo. Si prendono B ed f considerando: 1. Consistenza: il metodo se viene applicato alla soluzione esatta x deve restituire x (quindi non ci deve essere aggiornamento se si arriva alla soluzione esatta): La formula precedente dà una relazione tra B ed f in funzione dei dati e fornisce una con dizione necessaria e non sufficiente per la convergenza; 2. Convergenza: introduciamo l’errore �(k)= �− �(������) e una norma vettoriale ||∙||, attraverso diversi passaggi matematici si ottiene: se ||B||< 1 allora limk→∞|| �(k+1)||= 0 L’errore quindi va a zero se ||B||< 1, e questa è detta condizione di stabilità. Come norma di mat rice si è soliti prendere il raggio spettrale: con gli autovalori di B Convergenza dei metodi di Jacobi e GS Siano: D la diagonale di A; -E la sua parte triangolare inferiore; -F la sua parte triangolare superiore. E’ facile vedere che il metodo d i Jacobi può essere scritto come: D�(������+�)= (E+ F)�(������)+ � e quindi la matrice di iterazioni è data da: Per il metodo di Gauss -Seidel invece si ha: (D− E)�(������+�)= F�(������)+ � e matrice di iterazione: Entrambi i metodi sono consistenti perché D− E− F= A e A�= �. Valgono qu indi i seguenti risultati di convergenza: - Se A è strettamente dominante diagonale per r ighe, |aii|> ∑ |aij| j≠i ,i= 1,… ,n allora J eGS convergono ; - Se A è simmetrica e definita positiva, allora GS converge; - Se A è tridiagonale, allora GS e J o conv ergono entrambi o non convergono entrambi. Se convergono, allora GS è più veloce. Metodo di Rich ardson stazionario Il metodo di Richardson è basato sulla seguente legge di aggiornamento, assegnati �(0) e α∈R: L’idea è di aggiornare la soluzione nu merica aggiungendo una quantità proporzionale al residuo. Infatti, ci si aspetta che se il residuo è grande (piccolo) bisogna correggere di molto (poco) la soluzione al passo k . Poiché dalla sua definizione risulta �(k+1)= (I− αA)�(k)+ α�, segue che il metodo di Richardson stazionario è caratterizzato dalla matrice di iterazione Bα= I− αA e da �= α�. È un metodo consistente, infatti: �= �+ α(�− A�). pag. 15 Stabilità e convergenza Concentriamoci sul caso di matrice A simmetrica e definita positiva. Int roduciamo gli autovalori di A che sono reali e positivi: Teorema: il metodo di Richardson stazionario con matrice A simmetrica e definita positiva è convergente se e solo se: Sce lta parametro ottimale Ci chiediamo ora quale sia il valore del parametr o, fra quelli che garantiscono la convergenza, che massimizza la velocità di convergenza. Introduciamo la seguente norma indotta da A: ||������||������= √∑ ����������������� ��,�=1 Equivalentemente s i ha: ||������||������= √������������������� , ricordando che ||�(�+1)||������≤ ������(�������)||�(�)||������, cerchiamo: Tale che ������(�������) sia minimo, cioè che a parità di errore al passo k, m inimizzi l’errore al passo k+1: Mettiamo in ascissa i valori di α e in ordinata i valori |1− αλi(A)|. Da questo grafico, si nota come il valore di α ottimale sia dato dall’intersezione fra le curve |1− αλ1(A)| ascendente e |1− αλn(A)| discendente ; quindi si ottiene: I valori ottimali per cui si ha convergenza quindi sono tra 0 e 2 λ1 ⁄ . Massima velocità di convergenza Calcoliamo il valore del raggio spettrale in corrispondenza del valore del parametro ottimale: Poiché A è simmetrica e definita positiva, abbiamo ||A||= λmax (A); inoltre vale: Siccome anche A−1 è SDP, abbiamo che (attraverso semplici passaggi): Tecnica di precondizionamento Il valo re ottimale ρopt esprime la mass ima velocità di convergenza possibile ottenibile per una matrice A con il metodo di Richardson stazionario. Le matrici mal condizionate (K(A) grande) sono caratterizzate da una velocità di convergenza del metodo molto bassa. Per migliorare la velocità di c onvergenza introduciamo una matrice P simmetrica definita positiva e invertibile. Allora il sistema lineare di partenza è equivalente (nel senso che ammette la stessa soluzione) al seguente sistema precondizionato: Metodo di Richardson p recondizionato Applichiamo il metodo di Richardson stazionario al precedente sistema: pag. 16 Otteniamo gli stessi risultati di convergenza nel caso non precondizionato, a patto di sostituire A con P−1A (anch’essa di autovalori reali positivi): - Convergenza: - Valori otti mali: Quindi, se K(P−1A)≪ K(A) otteniamo una velocità di convergenza massima maggiore rispetto al caso non precondizionato. Non possiamo affermare che un tale metodo precondizionato è più veloce del caso non precondizionato. Perché per il momento possi amo solo dire che il numero delle i terazioni necessarie per soddisfare un certo criterio di arresto è minore, nel caso precondizionato K(P−1A)≪ K(A). L’effettivo tempo di calcolo per risolvere il sistema lineare è dato da: CPUtime = n.iter ∙C, dove C è il c osto della singola iterata. Introduci amo ora il residuo precondizionato: z(k)= P−1r(k). Un effettivo algoritmo è dato quindi dai seguenti passi ad ogni k: Al terzo passo dobbiamo risolvere un sistema lineare in P (assen te nel caso non precondizionato). Scelta precondizionatore 1. K(P−1A)≪ K(A) Come già visto, questo permette di ridurre il numero delle iterazioni. A tal fine, P deve poter ridurre di ampiezza l’intervallo che contiene gli autovalori di A. 2. Il sistema linear e in P deve essere faci lmente risolvibile Questo permette di avere C non troppo elevato. A tal fine, P deve avere una struttura speciale, ad esempio essere diagonale, triangolare o data dal prodotto di esse. Metodo del gradiente Si introduca la seguente funzione energia: Ф ∶ Rn → R: Se A è SDP, l’energia è una funzione convessa che ammette un unico punto di minimo. Poiché ∇Ф(�)= A�− � , abbiamo che il punto di minimo (gradiente nullo) coincide con la soluzione di A x=b . Dato un punto ������̅, il vettore −∇Ф(������̅) indivi dua la direzione di massima decrescita dell’energia (ortogonale all’isolinea con energia pari a Ф(������̅)). Quindi se traccio dei piani equipotenziali, intersecati alla funzione di energia, ottengo delle ellissi che hanno energia crescente da x a ������̅. Il m etodo del gradiente assume come direzione di aggiornamento quella di massima decrescita dell’energia nel pun to ������(�): dove ������� è un parametro dinamico (ossia cambia ad ogni iterazione k). Sappiamo che: e quindi il metodo del gradiente può ess ere interpretato come un metodo di Richardson con parametro dinamico. pag. 17 Il vantaggio di questo metodo è che si può scegliere il parametro ������� in maniera ottimale ad ogni iterata e soprattutto che non richiede più di conoscere gli autovalori di A (che spe sso vanno approssimati essi stessi per via numerica). Scelta del parame tro dinamico Consideriamo l’energia nel nuovo aggiornamento: Il parametro ottimale viene scelto in modo da minimizzare l’energia fra tutte le possibili ottenute muovendosi in dire zione del gradiente: Notiamo che: quindi il prodotto matrice -vettore nel calcolo di αk viene usato nel calcolo del residuo. [È vero che il prodotto matriciale viene fatto come passo in più rispetto a prima, però nel calcolo del residuo non va ripetuto perché è già stato fatto. Inoltre è pi ù conveniente di Richardson perché ha un alfa di namico contro uno statico]. Metodo del gradiente coniugato La velocità di convergenza del metodo del gradiente è la stessa di Richardson stazionario ottimale: dove, dato un vettore z, abbiamo utilizzato la norma in A che, ricordiamo, è data da: [La norma in A è utile rispetto al raggio spettrale perché tiene conto della matrice di A] La direzione di aggiornamento data dal residuo è ortogonale alla precedente per costruzion e. Infatti il gradiente è sempre ortogonale alla direzione precedente, cioè La scelta della direzione di massima decrescita è quindi ottimale fra un’iterata e l’altra ma potrebbe non esserlo a medio/lungo termine. Infatti, in generale si ha: Il carattere esplorativo di Rn da parte del gradient e non è quindi elevato. Per superare tale limite, si introducono nuove direzioni di aggiornamento �(k) tutte ortogonali (rispetto al prodotto scalare in A) fra di loro: (direzioni coniugate) Il metodo del gradiente non garantisce di essere perpendicolare a �(�) ma solo a �(�+1). In questo modo si esplora con più efficac ia lo spazio. Il gradiente è ottimale in locale, mentre il gradiente coniugato permette di trovare la direzione ottima. pag. 18 Se A è SDP, l’errore fra l’iterazione k e quella di partenza scala nel seguente modo : con per il metodo del gradiente s i ha: si deduce quindi che per matrici SDP la convergenza col metodo del gradiente coniugato risulta più veloce. Questo perché il numero di condizionamento viene smorzato dalla radice. [I costi computazionali della singola iterazione sono paragonabili a Richardson e al gradiente, ed è efficace per matrici sparse]. Metodi del gradiente e del gradiente coniugato precondizionati Entrambi i metodi soffrono in ogni caso per sistemi con matrici malcondizionate allora si introduce un precondizionatore con gli stessi requisiti discussi in precedenza, al fine di abbattere il numero di iterazioni . In entrambe i casi se A e P sono delle SDP , si considera il problema precondizionato �̂������̂= �̂ con: In questo caso, l’errore si comporta come quello del caso non precondizionato a patto di sostituire con ������−1. Ad esempio per il gradiente coniugato precondizionato , se il precondizionatore è scelto bene si ha: e quindi l’erro re nel caso precondizionato si riduce più velocemente rispetto al caso non precondizionato . Criteri di arresto Come discusso, servono criteri per arrest are il processo iterativo. Ciò deve avvenire quando stimiamo che l’errore �(�)= ������− ������(�) (che non possiamo calcolare non conoscendo x) sia sufficientemente piccolo. 1. Criterio sul residuo: è un buon criterio per bassi condizionamenti , ed è utilizzato s pesso per metodi precondizionati. 2. Criterio sull’incremento: è un buon metodo se ρ(B) è piccol o, e non molto buono se ρ(B) è molto vicino a 1. pag. 19 Cap 2c: Metodi numerici per sistemi non lineari Spesso i problemi ingegneristici sono non lineari. La loro approssimazione numerica porta a dover risolvere un sistema di equazioni non lineari. In gen erale abbiamo il seguente problema: Dove �∈Rn è il vettore delle incognite xj ,�∶ Rn→ Rn con le fi che esprimono delle relazioni non lineari fra le incognite. Indichiamo con α le radici di f. [2 difficoltà dal punto di vista num erico: 1. È un sistema, quindi ho equazioni accoppiate 2. È non lineare, quindi va linearizzato ] Il metodo di Newton Riprendiamo il metodo di Newton introdotto per il caso scalare, cioè di determinazione delle radici di una funzione non lineare: Questo può es sere scritto in maniera equivalente nel seguente modo: Dato un punto z∈Rn , introduciamo la matrice Jacobiana relativa a f: Per ogni z∈Rn sì ha quindi J(z)∈Rn×n In analogia con il metodo di Newton introdotto per trovare le radici della funzione f, introduciamo il metodo di Newton per la risoluzione numerica del s istema non li neare �(x)= �, cioè per la determinazione delle sue radici α: Dato x(0), se J(x(k)) è una matrice non singolare, per ogni k risolvo: Ad ogni iterazione k, il primo passo del metodo di Newton consiste nella risoluzione di un sistema lineare di dimensione n 1. sistema lineare nxn Infatti J(x(k)) è una matrice non singolare di coefficienti noti e −f(x(k)) è noto. Una volta risolto il sistema lineare e ottenuto quindi δ�(k) aggiorno la �(k+1) mediante il s econdo passo del metodo. 2. aggiorno tutto a d ogni iterazione k [È come risolvere A x = b dove A e b sono noti; ma in questo caso A e b cambiano ad ogni iterazione k, quindi si ha la concatenazione tra metodo di Newton e un metodo numerico. L’aggiornamento costa circa n] Analogament e al caso scalare, il metodo di Newton per sistemi è: - Un metodo locale: allora si ha convergenza pag. 20 - Un metodo del secondo ordine (scalatura secondo una parabola) : Se converge e J è derivabile, allora si ha: Criteri di arres to Essendo un metodo iterativo, bisogna introdurre opportuni criteri di arresto. Analogamente al caso scalare, si introducono: - Criterio sul residuo: - Criterio sull’incremento: Costi computazionali e varianti Il metodo di Newton richiede la risoluzi one di un sistema line are ad ogni iterazione. Il suo costo computazionale è quindi dato da: Dove ��������� è il tempo necessario per costruire la matrice J(x(k)) e ��� il tempo di risoluzione del corrispondente sistema lineare. Notiamo che J(x(k)) cambia ad ogni iterazione k e quindi va ricostruita ogni volta . In generale, #�������������� è mol to basso essendo Newton un metodo di ordine 2. Tuttavia, potrebbe essere molto gran de a causa di un tempo di costruzione e/o di risoluzione del sistema lineare molto alto . Per velocizzare il metodo di Newton si possono introdurre le seguenti varianti . Ag giornamento Jacobiana ogni p iterazioni Dato un numero �≥ 2, l’idea di aggiornare la matrice Jacobiana J(x(k)) solo ogni p volte. Il tempo di costruzione ��������� è quindi ridotto, dovendosi effettuare solo ogni p iterazioni. Mediamente si ha quindi di fatto ��������� = ��������� −������������������� /� ; dove ��������� −������������������� è il temo di costruzione della Jacobiana che si avrebbe con Newton. Se inoltre si utilizza la fattorizzazione LU per risolvere il sistema lineare, allora il costo della fattorizzazione di J(x(k)) è spalmato su p iterazioni. Di fatto a d ogni iterazione medi amente ho un costo di circa 2n3/3p . Il tutto al prezzo di un numero di iterazioni maggiori di quelle che si hanno con Newton. Infatti il metodo non è più del secondo ordine: #�������������� > #�������������� ������������������� Più p è grande più aument a il numero di iterazioni, n on deve essere però troppo grande altrimenti si perde la convergenza. La speranza è che il costo complessivo sia ridotto, specialmente nei casi dove la costruzione della matrice Jacobiana è molto oneroso: Inexact Newton L’ide a è di sostituire ad ogni iterazione k la matrice Jacobiana J(x(k)) con una sua approssimazione J̃k che sia in generale facilmente costruibile e tale che il sistema lineare associato sia facilmente risolvibile . Il tutto ancora al prezzo di numero di iterazioni maggiori di quell e che si hanno con Newton (ordine< 2). La speranza è ancora che il costo complessivo sia ridotto. Notiamo che nonostante il nome, questo è un metodo esatto: a convergenza la soluzione è α (inexact si riferisce alla Jacobiana). Quello che cambia è la veloc ità di convergenza. pag. 21 Uno dei più utilizzati metodi inexact -Newton è quello di Broyden, che generalizza il metodo delle secanti (come quest’ultimo, anche Broyden è un metodo super -lineare (ordine>1), ma non del secondo ordine): pag. 22 Cap 3: Approssimazione di dati e funzioni Supponiamo di conosceren+1 coppie di dati (xi,yi),i= 0,… ,n Queste in generale possono costituire: 1. Misure yi di una certa quantità fisica effettuate sperimentalm ente in corrispondenza dei punti xi (nodi) . 2. Valori di una funzione f(xi) in corrispondenza dei punti xi. Il problema dell’approssimazione di dati (o di una funzione) consiste nel determinare una funzion e f(xi) che abbia delle buone proprietà di approssimazione (in un senso da speci ficare) dei dati (xi,yi),i= 0,… ,n che permettano di avere carattere predittivo anche non in corrispondenza dei nodi . Interpolazione Lagrangiana Approssimante di tipo interpolatorio Date n+1 coppie di dati (xi,yi),i= 0,… ,n , diciamo che un ’approssimante f̃(x) è di tipo interpolatorio se vale la seguente relazione: Un interpolatore è quindi una funzione che assume il valore dei dati in corrispondenza dei nodi xi. Polinomi caratteristici di Lagr ange Si introducono di seguitole seguen ti n+1 funzioni, associate agli n+1 dati (xi,yi),i= 0,… ,n che prendono il nome di polinomi caratteristici di Lagrange: Essi sono dati dal prodotto di n termini di primo grado, perciò sono dei polinomi di gra do n. se inoltre li valutiamo in corris pondenza di uno dei nodi otteniamo: Interpolatore Lagrangiano Siamo pronti per introdurre l’interpolatore Lagrangiano degli n+1 dati (xi,yi),i= 0,… ,n . esso è dato dalla seguente espressione: Quando stia mo approssimando i valori di una funzio ne (xi,f(xi)), si è soliti indicare l’interpolatore Lagrangiano con ∏ f(x) n . Per l’interpolatore Lagrangiano valgono le seguenti proprietà: 1. ∏ f(x) n è un interpolatore. Infatti, valutandolo per un generico nodo, si ottiene: pag. 23 2. �������(������) è un polinomio di grado n. Infatti, esso è dato dalla somma dei polinomi di grado n 3. ∏ f(x) n è l’unico polinomio di grado n interpolante gli n+1 dati (xi,yi),i= 0,… ,n. Infatti, supponiamo esista un altro interp olatore polinomiale di g rado n Ψn(x) che interp ola i dati (xi,yi),i= 0,… ,n Introduciamo quindi il seguente polinomio di grado n: Vale che: Allora D(x) ha n+1 zeri. Essendo un polinomio di grado n si deve avere D(x)≡ 0, da cui segue l’unicit à. Accuratezza dell’int erpolatore Lagrangiano Si consideri il caso di approssimazione dei valori di una funzione f(x). Si ha il seguente risultato: Si noti come la funzione ��������(������) sì annulli in corrispondenza dei nodi �������. Il suo valore è invece in generale diver so da zero lontano dai nodi, cioè per ������≠ �������. Al fine di ottenere un’informazione più compatta sull’errore, nel seguito deriviamo una stima del suo massimo valore. Nel caso di nodi equi -spaziati con vale la disuguaglianz a: [Infatti si ha ch e la funzione ammette massimo in prossimità degli estremi] Dall’espressione dell’errore nella proposizione e introducendo la stima, si ottiene che nel caso di nodi equi -spaziati vale la seguente stima dell’errore mas simo: pag. 24 Co nvergenza dell’interpola tore Lagrangiano Vorremmo che aumentando le informazioni a disposizione (cioè n+1) l’accuratezza di una funzione approssimante migliori. Nello specifico, per l’interpolatore Lagrangiano, vorremmo quindi che : Tornando alla stima: Notiamo che il denomina tore e il moltiplicante della frazione sono “buoni” nel senso che vanno a 0 per ������→ ∞. Infatti si ha che h=(b -a)/n con (a, b) l’intervallo in cui sono presi i nodi. Tuttavia, non abbiamo alcuna certezza a priori che il numerat ore della frazione sia limitato per ������→ ∞, cosa che garantirebbe la convergenza. Esistono casi per cui: ovvero la non convergenza dell’interpolatore Lagrangiano. Fenomeno di Runge La possibile mancata convergenza dell’interpolatore Lagrangiano prende il nome di fenomen o di Runge. In particolare, in presenza di questo fenomeno, la funzione errore ������� presenta delle oscillazioni ai nodi estremi che crescono con il crescere di n. Interpolazione Lagrangiana compo sita Per superare i limiti d ella interpolazione Lagr angiana, una possibile strategia consiste nell’introdurre un interpolatore continuo dato dalla unione di tanti interpolatori Lagrangiani di basso ordine, diciamo k≪ n. Tali interpolatori locali sono costruiti sugli intervalli disgiun ti Ij ognuno composto da k+1 nodi e di lunghezza H=kh. Tale interpolatore globale viene indicato con ΠkH(x). Quando vogliamo enfatizzare che esso è stato costruito per approssimare una fu nzione useremo la notazione ΠkHf(x). Accuratezza dell ’iterpolatore Lagrangian o composito In questo caso ci aspettiamo che le cose funzionino bene, cioè che quando aumentiamo il numero di informazioni n+1 che abbiamo a disposizione, l’accuratezza dell’inte rpolatore Lagrangiano composito migliori sempre . Infat ti, questa volta, quando aumenta n faremo aumentare il numero degli intervalli Ij su cui costruire gli interpolatori locali, senza variare il grado locale k che rimane costante . Se perciò si lavora con k piccolo (di solito1, 2 o 3), non si incorrerà più quando n cresce nel fen omeno di Runge . pag. 25 Introducendo l’errore puntuale EkHf(x)= f(x)− ΠkHf(x) e notando che esso è dato dal massimo degli errori degli interpolatori Lagrangiani di grado k su ogni Ij, otteniamo nel caso di nodi equispaziati: Avendo usato h=H/k Que sto dimostra la convergenza dell’interpolatore Lagrangiano composito. Scelta di H e k 1. Nel caso di funzione f(x) nota sull’intervallo [a,b] che genera i valori yi= f(xi),i= 0,… ,n. a. Si decide il valore di k da utilizzare (al massimo k=3 per evitare i l fenomeno di Runge) b. Si sceglie il valore dell’ampiezza degli intervalli H in modo da avere l’errore desiderato in base alla stima dell’errore c. Si partiziona [a,b] in intervalli di ampiezza H e su ognuno di essi si considerano k+1 n odi 2. Nel caso di dati prov enienti da misure, il numero di questi n+1 è fissato . In questi casi di solit o si opta per: a. Un’interpolazione composita lineare (k=1, H=(b -a)/n) b. Un’interpolazione composita quadratica (k=2, 2(b -a)/n) Interpolazione con ubicazione specifica dei nodi Per ovviare ai possibili problemi di stabilità dell’interpolatore Lagrangiano all ’aumentare del numero di informazioni disponibili n+1, si può considerare di ubicare i nodi in precise posizioni che garantiscono la stabilità al crescere di n. Questo permette di usare sempre con successo il polinomio interpolatore Lagrangiano di grado n. In generale, queste scelte particolari prevedono l’addensamento dei nodi in prossimità degli estremi del dominio di analisi, dove, come visto, insorgono le oscillazioni che porta no a instabilità. Nodi di Chebyshev Consideriamo dapprima il caso in cui il dominio di interesse dato da [ -1,1]. Al solito abbiamo n+1 dati, ubicati in corrispondenza delle seguenti ascisse: Questi nodi sono ottenuti come proi ezioni sull’asse x di pun ti sulla circonferenza unitaria individuati da settori circolari ottenuti con lo stesso angolo π/n. Si nota come essi siano effettivamente più addensati verso le estremità . Interpolatore Lagrangiano sui nodi di Chebyshev Consideri amo ora l’interpolatore Lagrangiano di grado n costruito sui nodi di Chebyshev. Questo viene costruito esattamente come fatto in precedenza, a patto di prendere ������̂� come nodi su cui costruire gli n polinomi caratteristici di Lagrange . pag. 26 Esso viene indic ato usualmente con ΠnC(x). Indicheremo l’interpolatore con ΠnCf(x) quando applicato a dati ottenuti dalla valutazione di una funzione f, ovvero yi= f(xi). Facciamo ora un piccolo confronto dei polinomi caratteristici. Consideriamo 11 nodi (n= 10) e i corrispondenti polinomi caratteristici di Lag range. Nodi di Chebyshev Nodi equispaziati Si notano le assenze di oscillazioni nel caso Chebyshev in prossim ità degli estremi, dovute ad un infittimento dei nodi in quelle regioni . Convergenza dell’interpolazione sui nodi di Chebyshev Teorema: supponiamo che la funzione f(x) ammetta derivata continua fino all’ordine s+1 compreso, ovvero f∈Cs+1([−1,1]). Allora, si ha il seguente ris ultato di convergenza: (per un’opportuna costante C) Il precedente risultato ci dice che: i. Se s≥ 1 allora si ha sempre convergenza ii. La velocità di convergenza è tanto più alta tanto più è alto s iii. Se l’ipote si del Teorema vale per ogni s > 0, tale velocit à è esponenziale, che è la stima migliore che si può ottenere, essendo più veloce di qualsiasi 1/n5: Nodi di Chebyshev su un intervallo generale Nel caso in cui il dominio di interesse sia [a,b], è suff iciente mappare i nodi ������̂� da [ -1,1] ad [a ,b], ottenendo i nodi. A questo fine introduciamo la seguente mappa : Che agisce dal dominio [-1,1] nella variabile ������̂ al dominio corrente [a,b] . I nodi di Chebyshev su un intervallo generico [a,b] sono quindi dati da: pag. 27 I pol inomi caratteristici di Lagrange seguono al solito modo: Seguono in manier a analoga al caso di dominio [-1,1] la definizione dei polinomi caratteristici di Lagrange, dell’interpolatore Lagrangiano e delle sue proprietà di convergenza . Interpolazione trigonometrica Consid eriamo ora il caso di una funzione periodica, ad esempio di periodo 2π e con f(0) = f(2 π) , di cui si conoscano gli n+1 valori f(xj) nei nodi equispaziati : In questo caso ci aspettiamo che l’interpolatore Lagrangiano non app rossimi bene la funz ion e f lontano dai nodi, essendo un polinomio di grado n, che quindi non è in generale periodico . Si vuole quindi costruire un interpolatore periodico f̃(x). Interpolatore trigonometrico Supponiamo per semplicità che n sia pari e poni amo M=n/2 . L’idea è di scrivere l’interpolatore come combinazione lineare di seni e coseni: Per opportuni coefficienti complessi incogniti ak,per k= 0,… ,M e bk