Approssimazione di dati e funzioni -...

GeneralitaInterpolazione polinomiale

ConvergenzaInterpolazione polinomiale a tratti

Metodo dei minimi quadrati

Approssimazione di dati e funzioni

Stefano Berrone

Dipartimento di Matematicatel. 011 0907503

[email protected]

http://calvino.polito.it/~sberrone

Laboratorio di modellazione e progettazione materiali

Stefano Berrone Approssimazione di dati e funzioni




Generalita

Problema 1 Dati (xi , yi ) i = 0, ..., n (es. misurazioni) voglioapprossimarli con una funzione g in modo da poter stimarel’andamento dei dati anche in punti x 6= xi .Problema 2 Data f voglio approssimarla con una g piu facile dausare (ad esempio un polinomio). Es:∫ b

aexp(−x2)dx =?

Se exp(−x2) ' g e g so integrarla facilmente,∫ b

aexp(−x2)dx '

∫ b

ag(x)dx

NB Problema 2 si tratta come il Problema 1 con dati (xi , f (xi )),i = 0, ..., n





Step: data f ∈ F = C k [a, b]1 Individuare un sottospazio Fm di funzioni di dimensione finita

in cui scegliere la funzione fm che approssima f1 Pm, polinomi algebrici di grado m:

fm(x) = pm(x) =∑m

k=0 akxk

2 PTm,funzioni polinomiali a tratti di grado r (tipicamente basso):

fm|[xi ,xi+1] = pr (x)(i)

3 Sr , funzioni spline di ordine r (vedi dopo, caso particolare di 2)4 Πm, polinomi trigonometrici:πm(x) = a0 +

∑mk=1 (akcos(kωx) + bksin(kωx))

5 . . .

2 Individuare un criterio per determinare fm ∈ Fm

1 interpolazione: fm(xi ) = yi , i = 0, ..., n2 minimi quadrati: min

∑ni=0(fm(xi )− yi )

2 (ideale perapprossimazione di dati sperimentali)





Dimensioni...

1 dim(Pm) = m + 1

2 dim(PTm) = (r + 1)×(numero di intervalli)

3 dim(Sm) = ...

4 dim(Πm) = 2m + 1

5 . . .

NB Remind: Dimensione = # elementi della base di Fm = #parametri che individua fm ∈ Fm = # di condizioni da imporre(interpolazione: #nodi = n + 1 = m + 1 = grado polinomio + 1)





Misura dell’errore

Norma di funzione:

‖f ‖∞ = supx∈[a,b]

|f (x)| = maxx∈[a,b]

|f (x)|

Facendo crescere m mi aspetto che l’approssimazione migliori.

Definizione (convergenza)

Dati {Fm}m≥0, e fm ∈ Fm, se

limm→∞

‖f − fm‖∞ = 0

si dice che si ha convergenza (uniforme) di fm a f (i.e.convergenza dell’approssimazione)





Base di LagrangeBase di Newton

Dati: (xi , yi ), i = 0, . . . , n (xi : nodi)Target: pn(x) =

∑nk=0 akφk(x)

Pn = span{φ0(x), ...φn(x)}Esempio: Pn = span{1, x , x2, ..., xn}Criterio di interpolazione:pn(xi ) = yi i = 0, ..., n ⇒

∑nk=0 akφk(xi ) = yi

Posto:

G =

φ0(x0) φ1(x0) . . . φn(x0)φ0(x1) φ1(x1) . . . φn(x1)

...φ0(xn) φ1(xn) . . . φn(xn)

∈ R(n+1)×(n+1)

a = (a0, ..., an)T ∈ Rn+1, y = (y0, ..., yn)T ∈ Rn+1

=⇒ Ga = y G matrice di Grahm






Esistenza e unicita polinomio interpolante ⇔ invertibilita di G

Proposizione

