Elementi di Informatica e Applicazioni Numeriche T (2015-2016)

transcript

Elementi di Informatica eElementi di Informatica eApplicazioni Numeriche TApplicazioni Numeriche TInterpolazione ed Estrapolazione

Esempio: Profilo di Velocità di un FluidoEsempio: Profilo di Velocità di un Fluido

Consideriamo un problema realistico

■ Dato un fluido che scorre in una condotta di diametro cm...■ ...Vogliamo determinarne sperimentalmente il profilo di velocità

Metodo:

Utilizziamo un sensore per effettuare una serie di misurazioni

■ Ogni misurazione è nella forma ■ è la posizione del sensore■ è la velocità misurata

i ( , )xi yixiyi

Possiamo disegnare la misurazioni su un piano cartesiano...

Esempio: Velocità di ReazioneEsempio: Velocità di Reazione

Un secondo problema

■ Determinare come la velocità di una reazione chimica...■ ...dipende dalla temperatura a cui essa avviene

Metodo:

Effettuiamo una serie di esperimenti. Per ogni esperimento :i

■ Variamo la temperatura in modo regolare■ Misuriamo la velocità di reazione

Possiamo di nuovo disegnare la misurazioni...

Esempio: Velocità di ReazioneEsempio: Velocità di Reazione

Interpolazione/EstrapolazioneInterpolazione/Estrapolazione

Abbiamo le informazioni che ci interessano per alcuni punti

Ma come facciamo a conoscere la velocità in ogni punto?

Approccio 1: misuriamo la velocità in ogni punto

■ Ovviamente, questo approccio non è pratico

Approccio 1: misuriamo la velocità in ogni punto

■ Ovviamente, questo approccio non è pratico

Approccio 2: inferiamo la velocità

■ Nei punti intermedi■ Prima e dopo le posizioni estreme

I due problemi si chiamano interpolazione ed estrapolazione

Elementi di Informatica eElementi di Informatica eApplicazioni Numeriche TApplicazioni Numeriche T

Polinomio Interpolante

Funzione InterpolanteFunzione Interpolante

Proviamo ad impostare formalmente il problema:

■ Data una serie di di punti ■ HP: i punti sono distinti (i.e. )■ Vogliamo determinare una funzione tale che:

m ( , )xi yi≠ ∀i ≠ jxi xj

f (x)f ( ) = ∀i = 1..mxi yi

Se ci riusciamo, abbiamo risolto tutti i nostri problemi:

■ Possiamo valutare con (interpolazione)■ E possiamo estrapolare i dati nello stesso modo

f (x) < x <xi xj

Grossa difficoltà: come deve essere fatta ?f (x)

Polinomio InterpolantePolinomio Interpolante

Vediamo cosa succede se è un polinomio di grado f (x) n

In questo caso, abbiamo che:

f (x) = ∑j=0

nβj xj

Vediamo cosa succede se è un polinomio di grado f (x) n

In questo caso, abbiamo che:

f (x) = ∑j=0

nβj xj

Per ogni punto , la funzione deve soddisfare:i

f ( ) = =xi ∑j=0

nβj xj

■ Le incognite sono i coefficienti !■ Quindi le equazioni sono lineari

Possiamo quindi usare la forma matriciale:

⎜⎜⎜⎜⎜⎜

m−1xnm

xn−11

xn−12

⋮xn−1

m−1xn−1m

……⋱……

m−1x1m

11⋮11

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

βn−1…β1β0

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

ym−1…y2y1

⎟⎟⎟⎟⎟⎟

La matrice a sx è nota come matrice di VandermondeV

■ Ogni riga si riferisce ad un punto ■ La riga contiene , elevato a tutte la potenze

Proprietà: è non singolare (perché i valori di sono distinti)V xi

⎜⎜⎜⎜⎜⎜

m−1xnm

xn−11

xn−12

⋮xn−1

m−1xn−1m

……⋱……

m−1x1m

11⋮11

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

βn−1…β1β0

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

ym−1…y2y1

⎟⎟⎟⎟⎟⎟

