Velocità Media in una Condotta

Supponiamo di aver misurato la velocità di un fluido in condotta

• La condotta è larga 15 cm
• Abbiamo effettuato misurazioni per $x_0 = 1.5, x_1 = 7.8, x_2 = 14.1$
• I valori misurati sono $y_0 = 0.8, y_1 = 1.3, y_2 = 0.82$ (in m/s)

Velocità Media in una Condotta

Visualmente, abbiamo a disposizione i dati seguenti:

Velocità Media in una Condotta

Supponiamo di aver misurato la velocità di un fluido in condotta

• La condotta è larga 15 cm
• Abbiamo effettuato misurazioni per $x_0 = 1.5$, $x_1 = 7.8$, $x_2 = 14.1$
• I valori misurati sono $x_0 = 0.8$, $x_1 = 1.3$, $x_2 = 0.82$ (in m/s)

Vogliamo stimare la velocità media del fluido

Velocità Media

Innanzitutto dobbiamo definire cosa si intenda per velocità media:

• Se $f(x)$ rappresenta la velocità in funzione della posizione...
• ...Allora la velocità media può essere calcolata come

$$\overline{v} = \frac{1}{b-a} \int_{a}^{b} f(x) dx$$

• Dove $a$ e $b$ sono gli estremi dell'intervallo di calcolo

Focalizziamoci sul calcolo dell'integrale:

La nostra difficoltà è che $f(x)$ non è disponibile

...Però possiamo inferirla mediante interpolazione!

Integrazione via Interpolazione

Possiamo individuare il polinomio interpolante:

$$f(x) = \sum_{j = 0}^{n} \beta_j x^j$$

Poi otteniamo la funzione integrale per via analitica:

$$F(x) = \sum_{j = 0}^{n} \beta_j \frac{1}{j+1} x^{j+1} + C$$

E calcoliamo l'integrale come:

$$\int_a^b f(x) dx = F(b) - F(a)$$

Integrazione via Interpolazione

Visualmente, andiamo a calcolare questa area:

Integrazione via Interpolazione

Risolviamo l'esempio utilizzando Octave:

x = [1.50, 7.80, 14.10];
y = [0.80, 1.30,  0.82];
a = 0; b = 15.0;

n = length(x)-1; % grado del polinomio
p = polyfit(x, y, n) % polinonio approssimante
P = polyint(p) % funzione integrale
V = polyval(P, b) - polyval(P, a) % integrale
vm = V ./ (b-a) % velocità media

• La velocità media è: $1.0669$
• Per maggiori dettagli, si veda la registrazione della lezione

Integrazione via Interpolazione

Possiamo disegnare il risultato con:

xp = linspace(x_l, x_u, 200); % range
yp = polyval(p, xp); % valori del polinomio
hold on
plot(xp, yp, 'k', 'linewidth', 2) % polinomio
col_area = [0.9, 0.9, 0.9];
area(xp, yp, 'facecolor', col_area) % area
scatter(x, y, 8, 'r', 'o', 'filled') % punti
hold off
xlim([x_l, x_u])
ylim([y_l, y_u])
grid()

Elementi di Informatica e
Applicazioni Numeriche T

Formule di Newton-Cotes

Una Variante del Problema Iniziale

Come ci comportiamo se le misurazioni sono ancora da effettuare?

Una Variante del Problema Iniziale

In linea di principio, possiamo procedere come prima

• Misuriamo ("campioniamo") il valore della velocità in punti arbitrari
• Poi otteniamo ed integriamo il polinomio approssimante

Però conviene essere furbi!

• Campionando la velocità in punti opportuni...
• ...È possibile semplificare l'integrazione!

In particolare, è possibile ottenere delle formule di integrazione

• In particolare, noi vedremo le formule di Newton-Cotes

Una Variante del Problema Iniziale

• Lo stesso tipo di problema si presenta se la nostra $f(x)$ è nota...
• ...ma difficile da integrare per via analitica

Vediamo un esempio:

La pressione di un gas in condizioni isoentropiche è data da:

$$p = p_0 \left[ 1 + \left(1 - \frac{1}{\gamma}\right) \frac{M}{R T_0} g (z_0 - z) \right]^{\frac{\gamma}{\gamma - 1} }$$

