Introduzione al problema Sia assegnata una funzione f(x) di cui sono noti i valori f(x1), f(x2),…,f(xN) nei punti x1, x2, …, xN per semplicità supponiamo x1<x2<…<xN Si vuole calcolare il valore di f(x) in un punto x ≠ x1, x2, …, xN se x è compreso nell’intervallo [x1,xN] si parla di interpolazione se x non è compreso nell’intervallo [x1,xN] si parla di estrapolazione Con l’interpolazione (o l’estrapolazione) si costruisce un modello della funzione f(x) valido all’interno (o all’esterno) dell’intervallo [x1,xN] Un esempio (1) Come esempio consideriamo la seguente funzione: f ( x) 3x 2 1 4 ln x 1 2 Tale funzione è definita ovunque, eccetto che in x=π, dove si ha: lim f ( x) x Supponiamo di voler interpolare la funzione f(x) in x=π partendo dai valori di f in alcuni punti abbastanza “vicini” a x=π Un esempio (2) Un esempio (3) Anche utilizzando un programma di grafica non è facile mettere in evidenza la divergenza di f(x) in x=π Qualsiasi interpolazione basata sui valori di f(x) in punti “vicini” a x=π fornirà un risultato errato! Nell’effettuare un’interpolazione (o un’estrapolazione) solitamente si assume che la funzione in esame sia abbastanza “regolare” (“smooth”) la condizione di “regolarità” (“smoothness”) in sostanza equivale a richiedere che la f(x) sia continua e che esistano e siano continue le sue derivate f’(x), f’’(x), ..., f(N)(x) fino all’ordine N tanto più grande è N, tanto più la funzione f(x) è regolare (“smooth”) Strategie di interpolazione Interpolazione mediante fit: si sceglie la funzione interpolatrice g(x), che di solito dipende da un insieme di parametri si effettua il fit degli N punti conosciuti (xi, yi=f(xi)) con la funzione g(x) si calcola il valore della funzione interpolatrice g(x) nel punto x in cui si vuole effettuare l’interpolazione e si assume f(x)=g(x) Questo metodo solitamente è poco efficiente dal punto di vista del calcolo ed è soggetto ad errori di arrotondamento Interpolazione per approssimazioni successive: si parte dal valore di f(xi) nel punto xi più vicino a x si aggiungono a f(xi) una serie di termini correttivi, sempre più piccoli man mano che si va avanti, che tengono conto dell’informazione contenuta negli altri valori f(xj) forniti dal problema l’ultimo termine correttivo fornisce una stima dell’errore su f(x) Questa procedura richiede tipicamente un numero di operazioni ~O(N2) In generale la f(x) (e le derivate di ordine superiore) non è continua Ordine di interpolazione Ordine dell’interpolazione = N-1 (N = numero di punti usati nell’interpolazione) Non sempre aumentare l’ordine dell’interpolazione ne migliora la precisione In particolare, se i punti aggiunti sono lontani dal punto di interesse x, la funzione interpolante rischia di essere caratterizzata da forti oscillazioni, che non riproducono il comportamento della funzione vera f(x) tipicamente 3-4 punti sono sufficienti per una buona interpolazione Formula di Lagrange Dati N punti nel piano cartesiano, esiste un unico polinomio di grado N-1 che passa per tutti i punti assegnati Siano (x1,y1=f(x1)), (x2,y2=f(x2)), ..., (xN,yN=f(xN)) gli N punti assegnati. Il polinomio P(x) di grado N-1 che passa per gli N punti si ottiene con la formula di Lagrange: x x2 x x3 x xN y x x1 x x3 x xN y x1 x2 x1 x3 x1 xN 1 x2 x1 x2 x3 x2 xN 2 x x1 x x2 x xN 1 y xN x1 xN x2 xN xN 1 N P( x) Il termine i-esimo della somma presente nella formula di Lagrange è così composto: a numeratore la yi a numeratore il prodotto di N-1 termini del tipo (x-xj) con j≠i a denominatore il prodotto di N-1 termini del tipo (xi-xj) con j≠i Si noti che se x=xi il termine i-esimo della somma vale yi, mentre gli altri sono nulli Algoritmo di Neville (1) Indichiamo con P1(x) l’unico polinomio di grado zero che passa per il punto (x1,y1), con P2(x) l’unico polinomio di grado zero che passa per il punto (x2,y2) e così via. Nel punto x avremo così N valori: P1(x)=y1 , P2(x)=y2 , ... , PN(x)=yN Indichiamo quindi con P12(x) l’unico polinomio di primo grado che passa per i punti (x1,y1) e (x2,y2), con P23(x) l’unico polinomio di primo grado che passa per i punti (x2,y2) e (x3,y3), e così via. Nel punto x avremo N-1 valori: P12(x), P23(x), ..., P(N-1)N(x) In maniera analoga possiamo continuare a definire i polinomi di grado più alto, fino all’unico polinomio di grado N-1 che passa per tutti i punti (x1,y1),(x2,y2),..., (xN,yN), il cui valore in x sarà: P12...N(x) Algoritmo di Neville (2) Con tutti i valori così calcolati si può costruire una tabella. Per esempio, se N=4 la tabella avrà la forma seguente: y1=P1 y2=P2 y3=P3 y4=P4 P12 P123 P23 P34 P1234 P234 L’algoritmo di Neville permette di riempire la tabella sfruttando le relazioni ricorsive tra i polinomi: Pi (i 1)(i m ) x xi m Pi (i 1)(i m1) xi x P(i 1)(i 2)(i m) xi xi m Dimostrazione della formula di Neville (1) Partiamo dalla formula di Lagrange scritta in forma compatta: im im j i k i k j Pi ( i 1)( i m ) y j Pi ( i 1)( i m 1) P( i 1)( i 2 )( i m ) j i x j xk i m 1 i m 1 x xk x xk x yj k i k j im im x xk y x j i 1 xk j j k i 1 k j j xk Dimostrazione della formula di Neville (2) Sostituiamo quindi nel secondo membro della formula di Neville le espressioni di Pi(i+1)...(i+m-1) e P(i+1)(i+2)...(i+m): Pi ( i 1)( i m ) x xi m xi xi m i m 1 x xi xi m xi im j i i m 1 yj x xk x k i k j im j k i 1 k j xk x xk y x j i 1 j j xk Dobbiamo a questo punto verificare che vale effettivamente l’uguaglianza tra i due membri A tal fine isoliamo dalla prima sommatoria il termine con l’indice j=i e dalla seconda sommatoria il termine con l’indice j=i+m Dimostrazione della formula di Neville (3) Si ha quindi: Pi ( i 1)( i m ) x xi m i m 1 x xk yi xi xi m k i 1 xi xk x xi m xi xi m i m 1 x xi xi m xi i m 1 j i 1 j i 1 i m 1 yj x xk x k i k j im yj k i 1 k j i m 1 1 j xk 3 x xk x j xk x xi x xk yi m xi m xi k i 1 xi m xk 2 Dimostrazione della formula di Neville (4) I termini (1) e (2) nell’equazione precedente possono essere riscritti come segue: (1) (2) x xk yi k i 1 xi xk i m 1 x xk y im im x k i j xk x xk yi k i xi xk k i im x xk y im im x k i k i m j xk In sostanza i termini (1) e (2) sono il primo e l’ultimo termine nella sommatoria che compare nella formula di Lagrange per Pi(i+1)...(i+m) Dimostrazione della formula di Neville (5) Esaminiamo ora il termine (3): x xk x xi x xi m x xi m x xi yj x j xi xi xi m x j xi m xi m xi j i 1 k i 1 x j xk k j i m 1 (3) i m 1 Mettendo in evidenza i termini comuni si ha: x xi x xi m xi xi m 1 1 x x x x j i j im x xi x xi m x j xi m x j xi xi xi m x j xi x j xi m x xi x xi m x j xi x j xi m Dimostrazione della formula di Neville (6) Il termine (3) pertanto diventa: (3) i m 1 i m 1 j i 1 yj x xk x xi x xi m x k i 1 k j j xk x j xi x j xi m im i m 1 x xk y x j i 1 j k i k j Raggruppando i termini (1), (2) e (3) abbiamo quindi trovato tutti i termini nella sommatoria che compare nello sviluppo di Lagrange di Pi(i+1)...(i+m) j xk Algoritmo di Neville (3) Nel costruire la tabella dei valori con l’algoritmo di Neville possiamo definire le differenze tra “genitori” e “figli” nel modo seguente: Cm,i Pi (i 1)(i m) Pi (i 1)(i m1) Dm,i Pi (i 1)(i m) P(i 1)(i 2)(i m) I termini Cm,i e Dm,i rappresentano le correzioni che occorre apportare per aumentare di una unità l’ordine di interpolazione Valgono le seguenti relazioni ricorsive: Dm 1,i Cm 1,i xi m1 x Cm,i 1 Dm,i xi xi m 1 xi x Cm,i 1 Dm,i xi xi m 1 Dimostrazione delle formule ricorsive (1) La dimostrazione sfrutta la formula di Neville: Cm 1,i Pi (i 1)(i m 1) Pi ( i 1)(i m ) x xi m1 Pi (i 1)(i m) xi x P(i 1)(i 2)(i m1) xi xi m 1 xi x x xi m1 P xi xi m 1 i ( i 1)( i m ) xi x P(i 1)(i 2)(i m1) Pi (i 1)(i m) xi xi m 1 xi x P(i 1)(i 2)(i m1) P(i 1)(i 2)(i m) P(i 1)(i 2)(i m) Pi (i 1)(i m) xi x Cm,i 1 Dm,i xi xi m 1 xi xi m 1 Dimostrazione delle formule ricorsive (2) La dimostrazione sfrutta la formula di Neville: Dm 1,i Pi (i 1)(i m 1) P(i 1)( i 2)( i m 1) x xi m1 Pi (i 1)(i m) xi x P(i 1)(i 2)(i m1) xi xi m 1 xi x x xi m1 P xi xi m 1 ( i 1)( i 2 )( i m 1) x xi m1 Pi (i 1)(i m) P(i 1)(i 2)(i m1) xi xi m 1 xi m1 xi P(i 1)(i 2)(i m1) P(i 1)(i 2)(i m) P(i 1)(i 2)(i m) Pi (i 1)(i m) xi m1 x Cm,i 1 Dm,i xi xi m 1 xi xi m 1 Algoritmo di Neville (4) Ordine di interpolazione 1 0 y1=P1 y3=P3 3 C1,1 D1,1 y2=P2 2 P12 C2,1 D2,1 C1,2 D1,2 P23 C1,3 D1,3 P34 C2,2 P123 P234 C3,1 P1234 D3,1 D2,2 y4=P4 Partendo da uno dei valori Pi, e aggiungendo le correzioni date dai coefficienti Cm,i e Dm,i , seguendo la mappa si può arrivare all’ordine di interpolazione desiderato Un esempio di applicazione dell’algoritmo di Neville x y -2.0 -1.5 -1.0 0.2 0.0 1.0 1.0 1.4 2.0 0.2 Interpolazione con funzioni razionali Alcune funzioni possono non essere ben approssimate da polinomi In questi casi si può utilizzare una funzione razionale, cioè un rapporto tra polinomi Le funzioni razionali costituiscono una buona approssimazione quando la funzione f(x) che si vuole interpolare presenta dei poli (limite infinito) Indichiamo con Ri(i+1)...(i+m) una generica funzione razionale che passa per i gli m+1 punti (xi,yi),(xi+1,yi+1),...,(xi+m,yi+m): Ri (i 1)( i m ) ( x) P ( x) Q ( x) p0 p1 x p x q0 q1 x q x Il grado del polinomio a numeratore e di quello a denominatore, μ e ν, sono arbitrari, ma si deve tener conto di un vincolo: nell’espressione di Ri(i+1)...(i+m) ci sono μ+ν+1 incognite (q0 è arbitraria) Ri(i+1)...(i+m) deve passare per m+1 punti assegnati ne consegue che deve essere μ+ν+1 = m+1 Algoritmo di Burlisch e Stoer (1) Si tratta di un algoritmo simile a quello di Neville per costruire delle funzioni razionali che interpolino i dati Le funzioni razionali prodotte dall’algoritmo di Burlisch e Stoer hanno: numeratore e denominatore di grado uguale se m+1 è dispari grado del denominatore maggiore di una unità rispetto al grado del numeratore se m+1 è pari Come nell’algoritmo di Neville, le funzioni razionali Ri(i+1)...(i+m) sono generate a partire da una relazione ricorsiva in questo caso la relazione ricorsiva non deriva da alcuna formula (mentre nell’algoritmo di Neville deriva dalla formula di Lagrange) Anche in questo caso si costruisce una tabella che permette di passare da un ordine di interpolazione al successivo Algoritmo di Burlisch e Stoer (2) Come nell’algoritmo di Neville si parte dalle funzioni costanti R1(x)=y1, R2(x)=y2, ... , RN(x)=yN (Ri=yi) si pone inoltre R=[Ri(i+1)...(i+m) con m=-1]=0 la relazione ricorsiva che definisce le Ri(i+1)...(i+m) è: Ri (i 1)(i m ) R(i 1)( i 2)(i m ) R(i 1)( i 2)(i m ) Ri (i 1)(i m 1) R(i 1)( i 2)(i m ) Ri (i 1)(i m 1) x xi 1 1 x xi m R(i 1)( i 2)(i m ) R(i 1)( i 2)(i m 1) Analogamente all’algoritmo di Neville, si pone: Cm,i Ri (i 1)(i m) Ri (i 1)(i m1) Dm,i Ri (i 1)(i m) R(i 1)(i 2)(i m) e si verifica che: Cm 1,i Dm1,i Cm,i 1 Dm,i Dimostrazione della relazione tra i coefficienti C e D La dimostrazione segue dalle definizioni: Cm 1,i Dm 1,i R i ( i 1)( i m 1) Ri (i 1)(i m ) Ri ( i 1)(i m 1) R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m 1) Ri ( i 1)(i m ) R( i 1)( i 2 )( i m 1) R(i 1)( i 2 )(i m ) R(i 1)( i 2 )(i m ) Ri ( i 1)(i m ) Cm ,i 1 Dm ,i Algoritmo di Burlisch e Stoer (3) La relazione precedente permette di dimostrare le seguenti formule ricorsive: Dm 1,i Cm 1,i Cm,i 1 Cm ,i 1 Dm,i x xi Dm,i Cm,i 1 x xi m 1 x xi Dm,i Cm,i 1 Dm,i x xi m 1 x xi Dm,i Cm,i 1 x xi m 1 Dimostrazione delle formule ricorsive (1) Dm 1,i Ri ( i 1)( i m 1) R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m 1) Ri ( i 1)( i m ) R( i 1)( i 2 )( i m 1) Ri ( i 1)( i m ) x xi 1 R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m ) x xi m 1 1 R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m 1) R( i 1)( i 2 )( i m ) R( i 1)( i 2 )( i m ) Ri ( i 1)( i m ) R R( i 1)( i 2 )( i m ) R( i 1)( i 2 )( i m ) Ri ( i 1)( i m ) x xi 1 ( i 1)( i 2 )( i m 1) Cm ,i 1 x xi m 1 Cm ,i 1 Dm ,i x xi Cm ,i 1 Dm ,i 1 Cm ,i 1 x xi m 1 Cm ,i 1 Cm ,i 1 Dm ,i x xi Dm ,i Cm ,i 1 x xi m 1 1 C m ,i 1 Dm ,i x xi Cm ,i 1 Cm ,i 1 Dm ,i Cm ,i 1 x xi m 1 1 1 Dimostrazione delle formule ricorsive (2) Cm 1,i Dm 1,i Cm ,i 1 Dm ,i Cm ,i 1 Cm ,i 1 Dm ,i x xi Dm ,i Cm ,i 1 x xi m 1 Cm ,i 1 Dm ,i x xi Dm ,i Cm ,i 1 Cm ,i 1 x xi m 1 Cm ,i 1 Dm ,i x xi Dm ,i Cm ,i 1 x xi m 1 x xi Dm ,i x xi m 1 Cm ,i 1 Dm ,i x xi Dm ,i Cm ,i 1 x xi m 1 Algoritmo di Burlisch e Stoer (4) Ordine di interpolazione 1 0 y1=R1 y3=R3 3 C1,1 D1,1 y2=R2 2 R12 C2,1 D2,1 C1,2 D1,2 R23 C1,3 D1,3 R34 C2,2 R123 R234 C3,1 R1234 D3,1 D2,2 y4=R4 Come nell’algoritmo di Neville, si parte da uno dei valori Ri e, aggiungendo le correzioni date dai coefficienti Cm,i e Dm,i , seguendo la mappa si può arrivare all’ordine di interpolazione desiderato Un esempio di applicazione dell’algoritmo di Burlisch e Stoer x y -2.0 -1.5 -1.0 0.2 0.0 1.0 1.0 1.4 2.0 0.2 Il punto di partenza: la formula di Lagrange Sia data una funzione y(x) tabulata in N punti Indichiamo con (xi, yi) i punti tabulati Applicando la formula di Lagrange nell’intervallo [xj,xj+1], la funzione y(x) può essere approssimata linearmente nel modo seguente: y ( x) y j x x j 1 x j x j 1 y j 1 x xj x j 1 x j dove: A B x x j 1 x j x j 1 x xj x j 1 x j 1 A Ay j By j 1 La funzione “spline” cubica Indichiamo ora con y’’i i valori delle derivate seconde della y(x) nei punti xi e cerchiamo di migliorare l’approssimazione lineare tra xj e xj+1 ponendo: y Ay j By j 1 Cyj Dyj1 con: A B C D x j 1 x x j 1 x j 1 A 1 3 2 A A x j 1 x j 6 1 3 2 B B x j 1 x j 6 La funzione y(x) così costruita è un polinomio di terzo grado in x la x compare in A (e in B) al primo grado, e i termini di grado massimo in A e B sono di terzo grado Derivata prima della funzione “spline” Calcoliamo preliminarmente le derivate di A,B,C,D: A x j 1 x x j 1 x j B 1 A dA 1 dx x j 1 x j dB dA 1 dx dx x j 1 x j 1 3 dC 1 dA 3 A2 1 2 2 2 x j 1 x j C A Ax j 1 x j x j 1 x j 3 A 1 6 dx 6 dx 6 1 3 dD 1 dB 3B 2 1 2 2 2 x j 1 x j D B B x j 1 x j x j 1 x j 3B 1 6 dx 6 dx 6 Ne segue: dy 1 1 3 A2 1 3B 2 1 x j 1 x j yj1 x j 1 x j yj y j 1 yj dx x j 1 x j x j 1 x j 6 6 y j 1 y j 3 A2 1 3B 2 1 x j 1 x j yj x j 1 x j yj1 x j 1 x j 6 6 Derivata seconda della funzione “spline” Derivando ulteriormente si ha: d2y 6 A 1 6B 1 x j 1 x j yj1 Ayj Byj1 x j 1 x j y j 2 dx 6 x j 1 x j 6 x j 1 x j se x=xj si ha A=1 e B=0 e quindi y’’(xj)=y’’j se x=xj+1 si ha A=0 e B=1 e quindi y’’(xj+1)=y’’j+1 Questo risultato dimostra che i coefficienti A,B,C,D sono definiti in maniera corretta I valori delle y’’i in genere non sono noti, e potrebbero, in principio, essere assegnati arbitrariamente In generale, però, si preferisce assegnare i valori delle y’’i in modo da garantire la continuità della derivata prima della funzione “spline” Continuità della derivata prima (1) Derivata prima della “spline” tra xj e xj+1: dy y j 1 y j 3 A2 1 3B 2 1 x j 1 x j yj x j 1 x j yj1 dx x j 1 x j 6 6 in x=xj è A=1 e B=0 per cui: y j 1 y j 1 1 y( x j ) x j 1 x j yj x j 1 x j yj1 x j 1 x j 3 6 Derivata prima della “spline” tra xj-1 e xj: dy y j y j 1 3 A2 1 3B 2 1 x j x j 1 yj1 x j x j 1 yj dx x j x j 1 6 6 i coefficienti A e B sono diversi dai precedenti! in x=xj stavolta è A=0 e B=1 per cui: y( x j ) y j y j 1 x j x j 1 1 x j x j 1 yj1 1 x j x j 1 yj 6 3 Continuità della derivata prima (2) Imponendo che le due espressioni della derivata prima in xj siano uguali si ha: y j 1 y j x j 1 x j y j y j 1 x j x j 1 1 x j 1 x j yj 1 x j 1 x j yj1 3 6 1 x j x j 1 yj1 1 x j x j 1 yj 6 3 y j 1 y j y j y j 1 1 1 1 x j x j 1 yj1 x j 1 x j 1 yj x j 1 x j yj1 6 3 6 x j 1 x j x j x j 1 Si ottengono così N-2 equazioni nelle N incognite y’’i per cui occorre fissare due di queste incognite per risolvere il sistema si può porre y’’1=0 e y’’N=0 in questo caso si parla di “spline cubica naturale” alternativamente, si possono determinare y’’1 e y’’N se sono noti i valori della derivata prima negli estremi dell’intervallo Alcuni esempi di spline cubiche x y -2.0 -1.5 -1.0 0.2 0.0 1.0 1.0 1.4 2.0 0.2 Interpolazione multidimensionale Il problema generale è quello di valutare una funzione y(x1,x2,...xn)=y(x) a partire da una griglia di N valori assegnati y1=y(x1), y2=y(x2),..., yN=y(xN) Per semplicità studieremo il problema nel solo caso bidimensionale Assumeremo inoltre che siano assegnati i valori della funzione da interpolare in corrispondenza dei vertici di una griglia rettangolare Interpolazione in due dimensioni y yN y2 y1 x1 x2 x3 xM x Supponiamo che siano assegnati i valori della funzione z=f(x,y) in corrispondenza degli MN vertici di una griglia indichiamo con x1,x2,...,xM e con y1,y2,...,yN le ascisse e le ordinate dei vertici del reticolo indichiamo con zjk=z(xj,yk) i valori della funzione in corrispondenza dei vertici del reticolo Interpolazione bilineare Dato il punto (x,y), (1-u)(yk+1-yk) indichiamo con j e k gli indici per cui si ha xj<x<xj+1 e yk<y<yk+1 Poniamo: yk+1 y yk u(yk+1-yk) xj x xj+1 t(xj+1-xj) t x xj x j 1 x j y yk u yk 1 yk (1-t)(xj+1-xj) Interpolazione bilineare: z ( x, y ) (1 t )(1 u ) z j ,k t (1 u ) z j 1,k tuz j 1,k 1 (1 t )uz j ,k 1