Se abbiamo (più termini nel polinomio che punti):n > m − 1■ Il sistema è sottodeterminato■ Ci sono molte soluzioni possibili■ A noi basterebbe trovarne una, però...

⎜⎜⎜⎜⎜⎜

m−1xnm

xn−11

xn−12

⋮xn−1

m−1xn−1m

……⋱……

m−1x1m

11⋮11

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

βn−1…β1β0

⎟⎟⎟⎟⎟⎟

⎜⎜⎜⎜⎜⎜

ym−1…y2y1

⎟⎟⎟⎟⎟⎟

Se invece (tanti punti quanti termini del polinomio):n = m − 1■ Esiste una sola soluzione, quindi un solo polinomio interpolante■ Per due punti passa una sola retta, per tre punti una sola parabola...

Vediamo un esempio numerico:

Supponiamo di avere i seguenti punti per il nostro fluido in condotta:

(m) (m/s)xi yi

0.010 1.400.054 1.700.111 1.540.138 1.33

Dobbiamo risolvere il sistema:

Vβ = y

Vediamo un esempio numerico:

La matrice di Vandermonde è:

x3i x2

i x1i x0

⎜⎜⎜⎜

0.00100000.15746401.36763102.6280720

0.01000000.29160001.23210001.9044000

0.10000000.54000001.11000001.3800000

1.00000001.00000001.00000001.0000000

⎟⎟⎟⎟

Risolvendo il sistema otteniamo:

f (x) = 282.21 − 144.686 + 15.0746 + 1.26344x3 x2 x1

Vediamo come si comporta nel fare interpolazione (non male!)...

...E poi vediamo come va l'estrapolazione (bene anche qui!)

Polinomio di Lagrange

Metodi Espliciti e ImplicitiMetodi Espliciti e Impliciti

In questo corso facciamo spesso una distinzione tra:

■ Soluzioni in forma chiusa (o metodi espliciti)■ Descrivono la soluzione con una formula

■ Metodi numerici (o metodi impliciti)■ Permettono di calcolare la soluzione con un algoritmo

L'approccio presentato per ottenere il polinomio interpolante

■ Permette di calcolare i coefficienti con un algoritmo...■ ...Ma non di esprimerli in forma chiusa

Si tratta quindi di un metodo numerico

Proviamo ora a vedere una soluzione in forma chiusa

Un Approccio AlternativoUn Approccio Alternativo

Il nostro obiettivo è trovare una funzione tale che:f (x)f ( ) = ∀i = 1..mxi yi

Proviamo ora ad utilizzare una forma diversa, in particolare:

f (x) = (x)∑i=1

myi li

Dove i termini sono dei polinomi(x)li■ Conseguenza: è di fatto ancora un polinomio■ Solo che la descriviamo con una somma di funzioni...■ ...Anziché attraverso i coefficienti

Idea principale: I polinomi in(x)li

f (x) = (x)∑i=1

myi li

...sono definiti in modo che, quando :x = xi

■ Tutti gli con si annullano■ Il termine vale esattamente

( )lj xi j ≠ i( )li xi 1

In questo modo, abbiamo che (come desiderato):

f ( ) = + =xi yi ( )li xi⏟=1∑j≠i

yj ( )lj xi⏟=0

Vediamo i polinomi per il nostro esempio:(x)li

Come sono fatti i polinomi ?(x)li

Ogni deve annullarsi se con (x)li x = xj j ≠ i

■ Possiamo ottenerlo inserendo in l'espressione(x)li(x − )∏

j≠ixj

■ In altre parole il prodotto:

(x − ) … (x − )(x − ) … (x − )x1 xi−1 xi+1 xm

■ Che si annulla per , a meno che x = xj j = i

Come sono fatti i polinomi ?(x)li

Ogni deve valere per (x)li 1 x = xi

■ Possiamo ottenerlo modificando l'espressione precedente:

(x − ) ⟶∏j≠i

xj ∏j≠i

x − xj−xi xj

■ Ogni termine vale se ■ Conseguenza: l'intero prodotto vale

(x − )/( − )xj xi xj 1 x = xi1