• $p_0$ e $T_0$ sono la pressione e la temperature a quota $z_0$
• $M$ è la massa, $\gamma$ un parametro del gas
• $R$ è la costante dei gas perfetti

Una Variante del Problema Iniziale

Supponiamo di trovarci nella situazione:

• $M = 10\, \mathrm{Kg}$, $p_0 = 1\, \mathrm{atm}$, $T_0 = 295\, \mathrm{°K}$

Ci interessa calcolare la pressione media:

• Per $z_0 - z = -10\, {\rm m}$ (10 metri sopra $z_0$)
• Per valori di $\gamma$ tra $1.1$ e $2.0$

Integrare la funzione in $\gamma$ è abbastanza complesso:

$$f(\gamma) = \left[ 1 + \left(1 - \frac{1}{\gamma}\right) \frac{10}{295\, R} g (-10) \right]^{\frac{\gamma}{\gamma - 1} }$$

Una Variante del Problema Iniziale

...Anche se nel range considerato, la funzione è molto semplice!

Formule di Newton-Cotes

Le formule di Newton-Cotes si applicano ad integrali nella forma:

$$\int_a^b f(x) dx$$

E permettono di calcolarne il valore approssimato

• Richiedono solo la possibilità di valutare $f(x)$
• Si basano sull'idea di calcolare $f(x)$ per punti equispaziati
• Per i punti ottenuti, si deriva il polinomio di Lagrange
• Quindi si integra analiticamente il polinomio

Si ottengono così formule diverse a seconda del numero di punti

Derivazione delle Formule

• Le formule di Newton-Cotes sono piuttosto semplici
• La loro derivazione richiede però diversi passaggi

Innanzitutto, si normalizza l'intervallo di integrazione

$$[a, b] \longrightarrow [-1, 1]$$

• Questo passaggio semplifica considerevolmente i calcoli
• Per il resto la funzione resta invariata

Vediamo un esempio grafico...

Derivazione delle Formule: Normalizzazione

Partendo dalla funzione originale...

Derivazione delle Formule: Normalizzazione

Otteniamo una funzione identica, definita su $[-1, 1]$

Derivazione delle Formule: Normalizzazione

Formalmente, definiamo $x$ rispetto ad una nuova variabile $y$:

$$x = \phi(y) = \frac{a+b}{2} + \frac{b-a}{2} y$$

• La sostituzione assicura che:

$$\begin{align} y = -1 & \leftrightarrow x = a\\ y = 1 & \leftrightarrow x = b\\ y = 0 & \leftrightarrow x = (a+b)/2 \end{align}$$

• Invece di utilizzare $f$ direttamente, lavoreremo con:

$$g(y) = f(\phi(y))$$

Derivazione delle Formule: Normalizzazione

Per concludere, cambiamo la variabile di integrazione:

$$\int_a^b f(x)\,{\rm d}x = \int_{-1}^{1} f(\phi(y))\, \phi'(y)\, {\rm d}y = \frac{b-a}{2} \int_{-1}^{1} g(y)\,{\rm d}y$$

A questo punto procediamo come segue:

• Otteniamo $g(y)$ utilizzando il polinomio di Lagrange
• Otteniamo la funzione integrale $G(y)$
• Otteniamo una formula per $G(1) - G(-1)$

Il risultato dipende dal numero di punti di interpolazione

Regola del Rettangolo

Cominciamo usando un solo punto di interpolazione:

In particolare, scegliamo $y_0 = 0$ (punto intermedio)

• Otteniamo il polinomio di Lagrange:

$$g(y) = g(y_0) \frac{1}{1} = f\left(\frac{a+b}{2}\right)$$

Perché $g(0) = f(\phi(0))$ e $\phi(0) = (a+b)/2$

• In pratica, stiamo approssimando la funzione...
• ...con il suo valore nel punto intermedio

Regola del Rettangolo

A questo punto otteniamo la funzione integrale:

$$G(y) = f\left(\frac{a+b}{2}\right) y + C$$

A quindi abbiamo che:

$$\int_a^b f(x)\,{\rm d}x = \frac{b-a}{2} (G(1) - G(-1)) = \frac{b-a}{2} 2 f\left(\frac{a+b}{2}\right)$$

Da cui deriva la regola del rettangolo:

$$\int_a^b f(x)\,{\rm d}x \simeq (b-a)\, f\left(\frac{a+b}{2}\right)$$

Regola del Rettangolo

In pratica, calcoliamo l'area di questo rettangolo:

Regola del Rettangolo

Per l'esempio del profilo di velocità le cose vanno un po' peggio...

Regola del Trapezoide

Per migliorare l'accuratezza, possiamo campionare più punti:

In particolare, scegliamo $y_0 = -1$ e $y_1 = 1$

• Otteniamo il polinomio di Lagrange:

$$g(y) = g(y_0) \frac{y-y_1}{y_0-y_1} + g(y_1) \frac{y-y_0}{y_1-y_0}$$

$$g(y) = f(a) \frac{y-1}{-2} + f(b) \frac{y+1}{2}$$

$$g(y) = -\frac{f(a)}{2} y + \frac{f(a)}{2} + \frac{f(b)}{2} y + \frac{f(b)}{2}$$

Regola del Trapezoide

Da $g(y)$ otteniamo $G(y)$:

$$G(y) = -\frac{f(a)}{4} y^2 + \frac{f(a)}{2} y + \frac{f(b)}{4} y^2 + \frac{f(b)}{2} y$$

A quindi abbiamo che:

$$\int_a^b f(x)\,{\rm d}x = \frac{b-a}{2} (G(1) - G(-1)) = \frac{b-a}{2} 2 \frac{f(a) + f(b)}{2}$$

Da cui deriva la regola del trapezoide:

$$\int_a^b f(x)\,{\rm d}x \simeq \frac{1}{2}(b-a)(f(a) + f(b))$$

Regola del Trapezoide

In pratica, calcoliamo l'area di questo trapezio:

Regola del Trapezoide

Per il profilo di velocità, però, le cose vanno sempre malino...

Regola di Simpson

Proviamo a campionare un punto in più:

Scegliamo $y_0 = -1$, $y_1 = 0$ e $y_2 = 1$

• Otteniamo il polinomio di Lagrange:

$$\begin{align} g(y) & = g(y_0) \frac{y-y_1}{y_0-y_1}\frac{y-y_2}{y_0-y_2}\\ & + g(y_1) \frac{y-y_0}{y_1-y_0}\frac{y-y_2}{y_1-y_2}\\ & + g(y_2) \frac{y-y_0}{y_2-y_0}\frac{y-y_1}{y_2-y_1} \end{align}$$

Regola di Simpson

Proviamo a campionare un punto in più:

Scegliamo $y_0 = -1$, $y_1 = 0$ e $y_2 = 1$

• Otteniamo il polinomio di Lagrange:

$$\begin{align} g(y) & = f(a) \frac{y}{-1}\frac{y-1}{2}\\ & + f\left(\frac{a+b}{2}\right) \frac{y+1}{1}\frac{y-1}{-1}\\ & + f(b) \frac{y+1}{2}\frac{y}{1} \end{align}$$

Regola di Simpson

Proviamo a campionare un punto in più:

Scegliamo $y_0 = -1$, $y_1 = 0$ e $y_2 = 1$

• Otteniamo il polinomio di Lagrange:

$$\begin{align} g(y) & = - \frac{1}{2} f(a) (y^2 - y)\\ & - f\left(\frac{a+b}{2}\right) (y^2 - 1)\\ & + \frac{1}{2} f(b) (y^2 + y) \end{align}$$

Regola di Simpson

Da $g(y)$ otteniamo $G(y)$:

$$\begin{align} G(y) &= - \frac{1}{2} f(a) \left(\frac{y^3}{3} - \frac{y^2}{2} \right)\\ & - f\left(\frac{a+b}{2}\right) \left(\frac{y^3}{3} - y\right)\\ & + \frac{1}{2} f(b) \left(\frac{y^3}{3} + \frac{y^2}{2} \right) \end{align}$$

E quindi, con un po' di derivazioni:

$$\int_a^b f(x)\,{\rm d}x = \frac{b-a}{2} \left(\frac{1}{3}f(a) + \frac{4}{3}f\left(\frac{a+b}{2}\right) + \frac{1}{3} f(b) \right)$$

Regola di Simpson

Da $g(y)$ otteniamo $G(y)$:

$$\begin{align} G(y) &= - \frac{1}{2} f(a) \left(\frac{y^3}{3} - \frac{y^2}{2} \right)\\ & - f\left(\frac{a+b}{2}\right) \left(\frac{y^3}{3} - y\right)\\ & + \frac{1}{2} f(b) \left(\frac{y^3}{3} + \frac{y^2}{2} \right) \end{align}$$

Con qualche passo in più otteniamo la regola di Simpson:

$$\int_a^b f(x)\,{\rm d}x \simeq \frac{1}{6} (b-a) \left(f(a) + 4 f\left(\frac{a+b}{2}\right) + f(b) \right)$$

Regola di Simpson

La regola di Simpson calcola l'area di una parabola:

Regola di Simpson

In entrambi i nostri casi, fornisce una approssimazione accurata

Formule di Newton-Cotes

Riassumendo, abbiamo:

Formule di Newton-Cotes

Le formule di Newton-Cotes:

Funzionano solo per campionamenti equispaziati

• Se i punti di campionamento non lo sono...
• ...dobbiamo per forza costruire l'interpolazione separatamente

Richiedono solo l'abilità di valutare $f(x)$

• Dentro $f(x)$ si può "nascondere di tutto"!
• E.g. equazione non lineare, IPF...

Sono meno robuste ed accurate di metodi più moderni

• E.g. Gauss-Konrod, Clenshaw-Curtis...

Elementi di Informatica e
Applicazioni Numeriche T

Regole Composte

Gestione di Molti Campionamenti

Cosa fare se abbiamo molti campionamenti?

Gestione di Molti Campionamenti

Se applichiamo lo stesso metodo, potrebbe capitare questo!

Gestione di Molti Campionamenti

L'interpolazione esatta può divergere se il numeri di punti cresce

• Succede perché il grado del polinomio diventa troppo alto
• Ricordare il fenomeno Runge?

Conseguenza importante:

Che fare, allora se abbiamo molti punti?

Gestione di Molti Campionamenti

Una soluzione: usare il metodo dei minimi quadrati

Gestione di Molti Campionamenti

Problema: che grado usare per l'interpolazione?

• In alcuni casi è ovvio o quasi (e.g. velocità di fluidi)
• Però potrebbe non esserlo

Consideriamo questo esempio:

• Un serbatoio atmosferico contiene acqua
• Il serbatoio è alimentato naturalmente
• Il livello dell'acqua è stato misurato regolarmente per 31gg

Si vuole determinare il livello medio del serbatoio

Gestione di Molti Campionamenti

Supponiamo che le misurazioni siano queste:

Gestione di Molti Campionamenti

In questo caso il grado da utilizzare non è ovvio

Gestione di Molti Campionamenti

Peggio ancora: potrebbe essere necessario usare un grado alto

Regole Composte

Un approccio alternativo:

Primo passo:

• Dividiamo il range di integrazione in intervalli regolari...
• ...In modo che ogni intervallo contenga $n$ punti

Secondo passo:

• Per ogni intervallo interpoliamo con un polinomio di grado $n-1$
• Integriamo e sommiamo i risultati

Questo tipo di approccio si chiama regola composta

Regole Composte: Un Esempio

Per esempio, dividiamo in intervalli con 2 punti ciascuno

Regole Composte: Un Esempio

E poi integriamo delle approssimazioni lineari

Regole Composte: Un Esempio

Oppure, dividiamo in intervalli con 3 punti ciascuno

Regole Composte: Un Esempio

E poi integriamo delle approssimazioni di grado 2

Regole Composte: Un Esempio

Il valore dei due integrali è:

• Con l'approssimazione di grado 1: $70.343$
• Con l'approssimazione di grado 2: $70.517$

I valori della media sono quindi:

• Con l'approssimazione di grado 1: $70.343/31 = 2.269$
• Con l'approssimazione di grado 2: $70.517/31 = 2.274$

Alcune Osservazioni

• Sia $m$ il numero di campionamenti
• Sia $n$ il grado dell'approssimazione negli intervalli

Allora, per applicare una regola composta deve valere:

$$mod(m-1, n) = 0$$

Ossia, $n-1$ deve essere multiple di $n$

Conseguenze:

• L'approssimazione di grado $1$ si può sempre applicare
• Il grado $2$ si applica se $m$ è dispari

