Introduzione al problemaSia 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:
Tale funzione è definita ovunque, eccetto che in x=π, dove si ha:
Supponiamo di voler interpolare la funzione f(x) in x=π partendo dai valori di f in alcuni punti abbastanza “vicini” a x=π
1ln1
3)( 2
42 xxxf
)(lim xfx
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 interpolazioneInterpolazione 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 xsi 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 interpolazioneOrdine dell’interpolazione = N-1 (N =
numero di punti usati nell’interpolazione)Non sempre aumentare l’ordine
dell’interpolazione ne migliora la precisioneIn 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 LagrangeDati N punti nel piano cartesiano, esiste un unico
polinomio di grado N-1 che passa per tutti i punti assegnatiSiano (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:
Il termine i-esimo della somma presente nella formula di Lagrange è così composto:a numeratore la yia numeratore il prodotto di N-1 termini del tipo (x-xj) con j≠ia denominatore il prodotto di N-1 termini del tipo (xi-xj) con
j≠iSi noti che se x=xi il termine i-esimo della somma vale yi,
mentre gli altri sono nulli
N
NNNN
N
N
N
N
N
yxxxxxx
xxxxxx
yxxxxxx
xxxxxxy
xxxxxx
xxxxxxxP
121
121
223212
311
13121
32)(
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:
L’algoritmo di Neville permette di riempire la tabella sfruttando le relazioni ricorsive tra i polinomi:
y1=P1P12
P123
P1234
y2=P2P23
y3=P3P234P34
y4=P4
mii
miiiimiiimimiii xx
PxxPxxP
)()2)(1()1()1(
)()1(
Dimostrazione della formula di Neville (1)Partiamo dalla formula di Lagrange scritta in
forma compatta:
mi
ij
mi
jkik kj
kjmiii
mi
ij
mi
jkik kj
kjmiii
mi
ij
mi
jkik kj
kjmiii
xx
xxyP
xx
xxyP
xx
xxyP
1 1)()2)(1(
1 1
)1()1(
)()1(
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):
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
mi
ij
mi
jkik kj
kj
imi
i
mi
ij
mi
jkik kj
kj
mii
mimiii
xx
xxy
xx
xx
xx
xxy
xx
xxP
1 1
1 1
)()1(
Dimostrazione della formula di Neville (3)Si ha quindi:
1
1
1
1 1
1
1
1
1
1)()1(
mi
ik kmi
kmi
imi
i
mi
ij
mi
jkik kj
kj
imi
i
mi
ij
mi
jkik kj
kj
mii
mi
mi
ik ki
ki
mii
mimiii
xx
xxy
xx
xx
xx
xxy
xx
xx
xx
xxy
xx
xx
xx
xxy
xx
xxP
1
2
3
Dimostrazione della formula di Neville (4)I termini (1) e (2) nell’equazione precedente
possono essere riscritti come segue:
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)
mi
mikik kj
kmi
mi
ik kj
kmi
mi
ikik ki
ki
mi
ik ki
ki
xx
xxy
xx
xxy
xx
xxy
xx
xxy
1
1
)2(
)1(
Dimostrazione della formula di Neville (5)Esaminiamo ora il termine (3):
Mettendo in evidenza i termini comuni si ha:
1
1
1
1
)3(mi
ij imi
i
mij
mi
mii
mi
ij
imi
jkik kj
kj xx
xx
xx
xx
xx
xx
xx
xx
xx
xxy
mijij
mii
mijij
ijmij
mii
mii
mijijmii
mii
xxxx
xxxx
xxxx
xxxx
xx
xxxx
xxxxxx
xxxx
11
Dimostrazione della formula di Neville (6)Il termine (3) pertanto diventa:
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)
1
1
1
1
1
1
)3(mi
ij
mi
jkik kj
kj
mi
ij mijij
miimi
jkik kj
kj xx
xxy
xxxx
xxxx
xx
xxy
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:
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:
)()2)(1()()1(,
)1()1()()1(,
miiimiiiim
miiimiiiim
PPD
PPC
1
,1,,1
1
,1,1,1
mii
imimiim
mii
imimmiim
xx
DCxxC
xx
DCxxD
Dimostrazione delle formule ricorsive (1)La dimostrazione sfrutta la formula di
Neville:
1
,1,
1
)()1()()2)(1()()2)(1()1()2)(1(
1
)()1()1()2)(1(
)()1(1
1
1
)1()2)(1()()1(1
)()1()1()1(,1
mii
imimi
mii
miiimiiimiiimiiii
mii
miiimiiii
miiimii
mii
mii
miiiimiiimi
miiimiiiim
xx
DCxx
xx
PPPPxx
xx
PPxx
Pxx
xxxx
xx
PxxPxx
PPC
Dimostrazione delle formule ricorsive (2)La dimostrazione sfrutta la formula di
Neville:
1
,1,1
1
)()1()()2)(1()()2)(1()1()2)(1(1
1
)1()2)(1()()1(1
)1()2)(1(1
1
1
)1()2)(1()()1(1
)1()2)(1()1()1(,1
mii
imimmi
mii
miiimiiimiiimiiiimi
mii
miiimiiimi
miiimii
mii
mii
miiiimiiimi
miiimiiiim
xx
DCxx
xx
PPPPxx
xx
PPxx
Pxx
xxxx
xx
PxxPxx
PPD
Algoritmo di Neville (4)
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
y1=P1
y2=P2
y3=P3
y4=P4
P12
P23
P34
P123
P234
P1234
C1,1
D1,1
C1,2
D1,2
C1,3
D1,3
C2,1
D2,1
C2,2
D2,2
C3,1
D3,1
Ordine di interpolazione0 1 2 3
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 razionaliAlcune funzioni possono non essere ben approssimate da
polinomiIn questi casi si può utilizzare una funzione razionale, cioè un
rapporto tra polinomiLe 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):
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 assegnatine consegue che deve essere μ+ν+1 = m+1
xqxqq
xpxpp
xQ
xPxR miii
10
10)()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 datiLe funzioni razionali prodotte dall’algoritmo di
Burlisch e Stoer hanno: numeratore e denominatore di grado uguale se m+1 è
disparigrado del denominatore maggiore di una unità rispetto
al grado del numeratore se m+1 è pariCome nell’algoritmo di Neville, le funzioni razionali
Ri(i+1)...(i+m) sono generate a partire da una relazione ricorsivain 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) è:
Analogamente all’algoritmo di Neville, si pone:
e si verifica che:
11)1()2)(1()()2)(1(
)1()1()()2)(1(
)1()1()()2)(1()()2)(1()()1(
miiimiii
miiimiii
mi
i
miiimiiimiiimiii
RR
RR
xxxx
RRRR
)()2)(1()()1(,
)1()1()()1(,
miiimiiiim
miiimiiiim
RRD
RRC
imimimim DCDC ,1,,1,1
Dimostrazione della relazione tra i coefficienti C e DLa dimostrazione segue dalle definizioni:
imim
miiimiiimiiimiii
miiimiii
miiimiiimiiimiii
imim
DC
RRRR
RR
RRRR
DC
,1,
)()1()()2)(1()()2)(1()1()2)(1(
)()1()1()2)(1(
)1()2)(1()1()1()()1()1()1(
,1,1
Algoritmo di Burlisch e Stoer (3)La relazione precedente permette di
dimostrare le seguenti formule ricorsive:
1,,1
,1,,1
,1
1,,1
,1,1,,1
imimmi
i
imimimmi
i
im
imimmi
i
imimimim
CDxxxx
DCDxxxx
C
CDxxxx
DCCD
Dimostrazione delle formule ricorsive (1)
1,,1
,1,1,
1,
,1,1,
1
,1,
1,
,1,
1
,1,
1,
)()1()()2)(1()()2)(1()1()2)(1(
1
)()1()()2)(1()()2)(1()1()2)(1(
)1()2)(1(
)()2)(1()1()2)(1(
)()1()1()2)(1(
1
)()1()1()2)(1()1()2)(1(
)1()2)(1()1()1(,1
111
11
11
imimmi
i
imimim
im
imimim
mi
i
imim
im
imim
mi
i
imim
im
miiimiiimiiimiii
mi
i
miiimiiimiiimiii
miii
miiimiii
miiimiii
mi
i
miiimiiimiii
miiimiiiim
CDxxxx
DCC
C
DCC
xxxx
DC
C
DC
xxxx
DC
C
RRRR
xxxx
RRRR
R
RR
RR
xxxx
RRR
RRD
Dimostrazione delle formule ricorsive (2)
1,,
1
,1
,1,
1,,1
1,,1
1,
,1,
,1,
1,,1
,1,1,
,1,,1,1
imimmi
i
immi
i
imim
imimmi
i
imimmi
iim
imim
imim
imimmi
i
imimim
imimimim
CDxx
xx
Dxx
xx
DC
CDxx
xx
CDxx
xxC
DC
DC
CDxx
xx
DCC
DCDC
Algoritmo di Burlisch e Stoer (4)
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
y1=R1
y2=R2
y3=R3
y4=R4
R12
R23
R34
R123
R234
R1234
C1,1
D1,1
C1,2
D1,2
C1,3
D1,3
C2,1
D2,1
C2,2
D2,2
C3,1
D3,1
Ordine di interpolazione0 1 2 3
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 LagrangeSia data una funzione y(x) tabulata in N puntiIndichiamo con (xi, yi) i punti tabulatiApplicando la formula di Lagrange nell’intervallo
[xj,xj+1], la funzione y(x) può essere approssimata linearmente nel modo seguente:
dove:
11
11
1)(
jjjj
jj
jj
jj ByAy
xx
xxy
xx
xxyxy
Axx
xxB
xx
xxA
jj
j
jj
j
11
1
1
La funzione “spline” cubicaIndichiamo 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:
con:
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
11 jjjj yDyCByAyy
21
3
21
3
1
1
6
16
11
jj
jj
jj
j
xxBBD
xxAAC
ABxx
xxA
Derivata prima della funzione “spline”Calcoliamo preliminarmente le derivate di
A,B,C,D:
Ne segue:
jjjjjj
jjjjjj
jj
jjjj
j
xxB
dx
dBBxx
dx
dDxxBBD
xxA
dx
dAAxx
dx
dCxxAAC
xxdx
dA
dx
dBAB
xxdx
dA
xx
xxA
1
222
12
13
1
222
12
13
1
11
1
6
1313
6
1
6
1
6
1313
6
1
6
1
11
1
11
2
1
2
1
1
1
2
11
2
11
1
6
13
6
13
6
13
6
1311
jjjjjjjj
jj
jjjjjjjj
jjj
j
yxxB
yxxA
xx
yy
xxB
yxxA
yxx
yxx
ydx
dy
Derivata seconda della funzione “spline”Derivando ulteriormente si ha:
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”
1111
11
2
2 1
6
61
6
6
jjjjjjj
jjjjj
yByAyxxxx
Byxx
xx
A
dx
yd
Continuità della derivata prima (1)Derivata prima della “spline” tra xj e xj+1:
in x=xj è A=1 e B=0 per cui:
Derivata prima della “spline” tra xj-1 e xj:
i coefficienti A e B sono diversi dai precedenti!in x=xj stavolta è A=0 e B=1 per cui:
11
2
1
2
1
1
6
13
6
13
jjjjjjjj
jj yxxB
yxxA
xx
yy
dx
dy
1111
1
6
1
3
1)(
jjjjjjjj
jjj yxxyxx
xx
yyxy
jjjjjjjj
jj yxxB
yxxA
xx
yy
dx
dy
1
2
11
2
1
1
6
13
6
13
jjjjjjjj
jjj yxxyxx
xx
yyxy
111
1
1
3
1
6
1)(
Continuità della derivata prima (2)Imponendo che le due espressioni della derivata
prima in xj siano uguali si ha:
Si ottengono così N-2 equazioni nelle N incognite y’’i per cui occorre fissare due di queste incogniteper 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
1
1
1
1111111
1111
1
1111
1
6
1
3
1
6
1
3
1
6
1
6
1
3
1
jj
jj
jj
jjjjjjjjjjj
jjjjjjjj
jj
jjjjjjjj
jj
xx
yy
xx
yyyxxyxxyxx
yxxyxxxx
yy
yxxyxxxx
yy
Alcuni esempi di spline cubichex y
-2.0
-1.5
-1.0
0.2
0.0 1.0
1.0 1.4
2.0 0.2
Interpolazione multidimensionaleIl 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
Supponiamo che siano assegnati i valori della funzione z=f(x,y) in corrispondenza degli MN vertici di una grigliaindichiamo con x1,x2,...,xM e con y1,y2,...,yN le ascisse e le
ordinate dei vertici del reticoloindichiamo con zjk=z(xj,yk) i valori della funzione in
corrispondenza dei vertici del reticolo
y
xx1 x2 x3 xM
y1
y2
yN
Interpolazione bilineareDato il punto (x,y),
indichiamo con j e k gli indici per cui si ha xj<x<xj+1 e yk<y<yk+1
Poniamo:
Interpolazione bilineare:
xj xj+1
yk
yk+1
x
y
kk
k
jj
j
yy
yyu
xx
xxt
11
t(xj+1-xj)
1,1,1,1, )1()1()1)(1(),( kjkjkjkj uzttuzzutzutyxz
Top Related