Polinomio Interpolante in Forma di LagrangePolinomio Interpolante in Forma di Lagrange

Nel complesso, otteniamo il polinomio di Lagrange:

f (x) = (x)∑i=1

myi li

(x) =li ∏j≠i

x − xj−xi xj

Sarebbe meglio parlare di "polinomio interpolante in forma di Lagrange"

■ È un polinomio di grado (al più) , come nel caso precedente■ Di fatto, è lo stesso polinomio di prima!■ Semplicemente, adesso lo otteniamo con una formula■ In alcuni casi è un grosso vantaggio (esempi nella prossima lezione)!

m − 1

Elementi di Informatica eElementi di Informatica eApplicazioni Numeriche TApplicazioni Numeriche TDue Limitazioni dell'Interpolazione

con Polinomi

L1: Aggiunta di Punti da InterpolareL1: Aggiunta di Punti da Interpolare

Se volessimo migliorare la qualità dell'approssimazione?

Metodo più naturale: incrementare la misurazioni

■ Dividiamo il diametro in segmenti della stessa lunghezza■ Per ogni punto misuriamo la velocità

mxi yi

Idealmente, incrementando l'approssimazione dovrebbe convergere:m

f (x) g(x)− →−−−m→∞

■ Dove è il profilo di velocità realeg(x)In pratica, non è detto che succeda

L1: Funzione di RungeL1: Funzione di Runge

Supponiamo di volere interpolare la funzione di Runge:

g(x) = 11 + 25x2

Supponiamo di non conoscere direttamenteg(x)■ Ne vogliamo però misurare il valore...■ ...Per punti equispaziati nell'intervallo ■ E poi approssimarla con il polinomio interpolante

[−1, 1]

Intuitivamente:

■ Incrementando ...■ ...L'approssimazione dovrebbe migliorare

Questa è la funzione di Runge nell'intervallo [−1, 1]

Questo è il polinomio interpolante con n = 3

L1: Fenomeno di RungeL1: Fenomeno di Runge

■ Le oscillazioni peggiorano con il numero dei campioni■ Il polinomio interpolante diverge vicino agli estremi (i.e. )−1, 1Questo comportamento è noto come fenomeno di Runge

Aumentando il numero di punti, il comportamento ininterpolazione del polinomio può peggiorare

■ È apparentemente paradossale (vedremo poi il motivo)■ È un primo (grosso) limite del polinomio approssimante

L2: Difficoltà di EstrapolazioneL2: Difficoltà di Estrapolazione

Consideriamo ora il secondo dei nostri esempi iniziali

■ Vogliamo stimare come la velocità di una reazione chimica...■ ...Dipenda dalla temperatura a cui essa avviene

Innanzitutto dobbiamo fare degli esperimenti:

■ Facciamo variare la temperature da a ■ Usiamo intervalli di 10 gradi ( punti in tutto)■ Misuriamo la velocità di reazione ad ogni punto

x K300∘ K360∘

Supponiamo di avere ottenuto i punti seguenti:

(°K) (num. puro)xi yi

300 1.01310 1.02320 1.04330 1.08340 1.15350 1.30360 1.44

Questi sono i dati grezzi sul piano cartesiano:

Questo è il polinomio interpolante (in interpolazione)

Questo è il polinomio interpolante (in estrapolazione)

Il comportamento in estrapolazione è pessimo!

■ Estrapolare è più difficile che interpolare in generale■ C'è un margine di incertezza maggiore

In questo caso però la situazione è particolarmente brutta, perché:

Aumentando il numero di punti, il comportamento inestrapolazione del polinomio interpolante può peggiorare

■ Questo è un secondo limite del polinomio approssimante

Metodo Lineare dei Minimi Quadrati

Alcune ConsiderazioniAlcune Considerazioni

Abbiamo osservato che, aumentando il numero di punti :( , )xi yi

■ Il comportamento in interpolazione può peggiorare■ Il comportamento in estrapolazione può peggiorare

La radice dei due problemi è comune.

Alcune ConsiderazioniAlcune Considerazioni

