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