Se φk(x), k = 0, ..., n linearmente indipendenti e xi 6= xj per i 6= j(nodi distinti) una matrice della forma di G e invertibile.

Poiche φk(x) formano una base, nodi distinti garantisconoesistenza e unicita di pn(x).






Base di Lagrange

Come trovo ak? NON risolvo il sistema lineare:1 non conviene computazionalmente2 puo essere molto mal condizionato (es. G = matrice di

Vandermonde con base monomiale)

Basi diverse =⇒ metodi diversi per determinare coefficienti dellacombinazione lineare Base di Lagrange: φk(x) = `k(x) t.c.:

`k(xi ) = δik =

{1 i = k

0 i 6= k

yi = pn(xi ) =n∑

k=0

ak`k(xi ) = ai`i (xi ) = ai

⇓ai = yi !!!






Come e fatto `k(x)?

`k(x) =

∏i 6=k(x − xi )∏i 6=k(xk − xi )

Esempio

(x0, y0) = (−2, 3), (x1, y1) = (1,−7), (x2, y2) = (3,−5):

`0(x) =(x − 1)(x − 3)

(−2− 1)(−2− 3), `1(x) =

(x + 2)(x − 3)

(1 + 2)(1− 3)

`2(x) =(x + 2)(x − 1)

(3 + 2)(3− 1)

p2(x) = 3`0(x)− 7`1(x)− 5`2(x)






Sia pn(x) interpolante (xi , yi ), i = 0, ..., nAggiungo (xn+1, yn+1) e cerco di scrivere pn+1(x) sfruttando illavoro gia fatto

pn+1(x) = pn(x) + qn(x), q(x) =???

Che sappiamo di qn(x)?

1 e un polinomio di grado n + 1

2 pn+1(xi ) = pn(xi ) + qn+1(xi ), i = 0, ..., n ⇓yi = yi + qn+1(xi ) =⇒ qn+1(xi ) = 0 ⇓qn+1(x) ha come radici xi , i = 0, ..., n

qn+1(x) = an+1(x − x0)(x − x1) . . . (x − xn) = an+1ωn+1(x)

an+1 =yi+1 − pn(xi+1)

ωn+1(xi+1)






pn+1(x) = pn(x) + an+1ωn+1(x)

pn+1(x) = pn−1(x)+anωn(x) + an+1ωn+1(x)

pn+1(x) =n+1∑k=0

akωk(x)

Base di Newton:

ω0(x) ≡ 1

ω1(x) = x − x0

ω2(x) = (x − x0)(x − x1)...

ωn+1(x) = Πnk=0(x − xk)

Calcolo dei coefficienti ak in modo semplificato?Stefano Berrone Approssimazione di dati e funzioni





Differenze divise

Siano xi , i = 0, ..., n nodi distinti

Definizione

Differenza divisa di ordine 0 relativa al nodo xk

f [xk ] := f (xk)

Definizione

Differenza divisa di ordine k relativa ai nodi xi0 , xi1 , ..., xik

f [xi0 , xi1 , ..., xik ] :=f [xi1 , ..., xik ]− f [xi0 , xi1 , ..., xik−1

]

xik − xi0

Proprieta

ak = f [x0, x1, ..., xk ]






Tavola delle differenze divise

Esempio

Siano (x0, y0) = (−2, 2), (x1, y1) = (1,−7), (x2, y2) = (3,−5),(x3, y3) = (4,−7)

xk f [xk ] f [xk , xk+1] f [xk , xk+1, xk+2] f [xk , xk+1, xk+2, xk+3]−2 21 −7 −7−2

1+2 = −3

3 −5 −5+73−1 = 1 1+3

3+2 = 45

4 −7 −7+54−3 = −2 −2−1

4−1 = −1−1− 4

54+2 = − 3

10

p3(x) = 2− 3(x + 2) +4

5(x + 2)(x − 1)− 3

10(x + 2)(x − 1)(x − 3)






Proprieta delle differenze divise

Le differenze divise non dipendono dall’ordine dei nodi

f [xi0 , ..., xik ] = f [P(xi0 , ..., xik )]

dove P(xi0 , ..., xik ) e una qualunque permutazione dei nodi.

Osservazione

Se i nodi sono distinti, i denominatori delle differenze divise sonosempre diversi da zero tavola ben definita






Interpolazione di Hermite

Dati: (xi , f (xi )), i = 0, ..., n e (xi , f′(xi )), i = 0, ..., n La forma di

Newton e particolarmente efficiente

xk f [xk ] f [xk , xk+1] f [xk , xk+1, xk+2]x0 f (x0)x0 f (x0) f ′(x0)

x1 f (x1) f [x0, x1] f [x0,x1]−f ′(x0)x1−x0

x1 f (x1) f ′(x1)...

...

Generalizzabile a piu derivate e non su tutti i nodi le stesse derivate

Esercizio proposto

Esempi 5.1 e 5.2 libro di testo (Monegato)





Convergenza dell’interpolazione polinomiale

Errore di interpolazione:

En(x) = f (x)− pn(x)

Sicuramente:1 En(xi ) = 02 En(x) ≡ 0 se f e un polinomio di grado ≤ n

Teorema

Sia f ∈ Cn+1([a, b]), En(x) = f n+1(ξ)(n+1)! ωn+1(x), con ξ ∈ [a, b].

Ma in generale che accade a

limn→∞

‖En(x)‖∞?

Dipende dalla scelta nodi e dalle caratteristiche di f !Stefano Berrone Approssimazione di dati e funzioni




Creiamo una successione di nodi:

x(1)0

x(2)0 x

(2)1

...

x(n+1)0 x

(n+1)1 . . . x

(n+1)n

...Teorema

Sia f ∈ C 0([a, b]),||En||∞ ≤ (1 + Λn) min

q∈Pn

||f − q||∞, Λn costante di Lebesgue.

1 nodi equidistanti: Λn ≥ en2

2 nodi di Chebichev xi = cos( 2i+12(n+1)π) ∈ (−1, 1), i = 0, ..., n

Λn ∼ π2 log n





Esempio (Fenomeno di Runge)

f (x) = 11+x2 , x ∈ [a, b]

x ∈ [1, 2], nodi equidistantiStefano Berrone Approssimazione di dati e funzioni




x ∈ [−5, 5], nodi equidistanti





x ∈ [−5, 5], nodi Chebichev





Teorema ( ¨ )

Se f ∈ C 0[a, b] esiste una successione di nodi di interpolazione in[a, b] tale che la successione di polinomi interpolanti costruita suessi converge uniformemente a f

Teorema (Faber, _)

Data una qualunque successione di nodi di interpolazione distinti in[a, b] esiste sempre f ∈ C 0[a, b] per cui la successione di polinomiinterpolanti costruita su essi non converge uniformemente a f

Teorema (Bernstein, ¨ )

Se f ∈ C 1[a, b] la successione di polinomi interpolanti costruiti sunodi di Chebycev converge uniformemente a f in [a, b].Se f ∈ C 2[a, b], si ha inoltre ‖En‖∞ = O( 1√

n)





Teorema ( ¨ )

Se f (z) e analitica in un dominio del piano complesso che contiene[a, b] e le singolarita di f (z) distano da [a, b] piu di b − a, alloralimn→∞ ‖En‖∞ = 0.





Interpolazione polinomiale a tratti (non raccordata)

Interpolazione polinomiale: infittire nodi non garantisce migliorcomportamento di pn(x)Interpolazione polinomiale a tratti:

1 fissiamo a priori il grado r (piccolo) del polinomio2 partizioniamo [a, b] in tanti intervallini tramite n + 1 nodi xi :

a = x0 ≤ x1 ≤ . . . ≤ xn = b

3 ogni r + 1 nodi, costruiamo unpolinomio diverso

4 infittire i nodi non significaaumentare r ma il numero di trattidi polinomio

5 se f e regolare convergeuniformemente





Interpolazione tramite spline

Definizione (spline)

Siano xi nodi che formano una partizione di [a, b]. Fissato d ≥ 1,sd(x) e una funzione spline relativa alla partizione di [a, b] se:

1 sd(x) e un polinomio di grado (al piu) d su ogni [xi , xi+1]

2 s(k)d (x) e una funzione continua in [a, b] per ogni

k = 0, ..., d − 1

Definizione (spline interpolante)

Dati (xi , yi ), i = 0, ..., n, sd(x) e una spline interpolante i dati se:

1 sd(x) e una funzione spline relativa alla partizione {xi}2 sd(xi ) = yi , i = 0, .., n





Gradi di liberta

# parametri: (d + 1)× n = dn + n# condizioni: n + 1 + d × (n − 1) = dn + d + n − 1n + 1:interpolazione + (n − 1) ∗ d :raccordo nei nodi interniDifferenza: d − 1 Nel caso piu utilizzato (d = 3, spline cubiche)mancano 2 condizioni Le spline cubiche vengono dette...

1 naturali: s ′′3 (x0) = 0, s ′′3 (xn) = 0

2 periodiche: s ′3(x0) = s ′3(xn), s ′′3 (x0) = s ′′3 (xn)

3 vincolate: s ′3(x0) = y ′0, s ′3(xn) = y ′n, con y ′0, y ′n dati

4 not–a–knot: s ′′′3 (x) continua anche in x1 e xn−1





Convergenza delle spline cubiche

Teorema

Sia s3(x) la spline cubica interpolante i dati (xi , yi ) per i = 0, ..., ncon condizioni di tipo 1,2,3. Sia h = maxi hi , hi = xi − xi−1. Sef ∈ C 2[a, b] per h→ 0 si ha

‖s(p)3 − f (p)‖∞ = o(h2−p), p = 0, 1, 2

Se f ∈ C k [a, b], k = 3, 4 e hhi≤ cost. per h→ 0 si ha

‖s(p)3 − f (p)‖∞ =

{o(h3−p) k = 3

O(h4−p) k = 4p = 0, 1, 2, 3

Per regolarita di f maggiori la velocita di convergenza non migliora.





Proprieta delle spline cubiche

Teorema

Tra tutte le funzioni f ∈ C 2[a, b] che assumono valori yi nei nodixi e soddisfacenti condizioni di tipo 1,2,3, le spline cubiche sono lesole funzioni che minimizzano l’integrale

E (f ) =

∫ xn

x0

[f ′′(x)]2dx

Le spline naturali godono di una proprieta di minimo assoluto.

f ′′(x)

(1 + f ′(x)2)3/2= curvatura di f nel punto x

⇒ se |f ′(x)| e sufficientemente piccolo E (f ) e unaapprossimazione della curvatura globale di f in [a, b].





Comandi MATLAB

1 a = polyfit(x,y,n)

2 yy = polyval(a,xx)

3 yy = spline(x,y,yy)





Limiti dell’interpolazione:

interpolazione polinomiale globale non adatta per n grandi

interpolazione in generale non adatta per dati sperimentali(molti dati + affetti da errore)

Fissiamo a priori un modello per i nostri dati in uno spazio didimensione m + 1 (tipicamente m� n): se φi (x), i = 0, ...,msono m funzioni linearmente indipendenti cerchiamo

fm(x) = c0φ0(x) + c1φ1(x) + . . .+ cmφm(x)

Esempi:

f1(x) = c0 + c1x

f2(x) = c0 + c1x + c2x2

f3(x) = c0 + c1x + c2 cos(x) + c3 sin(x)

f1(x) = c0 + c1exp(x)

. . .Stefano Berrone Approssimazione di dati e funzioni




Come si determinano i coefficienti ck?

fm(xi ) = yi NO

minc∈Rm+1

n∑i=0

(fm(xi )− yi )2

minc∈Rm+1

n∑i=0

(m∑

k=0

ckφk(xi )− yi

)2

Definizione

fm(x) e la migliore approssimazione nel senso dei minimiquadrati





Notazione matriciale: posto

A =

φ0(x0) φ1(x0) . . . φm(x0)φ0(x1) φ1(x1) . . . φm(x1)

...φ0(xn) φ1(xn) . . . φm(xn)

∈ Rn+1×m+1

y = (y0, . . . , yn)T ∈ Rn+1, c = (c1, . . . , cm)T ∈ Rm+1

minc∈Rm+1

n∑i=0

m∑

k=0

ckφk(xi )︸︷︷︸Ac

−yi

2

= minc∈Rm+1

‖Ac − y‖22





‖Ac − y‖22 funzionale quadratico convesso: minimo assunto in

corrispondenza degli zeri del gradiente, quindi per c che soddisa

∇‖Ac − y‖22 = AT Ac − AT y = 0

AT Ac = AT y e detto sistema delle equazioni normali

Proprieta

AT A simmetrica semidefinita positiva. Se colonne di A linearmenteindipendenti, AT A simmetrica definita positiva. Se φk(x)linearmente indipendenti e nodi distinti, sicuramente colonne di Alinearmente indipendenti.





Esempio

(−1, 2), (1, 3), (3,−5), (2, 7), f2(x) = c0 + c1x + c2x2

A =

1 −1 11 1 11 3 91 2 4

, y =

23−57

AT A =

4 5 155 15 35

15 35 99

, AT y =

70−12

Soluzione sistema lineare AT Ac = AT y: c = ( 11

2 ,94 ,−

74 )T .





Approssimazione polinomiale:

AT A =

∑

i 1∑

i xi . . .∑

i xmi∑

i xi∑

i x2i . . .

∑i xm+1

i...∑

i xmi

∑i xm+1

i . . .∑

i x2mi

∈ Rm+1×m+1

AT y =

∑

i yi∑i yixi...∑

i yixmi

∈ Rm+1





Svantaggi dell’approccio precedente:

costoso costruire AT A

tipicamente mal condizionata

Alternativa computazionale: uso della fattorizzazione QR





Proprieta

Matrici ortogonali (QT Q = QQT = I ) non cambiano normaeuclidea di vettore:

‖Qx‖22 = (Qx)T Qx = xT QT Qx = xT x = ‖x‖2

2

Sia A = QR ∈ Rn+1×m+ R = QT A ∈ Rn+1×m+1,Q ∈ Rn+1×n+1

R =

(R1

0

)Rm+1×m+1

R(n−m)×m+1 QT y =

(z1

z2

)Rm+1

Rn−m

‖Ac − y‖22 = ‖QT (Ac − y)‖2

2 = ‖QT Ac − QT y‖22 =

= ‖Rc − QT y‖22 = ‖Rc − z1‖2

2 + ‖z2‖22

min ‖Ac−y‖22 ⇔ min ‖Rc−z1‖2

2+‖z2‖22 ⇔ min ‖Rc−z1‖2

2 ⇔ Rc = z1





Comandi MATLAB

polyfit se φk(x) polinomiali





Un esempio nonlineare

Modello: a1ea2x

NON e nella forma∑

ckφk(x)! Ma con una trasformazione divariabili...

y = a1ea2x

log(y)︸︷︷︸z

= log(a1ea2x) = log(a1)︸︷︷︸c1

+ a2︸︷︷︸c2

x

z = c1 + c2x

Quindi con i dati (xi , log(yi )) posso ricostruire c1 e c2 e quindirisalire a a1 = ec1 e a2 = c2


Approssimazione di dati e funzioni -...

Documents

Transcript of Approssimazione di dati e funzioni -...