Abbiamo osservato che, aumentando il numero di punti :( , )xi yi

■ Il comportamento in interpolazione può peggiorare■ Il comportamento in estrapolazione può peggiorare

La radice dei due problemi è comune. In particolare:

I due problemi nascono dal fatto che il gradodel polinomio interpolante diventa troppo alto

■ I polinomi di grado elevato possono variare in modo brusco...■ ...E si prestano a causare errori dovuti alla precisione finita

Utilizzo di Polinomi di Grado InferioreUtilizzo di Polinomi di Grado Inferiore

Cosa succede se utilizziamo un polinomio di grado ?n < m − 1Il nostro sistema:

Vβ = y

■ è la matrice di Vandermonde corrispondente ad ■ è il vettore dei coefficienti

Diventa sovradeterminato (e probabilmente impossibile)

■ Se vogliamo utilizzare un polinomio di grado inferiore...■ ...Dobbiamo cambiare la formulazione del problema

Utilizzo di Polinomi di Grado InferioreUtilizzo di Polinomi di Grado Inferiore

Vediamo le idee principali del nuovo approccio:

■ Stabiliamo il grado del polinomio a priori (HP: )■ Non richiediamo che il polinomio passi per tutti i punti■ Cerchiamo invece di minimizzare una stima di errore

n n ≤ m − 1

In particolare, introduciamo il concetto di residuo

= − f ( )ri yi xi

■ Dove, ancora una volta:

f (x) = ∑j=0

nβj xj

Residui e Stima dell'ErroreResidui e Stima dell'Errore

In sostanza il residuo è il nostro errore di stima per ri ( , )xi yi

Come ottenere una stima dell'errore globale?

Idea 1: utilizziamo una somma

e = ∑i=1

■ Una pessima idea in pratica■ Residui con segno diverso si possono cancellare

Idea 2: somma dei valori assoluti

e = ∑i=1

m∣∣ri∣∣

■ Già meglio!■ Però è difficile da trattare (derivata discontinua)

Idea 3: somma dei residui al quadrato

e = ∑i=1

■ La derivata è continua■ Vantaggio addizionale: gli errori più grandi contano di più!

Minimizzazione ai Minimi QuadratiMinimizzazione ai Minimi Quadrati

Vogliamo che l'errore sia più piccolo possibile

Quindi si tratta di risolvere il problema:

min e(β) = (β∑i=1

mri )2

Noto come problema di minimizzazione ai minimi quadrati

■ Nel problema abbiamo evidenziato la dipendenza di da ■ D'ora in poi, lo faremo sempre■ La soluzione del problema definisce la funzione approssimante■ Il corrispondente errore indica la qualità dell'approssimazione

Metodo (Lineare) dei Minimi QuadratiMetodo (Lineare) dei Minimi Quadrati

Nel nostro caso, per ogni punto , abbiamo:( , )xi yi

(β) = −ri yi ∑j=0

βj xji

■ Conseguenza: i residui dipendono linearmente da ■ Per questa ragione, il metodo di soluzione che vedremo...■ ...Si chiama metodo lineare dei minimi quadrati

Per via della linearità, possiamo definire i residui in forma matriciale:

r = y − βV⏟m×(n+1)

■ Dove è la matrice di VandermondeV

Pertanto il nostro problema diventa:

min e(β) = r(β r(β))T

■ è un vettore riga e un vettore colonna■ Il prodotto scalare è la somma dei residui al quadrato

rT rrrT

Di qui otteniamo:

e(β) = (y − Vβ (y − Vβ))T

E quindi:

e(β) = y − Vβ − (Vβ y − (Vβ VβyT yT )T )T

A questo punto osserviamo che, nell'equazione:

e(β) = y − Vβ − (Vβ y − (Vβ VβyT yT )T )T

■ Tutti i termini sono scalari e possono essere trasposti a piacimento■ Per il prodotto matriciale, vale che (AB =)T BT AT

Possiamo quindi ottenere:

e(β) = y − ( Vβ − (Vβ y + (Vβ VβyT yT )T )T )T

e(β) = y − (Vβ y − (Vβ y + (Vβ VβyT )T )T )T

e(β) = y − 2 y + VβyT βT V T βT V T

A questo punto per trovare un punto estremo possiamo risolvere

∇e(β) = 0■ Poi dovremmo valutare il numero di soluzioni possibili...■ ...E verificare se corrispondano a massimi o minimi

A questo punto per trovare un punto estremo possiamo risolvere

∇e(β) = 0■ Poi dovremmo valutare il numero di soluzioni possibili...■ ...E verificare se corrispondano a massimi o minimi

Per poter calcolare in forma matriciale, è utile:∇e(β)■ Ottenere una formula per la "derivata" di un prodotto scalare ■ Ottenere una formula per la "derivata" di una forma quadratica

Nelle prossime due slides vedremo come farlo...

Digressione: Derivata di un Prodotto ScalareDigressione: Derivata di un Prodotto Scalare

Consideriamo un prodotto scalare xvT

■ Tecnicamente, si tratta di una funzione :f : → ℝℝn

x =vT ∑i=1

■ Il gradiente della funzione è dato da:

∇( x) = = = vvT

⎜⎜⎜⎜

∂( x)vT

⋮∂( x)vT

⎟⎟⎟⎟

⎝⎜⎜⎜

⎠⎟⎟⎟

Digressione: Derivata di una Forma QuadraticaDigressione: Derivata di una Forma Quadratica

Consideriamo una forma quadrativca AxxT

■ Corrisponde ancora ad una funzione :f : → ℝℝn

Ax =xT ∑j=1

∑i=1

naijxixj

■ Derivando rispette ad otteniamo:xk

= +∂( Ax)xT

xk ∑j=1

nakjxj ∑

naikxi

Digressione: Derivata di una Forma QuadraticaDigressione: Derivata di una Forma Quadratica

Nell'espressione:

= +∂( Ax)xT

xk ∑j=1

nakjxj

∑i=1

naikxi

I due termini (A) e (B) corrispondono a:

(A) Il prodotto scalare di una colonna di per (B) Il prodotto scalare di una riga di per

A xA x

Quindi, in forma matriciale abbiamo:

∇( Ax) = ( + A)xxT AT

A questo punto possiamo riprendere la discussione

Per trovare un polinomio approssimante, dobbiamo risolvere:

∇e(β) = 0Dove:

e(β) = y − 2 y + VβyT βT V T βT V T

e(β) = y − 2 β + βyT ( yV T )T

⏟vettore

βT VV T⏟matrice

Quindi, applicando le due regole di differenziazione matriciale:

−2( y) + ( V + V )β = 0V T V T V T

La matrice è simmetrica per costruzione, infatti:VV T

( V = VV T )T V T

Quindi a partire da:

−2( y) − ( V + V )β = 0V T V T V T

Possiamo ottenere:

−2( y) − 2( V)β = 0V T V T

E finalmente otteniamo un sistema lineare!

Vβ = yV T V T

■ Poiché è non-singolare...■ ...Anche e quindi sono non singolari

VV T VV T

Quindi il sistema lineare:

β =VV T⏟(n+1)×nV T

⏟(n+1)×my⏞

■ Ha una sola soluzione...■ ...Che corrisponde ad un estremo per e(β)Ci rimane da capire se questa corrisponda a un minimo...

Per capire se la soluzione del sistema corrisponda al minimo...

...Usiamo le normali regole per le funzioni multivariate:

■ Calcoliamo la matrice Hessiana di ...■ ...E ci assicuriamo che sia definita positiva

Per capire se la soluzione del sistema corrisponda al minimo...

...Usiamo le normali regole per le funzioni multivariate:

■ Calcoliamo la matrice Hessiana di ...■ ...E ci assicuriamo che sia definita positiva

Per ottenere possiamo derivare nuovamente:H(e(β))∇e(β) = Vβ − yV T V T

■ Applicando la regola per alle righe di otteniamo∇ xvT VV T

H(e(β)) = VV T

Si può dimostrare (abbastanza facilmente) che:

■ Se è non singolare (come nel nostro caso)...■ ...Allora è definita positiva

Finalmente, concludiamo che risolvendo il sistema lineare:

Vβ = yV T V T

Otteniamo un vettore di terminiβ n + 1■ Questi sono i coefficienti di un polinomio di grado ...■ ...Che minimizza la somma dei residui al quadrato

Questo è il metodo lineare dei minimi quadrati

Matrice Pseudo-InversaMatrice Pseudo-Inversa

A volte il sistema viene scritto nella forma:

β = ( V yV T )−1V T

■ Dove è una matrice ■ È nota come matrice pseudo-inversa di

( VV T )−1V T (n + 1) × mV

Perché parlarne? Perché se in Octave eseguite:

b = V \ y % divisione sinistra

■ Se è quadrata, viene risolto il sistema lineare■ Se è quadrata, viene risolto problema ai minimi quadrati

Di fatto eseguire questo metodo in Octave è molto facile!

Metodo dei Minimi Quadrati per l'Esempio 1Metodo dei Minimi Quadrati per l'Esempio 1

Approssimazione di grado del profilo di velocità:2

Metodo dei Minimi Quadrati per l'Esempio 2Metodo dei Minimi Quadrati per l'Esempio 2

Approssimazione di grado della velocità di reazione:3

Minimi Quadrati: Funzione di RungeMinimi Quadrati: Funzione di Runge

Approssimazione di grado della funzione di Runge:2

Minimi Quadrati: Funzione di RungeMinimi Quadrati: Funzione di Runge

In questo caso non funziona molto bene...

Oltre i Polinomi

Approssimazione con Funzioni Non LineariApprossimazione con Funzioni Non Lineari

Nel metodo lineare dei minimi quadrati:

■ Approssimiamo una funzione ...■ ...Con una funzione approssimante

g(x)f (x, β)

Nel nostro caso era un polinomio e i coefficientif (x) β

Ma potremmo usare anche qualcosa diverso dai polinomi!

Possiamo usare qualunque funzione

■ Purché dipenda linearmente da ...■ ...dove i termini di si chiamano genericamente parametri

Per la precisione, possiamo usare qualunque nella forma:f (x)

f (x, β) = (x)∑j=1

nβjϕj

■ Dove è il numero dei parametri■ E ogni è una funzione arbitraria di

n(x)ϕj x

In questo caso, i residui prendono la forma:

(β) = − ( )ri yi ∑j=1

nβjϕj xi

■ Dove può essere pre-calcolato( )ϕj xi

In forma matriciale abbiamo:

r(β) = y − Φβ

Dove ogni termine della matrice corrisponde a Φi,j ( )ϕj xi

■ La formulazione è analoga a quella dei polinomi■ Semplicemente abbiamo invece di Φ V

Possiamo ottenere risolvendo:β

Φβ = yΦT ΦT

Esattamente come nel caso dei polinomi!

■ Condizione: deve essere non singolareΦ

Esempio: Funzione di RungeEsempio: Funzione di Runge

Abbiamo avuto difficoltà ad approssimare la funzione di Runge

Proviamo ad utilizzare una funzione approssimante diversa

Per esempio:

f (x, β) = + + +β1x1 β2 ( )x0⏟=1β3

11 + |x| β4

11 + x2

■ La funzione non è un polinomio■ Però dipende linearmente dai parametri β

Quindi possiamo applicare il metodo lineare dei minimi quadrati

Assumendo di avere x = [−1, −0.33, 0.33, 1]■ Calcoliamo la matrice Φ

x1i x0

1+∣∣xi ∣∣1

⎜⎜⎜⎜

−1.00000−0.333330.333331.00000

1.000001.000001.000001.00000

0.500000.750000.750000.50000

0.500000.900000.900000.50000

⎟⎟⎟⎟

■ Poi calcoliamo:

β = ( Φ) yΦT ΦT

Otteniamo questa funzione approssimante (non un gran che):

Aumentando il numero di punti, però...

...L'approssimazione migliora considerevolmente!

Elementi di Informatica e Applicazioni Numeriche T (2015-2016)

Documents