In ogni caso, nessuno ci vieta di "mescolare" più approssimazioni

Alcune Osservazioni

• Se i campionamenti sono fatti ad intervalli regolari...
• ...Allora si possono usare le formula di Newton-Cotes

Nell'esempio appena visto:

• Per il grado 1 si può usare la regola del trapezoide
• Per il grado 2 si può usare la regola di Simpson

Con i dati sperimentali, però, non capita spesso

Regole Composte per Funzioni Note

Se $f(x)$ è nota, ma difficile da integrare:

• Possiamo scegliere i punti da campionare a piacimento
• Le formule di Newton-Cotes si possono sempre applicare

Inoltre:

• Se l'approssimazione non è adeguata...
• ...Si possono campionare più punti

Molti metodi di integrazione avanzati sfruttano questa possibilità

• Quelli che lo fanno si chiamano "metodi adattativi"

Elementi di Informatica e
Applicazioni Numeriche T

Integrazione Numerica in Octave

Integrazione Numerica in Octave

Octave fornisce diverse funzioni per fare integrazione:

Se la funzione è nota:

% Quadratura gaussiana
[q, ier] = quad(f, a, b)
% Regola (composta) di Simpson adattativa
q = quadv(f, a, b) 
% Metodo di Clenshaw-Curtis adattativo
[q, err] = quadcc(f, a, b)

Ingressi (per dettagli usare help):

• $f$ la funzione da integrare (passata come valore)
• $a$, $b$ gli estremi dell'intervallo di integrazione

Integrazione Numerica in Octave

Octave fornisce diverse funzioni per fare integrazione:

Se la funzione è nota:

% Quadratura gaussiana
[q, ier] = quad(f, a, b)
% Regola (composta) di Simpson adattativa
q = quadv(f, a, b) 
% Metodo di Clenshaw-Curtis adattativo
[q, err] = quadcc(f, a, b)

Uscite (per dettagli usare help):

• $q$ il valore dell'integrale
• $ier$ codice di errore (0 = ok), $err$ stima dell'errore

Integrazione Numerica in Octave

Octave fornisce diverse funzioni per fare integrazione:

Se la funzione è nota:

% Quadratura gaussiana
[q, ier] = quad(f, a, b)
% Regola (composta) di Simpson adattativa
q = quadv(f, a, b) 
% Metodo di Clenshaw-Curtis adattativo
[q, err] = quadcc(f, a, b)

• quad: buono per funzioni con derivate continue
• quadv: relativamente simile ai metodi visti
• quadcc: il metodo più robusto

Integrazione Numerica in Octave

Octave fornisce diverse funzioni per fare integrazione:

Se abbiamo dei dati sperimentali:

q = trapz(x, y)

Integra utilizzando la regola (composta) dei trapezi

• Quando si lavora con dati sperimentali...
• ...Octave è un po' limitato

Per fortuna, i metodi che abbiamo visto sono abbastanza semplici:

• Se volete utilizzare approssimazioni diverse...
• ...Potete sempre implementarle direttamente

-

Elementi di Informatica e
Applicazioni Numeriche T

Qualche Esercizio di Esempio

Esercizio 1

Integrare sull'intervallo $[0,10]$ la funzione:

$$f(x) = \frac{x}{1 + |sin(x)|}$$

• Si disegni la funzione (per farsene una idea)
• Si provino ad utilizzare i metodi di Octave: quad, quadv, quadcc
• Poi si provi ad implementare la regola (composta) dei rettangoli
• Si verifichi cosa succede aumentando il numero di campioni

La soluzione è disponibile sul sito del corso

• Per dettagli, si rimanda alla registrazione della lezione

Esercizio 2

Si consideri il serbatoio dell'esempio della lezione

Si assuma che il serbatoio scarichi acqua attraverso una condotta:

• La condotta è lunga $15 m$
• Il diametro è di $5 cm$
• La scabrezza è di $1 mm$
• La costante di perdita di carico per l'imbocco è $k = 0.5$
• Lo sbocco della condotta è $1 m$ sotto il serbatoio

I dati sul livello del serbatoio sono nel file es2.csv

Si carichi la portata media nei 31 giorni considerati

Esercizio 2

La situazione del serbatoio in sintesi: