Interpolazione Polinomiale Esatta

Per il composto acido formico:

• Vogliamo determinare come varia la pressione di saturazione \(P^*\)…
• …In funzione della temperatura \(T\), sopra il punto di fusione

Abbiamo a disposizione i seguenti dati sperimentali:

• Vogliamo determinare la pressione per ogni temperatura

Un possibile approccio:

• Possiamo provare a tracciare una curva…
• …Che soddisfi tutte le osservazioni sperimentali

Interpolazione Polinomiale Esatta

Innanzitutto preparo i dati:

% Temperatura (°C) e Pressione (Torr)
T = [-20,-5,2.1,10.3,24,32.4,43.8,61.4,80.3,100.6];
P = [1, 5, 10, 20, 40, 60, 100, 200, 400, 760];
Tf = 8.2; % Temperatura di fusione (°C)
% Trasformo le temperature in °K
T = T + 273.15;
Tf = Tf + 273.15;
% Considero i punti sopra la temperatura di fusione
P = P(T > Tf) % Uso T per selezione una parte di P
T = T(T > Tf) % Quindi rimpiazzo T

• Nota: per restringere il campo ai valori sopra la temperature di fusione…
• …Abbiamo usato una indicizzazione con vettore di valori logici

Interpolazione Polinomiale Esatta

Che tipo di curva andremo a costruire?

• Potremmo usare un polinomio…
• …Con un parametro per ogni osservazione sperimentale

Ragionando così, dovremmo usare un polinomio di 7° grado

\[\alpha_7 t_0^7 + \alpha_6 t_0^6 + \alpha_5 t_0^5 + \alpha_4 t_0^4 + \alpha_3 t_0^3 + \alpha_2 t_0^2 + \alpha_1 t_0 + \alpha_0 = p_0 \\ \alpha_7 t_1^7 + \alpha_6 t_1^6 + \alpha_5 t_1^5 + \alpha_4 t_1^4 + \alpha_3 t_1^3 + \alpha_2 t_1^2 + \alpha_1 t_1 + \alpha_0 = p_1 \\ \alpha_7 t_2^7 + \alpha_6 t_2^6 + \alpha_5 t_2^5 + \alpha_4 t_2^4 + \alpha_3 t_2^3 + \alpha_2 t_2^2 + \alpha_1 t_2 + \alpha_0 = p_2 \\ \ldots\]

• Un parametro \(\alpha_j\) per ogni monomio
• Una equazione per ogni coppia \(t_i, p_i\) di misurazioni

Interpolazione Polinomiale Esatta

A questo punto possiamo impostare un sistema lineare:

• Per costruire la matrice dei coefficienti…
• …Sfruttiamo il fatto che tutte le equazioni hanno la stessa struttura

Possiamo costruirla per concatenazione di colonne!

T = T' % T diventa un colonna
A = [T.^7, T.^6, T.^5, T.^4, T.^3, T.^2, T, T.^0]

A questo punto possiamo risolvere il sistema:

P = P' % il vettore delle pressioni diventa una colonna
x = A \ P

Interpolazione Polinomiale Esatta

Questo tipo di problema è molto comune

Si parla di interpolazione polinomiale esatta

• Trovare una curva polinomiale…
• …Che intercetti una serie di punti noti…
• …E non debba soddisfare alcuna condizione addizionale

Matlab mette a disposizione una funzione per risolverli:

P = polyfit(X, Y)

• X è il vettore con le coordinate \(x\)
• Y è il vettore con le coordinate \(y\)
• P contiene i coefficienti del polinomio interpolante

Interpolazione Polinomiale Esatta

Proviamo ad usarla:

p = polyfit(T, P, length(T)-1);

figure();
t = linspace(260, 400);
plot(t, polyval(p, t), 'linewidth', 2);
hold on
scatter(T, P, 'linewidth', 2);
hold off
grid()
xlabel('T')
ylabel('P')

Interpolazione Polinomiale Esatta

Purtroppo, il risultato lascia desiderare…

• All’“esterno” dei punti noti, l’andamento non ha senso fisico
• Il polinomio è di grado alto e può variare in modo troppo brusco

Interpolazione Lineare a Tratti

Come possiamo ovviare il problema?

Un primo approccio semplicissimo:

• Tra ogni coppia di misurazioni…
• …Assumiamo che la funzione sia lineare

Così otteniamo una approssimazione lineare a tratti

In Matlab la possiamo calcolare con:

y = interp1(X, Y, x)

• X, Y sono le coordinate \(x\) ed \(y\) dei punti noti
• x contiene i valori di \(x\) per cui si desidera la valutazione
• y è il valore dell’approssimazione nei punti desiderati

Interpolazione Lineare a Tratti

Nel nostro caso, potremmo avere:

% Funzione che calcola l'approssimazione
f = @(t) interp1(T, P, t)

figure();
t = linspace(260, 400);
plot(t, f(t), 'linewidth', 2);
hold on
scatter(T, P, 'linewidth', 2);
hold off
grid()
xlabel('T')
ylabel('P')
title('Interpolazione lineare a tratti');

Interpolazione Lineare a Tratti

Va molto meglio!

• Ma cosa succede all’esterno del range?
• E se abbiamo poche misurazioni? O se sono inaffidabili?

Metodo Lineare dei Minimi Quadrati

Idealmente, dovremmo limitare il grado dell’approssimazione

Per esempio, potremmo impostare:

\[\alpha_2 t_0^2 + \alpha_1 t_0 + \alpha_0 = p_0 \\ \alpha_2 t_1^2 + \alpha_1 t_1 + \alpha_0 = p_0 \\ \alpha_2 t_2^2 + \alpha_1 t_2 + \alpha_0 = p_0 \\ \alpha_2 t_3^2 + \alpha_1 t_3 + \alpha_0 = p_0 \\ \ldots\]

Alle equazioni corrisponde un sistema:

\[A x = b\]

• Purtroppo il sistema che si ottiene è tipicamente sovradeterminato
• Troppe equazioni, pochi gradi di libertà \(\Rightarrow\) nessuna soluzione

Metodo Lineare dei Minimi Quadrati

Metodo dei minimi quadrati

Potremmo però risolvere il sistema in senso approssimato:

• Misuriamo l’errore in base ai residui del sistema lineare

\[E = Ax - b\]

• Minimizziamo il quadrato della norma dell’errore:

\[\min_x \|E\|^2 = \min_x \|Ax-b\|^2\]

Perché il quadrato dell’errore?

• Considera come errore sia i residui positivi che quelli negativi
• Penalizza i grandi errori in modo maggiore
• Ma soprattutto, si deriva facilmente!

Metodo Lineare dei Minimi Quadrati

Il problema di minimizzazione:

\[\min_x \|Ax-b\|^2\]

Si riduce ad un secondo sistema lineare:

\[A^TA x = A^T b\]

Da cui si ottengono le cosiddette equazioni normali:

\[x = (A^TA)^{-1} A^T b\]

• Il risultato è il vettore \(x\) di coefficienti…
• …Che minimizza la somma dei quadrati dei residui

Intuitivamente: i coefficienti di una buona approssimazione

Metodo Lineare dei Minimi Quadrati

Quindi, possiamo impostare il problema come al solito

• Otteniamo una matrice A con più righe che colonne
• Otteniamo una colonna b dei termini noti

A questo punto, in teoria ci basta risolvere le equazioni normali:

x = (A' * A) \ (A * b)

In pratica, è ancora più semplice. Basta calcolare:

x = A \ b

• Matlab si accorge che A non è quadrata…
• …E risolve le equazioni normali (senza dirci niente)

Metodo Lineare dei Minimi Quadrati

Quindi, possiamo impostare il problema come al solito

• Otteniamo una matrice A con più righe che colonne
• Otteniamo una colonna b dei termini noti

A questo punto, in teoria ci basta risolvere le equazioni normali:

x = (A' * A) \ (A * b)

In pratica, è ancora più semplice. Basta calcolare:

x = A \ b

• Attenzione: La sintassi utilizzata è la stessa della divisione sinistra
• Può trarre in inganno (e.g. e se il sistema è davvero sovradeterminato?)

Metodo Lineare dei Minimi Quadrati

In realtà, il metodo funziona se il sistema di partenza è lineare

Questo è vero se la funzione approssimante è nella forma:

\[f(x) = \sum_{j=0}^{n-1} \alpha_j g_j(x)\]

• I.e. una somma pesata di “funzioni base” \(g_j(x)\)

In generale, la nostra matrice dei coefficienti avrà questa struttura:

\[\left(\begin{array}{cccc} g_{n-1}(x_0) & g_{n-2}(x_0) & \ldots & g_{0}(x_0) \\ g_{n-1}(x_1) & g_{n-2}(x_1) & \ldots & g_{0}(x_1) \\ \vdots & \vdots & \ddots & \vdots \\ g_{n-1}(x_{m-1}) & g_{n-2}(x_{m-1}) & \ldots & g_{0}(x_{m-1}) \\ \end{array}\right)\]

• Una colonna per funzione base, una riga per data point

Metodo Lineare dei Minimi Quadrati

Torniamo a noi

Nel nostro esempio, sappiamo che vale l’equazione di Antoine:

\[P^* = e^{\left(a - \frac{b}{T}\right)}\]

• Dove \(a\) e \(b\) sono parametri da determinare empiricamente
• \(a\) è un coefficiente lineare, ma \(b\) no
• Quindi l’equazione non è nella forma corretta!

Metodo Lineare dei Minimi Quadrati

Torniamo a noi

Nel nostro esempio, sappiamo che vale l’equazione di Antoine:

\[P^* = e^{\left(a - \frac{b}{T}\right)}\]

• Dove \(a\) e \(b\) sono parametri da determinare empiricamente
• \(a\) è un coefficiente lineare, ma \(b\) no
• Quindi l’equazione non è nella forma corretta!

Fortunatamente, è facile adattarla:

\[\begin{align} & \ln P^* = a - \frac{b}{T} & \longrightarrow & & \ln P^* = a + b \left(-\frac{1}{T}\right) \end{align}\]

• In questo modo sia \(a\) che \(b\) sono coefficienti lineari
• Quindi possiamo applicare il metodo dei minimi quadrati!

Metodo Lineare dei Minimi Quadrati

Impostiamo un sistema come per l’interpolazione esatta

Quindi avremo:

\[a + b \left(-\frac{1}{T_0}\right) = \ln P^*_0 \\ a + b \left(-\frac{1}{T_1}\right) = \ln P^*_1 \\ \ldots\]

In forma matriciale:

\[\underbrace{\left(\begin{array}{cc} 1 & -\frac{1}{T_0} \\ 1 & -\frac{1}{T_1} \\ \ldots & \ldots \end{array}\right)}_{A} \left(\begin{array}{c} a \\ b \end{array}\right) = \underbrace{\left(\begin{array}{c} \ln P^*_0 \\ \ln P^*_1 \\ \ldots \end{array}\right)}_{b}\]

Metodo Lineare dei Minimi Quadrati

Nel codice:

Costruiamo la matrice \(A\) e la colonna \(b\) dei termini noti:

T = T' % Vettore delle temperature come colonna
A = [T.^0, -1./T]
b = log(P)' % Logaritmo delle temp. come colonna

Risolviamo le equazioni normali:

x = A \ b

• x(1) corrisponde ad a, x(2) a b
• Possiamo definire la funzione e disegnarla, e.g.

f = @(t) exp(a - b ./ t)
plot(T, f(T))

Metodo Lineare dei Minimi Quadrati

Nel nostro caso otteniamo:

Esercizio: Comp. Reologico di un Fluido

Si possono dividere i fluidi in diverse categorie…

• …In base alla relazione tra la velocità di deformazione \(\dot{\gamma}\)…
• …E lo sforzo tangenziale \(\tau\)

Esercizio: Comp. Reologico di un Fluido

In particolare, distinguiamo:

Fluidi Newtoniani, per cui:

\[\tau = a \dot{\gamma}\]

Fluidi di Bingham, per cui:

\[\tau = a \dot{\gamma} + b\]

Fluidi pseudoplastici e dilatanti, per cui:

\[\tau = b \dot{\gamma}^a\]

• Se \(a > 1\) il fluido è dilatante
• Se \(a < 1\) il fluido è pseudoplastico

Esercizio: Comp. Reologico di un Fluido

Per due fluidi, sono stati misurati i valori di \(\dot{\gamma}\) e \(\tau\)

• I dati sono disponibili nei file fluid1.xlsx e fluid2.xlsx
• In Matlab si possono importare matrici da file Excel!
• Il file di script es_fluidbehavior.m contiene già il codice di lettura

Vogliamo determinare il tipo dei due fluidi. Per farlo:

Otterremo delle funzioni approssimanti \(f\) (via minimi quadrati)

• Verificheremo il valore del Mean Squared Error:

\[MSE = \frac{1}{n} \sum_{i = 1}^n (\tau_i - f(\dot{\gamma}_i))^2\]

È semplicemente la media dei residui al quadrato

• L’approssimazione con il minor MSE sarà tendenzialmente corretta

Esercizio: Comp. Reologico di un Fluido

Prima parte: verifica se il fluido sia di Bingham

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_Bingham(DGAMMA, TAU)

Che prenda come parametri:

• Un vettore riga DGAMMA con i valori noti di \(\dot{\gamma}\)
• Un vettore riga TAU con i valori noti di \(\tau\)

La funzione assume che il fluido sia di Bingham e restituisce:

• In MSE il Mean Squared Error
• Una funzione f che calcoli \(\tau = f(\dot{\gamma}) = a \dot{\gamma} + b\)

Esercizio: Comp. Reologico di un Fluido

Prima parte: verifica se il fluido sia di Bingham

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_Bingham(DGAMMA, TAU)

• Si esegua la funzione per i dati dei due fluidi
• Si utilizzi il parametro di uscita f per disegnare:
  • La funzione approssimante, tra min(DGAMMA) e max(DGAMMA)
  • I punti noti originali

È possibile farlo utilizzando la funzione:

function plot_results(DGAMMA, TAU, f, t)

• Già definita nello start-kit

Esercizio: Comp. Reologico di un Fluido

Seconda parte: verifica se il fluido sia Newtoniano

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_Newton(DGAMMA, TAU)

Che prenda come parametri:

• Un vettore riga DGAMMA con i valori noti di \(\dot{\gamma}\)
• Un vettore riga TAU con i valori noti di \(\tau\)

La funzione assume che il fluido sia Newtoniano e restituisce:

• In MSE il Mean Squared Error
• Una funzione f che calcoli \(\tau = f(\dot{\gamma}) = a \dot{\gamma}\)

Esercizio: Comp. Reologico di un Fluido

Seconda parte: verifica se il fluido sia Newtoniano

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_Newton(DGAMMA, TAU)

• Si esegua la funzione per i dati dei due fluidi
• Si utilizzi il parametro di uscita f per disegnare:
  • La funzione approssimante, tra min(DGAMMA) e max(DGAMMA)
  • I punti noti originali

È possibile farlo utilizzando la funzione:

function plot_results(DGAMMA, TAU, f, t)

• Già definita nello start-kit

Esercizio: Comp. Reologico di un Fluido

Terza parte: verifica se il fluido sia pseudoplastico o dilatante

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_PowerLaw(DGAMMA, TAU)

Che prenda come parametri:

• Un vettore riga DGAMMA con i valori noti di \(\dot{\gamma}\)
• Un vettore riga TAU con i valori noti di \(\tau\)

La funzione assume che il fluido segua la legge \(\tau = b \dot{\gamma}^a\) e restituisce:

• In MSE il Mean Squared Error
• Una funzione f che calcoli \(\tau = f(\dot{\gamma}) = b \dot{\gamma}^a\)

Esercizio: Comp. Reologico di un Fluido

Terza parte: verifica se il fluido sia pseudoplastico o dilatante

Nel file es_fluidbehavior.m nello start-kit, definite la funzione:

function [MSE, f] = check_PowerLaw(DGAMMA, TAU)

• Si esegua la funzione per i dati dei due fluidi
• Si utilizzi il parametro di uscita f per disegnare:
  • La funzione approssimante, tra min(DGAMMA) e max(DGAMMA)
  • I punti noti originali

È possibile farlo utilizzando la funzione:

function plot_results(DGAMMA, TAU, f, t)

• Già definita nello start-kit

Esercizio: Qualità del Vino

I due file winequality-red.xlsx e winequality-white.xlsx…

…Contengono informazioni su una serie di vini portoghesi

• L’ultima colonna contiene una stima della qualità (range 1-10)
• Le altre colonne contengono il valore di alcune proprietà fisico/chimiche

Formalmente, ogni data point è una coppia \((x_i, y_i)\) in cui:

• \(y_i\) è la qualità del vino
• \(x_i\) un vettore con le proprietà del vino (i.e. \(x_{i,0}\), \(x_{i,1}\), etc.)

Vogliamo ottenere una approssimazione \(f(x)\) della qualità \(y\)

• A differenza del caso precedente, \(f(x\)) in questo caso è multivariata…
• …Semplicemente perché \(x\) è un vettore!

Esercizio: Qualità del Vino

Come gestire funzioni multivariate? Nel solito modo!

Il metodo lineare dei minimi quadrati si può applicare a funzioni nella forma:

\[f(x) = \sum_{j=0}^{n-1} \alpha_j g_j(x)\]

• Ogni \(\alpha_j\) è un peso (da determinare), ogni \(g_j(x)\) è una funzione base

Una funzione basa può utilizzare solo alcune componenti del vettore

Una scelta molto frequente:

\[f(x) = \sum_{j=0}^{n-1} \alpha_j x_j\]

• I.e. \(f(x)\) è una combinazione lineare delle componenti di \(x\)

Esercizio: Qualità del Vino

A partire dal file es_wine.m:

Approssimare la qualità del vino (rosso e bianco) con una funzione \(f(x)\) tale che:

\[f(x) = \sum_{j=0}^{n-1} \alpha_j x_j\]

• Andranno determinati la matrice dei coefficienti e la colonna dei termini noti
• Poi si potranno determinare i pesi \(\alpha_j\) con il metodo dei minimi quadrati

Una volta determinati i pesi:

• Disegnare (con histogram) l’istogramma dell’errore di stima, i.e.:

\[err(x_i) = y_i - f(x_i)\]

• Disegnare (con bar) un grafico a barre che mostri i valori dei pesi \(\alpha\)
• …Le barre più grandi corrisponderanno alle proprietà più importanti!

Esercizio: Deformazione Meccanica

Sono date delle misure di deformazione di un pezzo meccanico

• Alle ascisse c’è la posizione originaria \(X\) di alcuni punti sul pezzo
• Alle ordinate c’è la posizione \(X2\) dopo la deformazione

Esercizio: Deformazione Meccanica

Le misure non sono molto precise (sono “rumorose”)

A partire dal file es_strain.m:

• Approssimare relazione \(X2 = F(X)\) con il metodo dei minimi quadrati
• Determinate (visivamente) il grado minimo necessario…
• …Per avere una approssimazione polinomiale adeguata

Una possibile definizione dello sforzo meccanico è data da:

\[S(X) = 1 - F'(X)\]

• Dove \(F'(X)\) è la derivata di \(F(X)\)

Si calcoli lo sforzo meccanico per \(X = 1.5\)

Esercizio: Velocità di un Fluido in Condotta

La velocità di un fluido in una condotta varia lungo la sezione

• Tendenzialmente, la velocità è maggiore al centro
• Il profilo di velocità è simile ad una parabola

Supponiamo di avere delle misurazioni di velocità

Le misurazioni sono nella forma \((x_i, y_i)\)

• \(x_i\) rappresenta la posizione sulla sezione della condotta
• \(y_i\) è il valore della velocità misurata
• Le posizioni degli estremi della condotta sono note

Vogliamo ricostruire l’intero profilo di velocità

Esercizio: Velocità di un Fluido in Condotta

Partite dal file es_waterspeed.m nello start-kit

Approssimate il profilo di velocità:

• Come funzione approssimante, utilizzate una parabola…
• …Quindi applicate il metodo dei minimi quadrati

Procedete poi a disegnare:

• Il valore del polinomio interpolante su tutta la condotta
• Utilizzando scatter disegnate anche le misurazioni

Provate quindi ad utilizzare una interpolazione polinomiale esatta

• Determinate il grado del polinomio necessario e disegnatelo
• Quale delle due curve vi sembra più realistica?

Esercizio: Prestazione di una Pompa Centrifuga

Il diagramma di prestazione di una pompa centrifuga:

• Mostra la relazione tra la portata \(Q\) (alle ascisse)…
• …E la prevalenza \(H\) (una grandezza analoga alla pressione)

Esercizio: Prestazione di una Pompa Centrifuga

Spesso il diagramma viene definito a partire da valori noti, e.g.:

Nel file es_pump.m dello start-kit:

• Approssimare la curva utilizzando il metodo dei minimi quadrati
• Determinate (visivamente) il grado minimo necessario…
• …Per avere una approssimazione polinomiale adeguata
• Disegnate il polinomio approssimante ed i punti originari
• Determinate il valore di \(H\) per \(Q = 90\)
• Determinate lo stesso valore con una approssimazione lin. a tratti