Esercizio: Equilibrio Termodinamico

Consideriamo una reazione chimica nella forma:

\[\nu_1 C_1 + \nu_2 C_2 + ... \longrightarrow \ldots \nu_{n-1} C_3 + \nu_{n} C_4\]

• Ogni \(C_i\) rappresenta una sostanza (reagente o prodotto)
• \(\nu_i\) è il corrispondente coefficiente stechiometrico
• I reagenti sono considerati con coefficiente negativo

Le quantità \(n_i\) delle sostanze all’equilibrio devono soddisfare:

\[\ln K = \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \prod_{i=1}^n n_i^{\nu_i}\]

• Dove \(P\) è la pressione e \(K\) è la costante di equilibrio termodinamico
• \(K\) dipende da pressione e temperatura
• \(n_{tot}\) è la somma di tutti gli \(n_i\)

Vedrete (avete visto?) questi argomenti nel corso di termodinamica

Esercizio: Equilibrio Termodinamico

La quantità \(n_i\) di ogni sostanza all’equilibrio…

…dipende dal “grado di avanzamento” \(\xi\) della reazione:

\[n_i = n_i^{(0)} + \xi \nu_i\]

• Dove \(n_i^{(0)}\) è la quantità iniziale della sostanza

Quindi, dobbiamo trovare uno \(\xi\) tale che valga:

\[\begin{align} & \ln K = \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \prod_{i=1}^n n_i^{\nu_i} \\ \text{con: } & n_i = n_i^{(0)} + \xi \nu_i \\ & n_{tot} = \sum_{i=1}^n n_i \end{align}\]

Si tratta quindi di risolvere una equazione non lineare

Esercizio: Equilibrio Termodinamico

Come caso specifico, consideriamo la reazione:

\[CH_4 + H_2O \longrightarrow CO + 3H_2\]

• A 25 atm e 800 K, la costante \(K\) vale 163
• Inizialmente abbiamo 4 moli di metano e 10 d’acqua

Trattiamo la nostra equazione come una funzione da azzerare:

\[\ln K - \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \prod_{i=1}^n n_i^{\nu_i} = 0\]

• Con \(n_i = n_i^{(0)} + \xi \nu_i\) e \(n_{tot} = \sum_{i=1}^n n_i\)

Per prima cosa, dovremo incapsularla in una funzione Matlab

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
end

I parametri sono:

• Il vettore n0 con le quantità iniziali
• Il grado di avanzamento xi
• Il vettore dei coefficienti stechiometrici nu
• La costante di equilibrio K
• La pressione P

Nota: usando dei vettori il codice potrà funzionare per qualsiasi reazione!

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
end

Corrisponde a:

\[n_i = n_i^{(0)} + \nu_i \xi\]

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
  ntot = sum(n);
end

Corrisponde a:

\[n_{tot} = \sum_{i=1}^n \nu_i\]

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
  ntot = sum(n);
  A = log(K);
end

Corrisponde a:

\[\underbrace{\ln K}_{A} - \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \prod_{i=1}^n n_i^{\nu_i} = 0\]

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
  ntot = sum(n);
  A = log(K);
  B = (P / ntot)^sum(nu);
end

Corrisponde a:

\[\ln K - \underbrace{\left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i}}_{B} \prod_{i=1}^n n_i^{\nu_i} = 0\]

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
  ntot = sum(n);
  A = log(K);
  B = (P / ntot)^sum(nu);
  C = prod(n.^nu);
end

Corrisponde a:

\[\ln K - \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \underbrace{\prod_{i=1}^n n_i^{\nu_i}}_{C} = 0\]

Esercizio: Equilibrio Termodinamico

Ci conviene fare i calcoli un passo per volta:

function z = thermoeq(n0, xi, nu, K, P)
  n = n0 + xi .* nu;
  ntot = sum(n);
  A = log(K);
  B = (P / ntot)^sum(nu);
  C = prod(n.^nu);
  z = A - B * C;
end

Corrisponde all’intera equazione:

\[\ln K - \left(\frac{P}{n_{tot}}\right)^{\sum_{i=1}^n \nu_i} \prod_{i=1}^n n_i^{\nu_i} = 0\]

Esercizio: Equilibrio Termodinamico

A partire dal file es_thermoeq.m

Codificare (come visto nelle slides precedenti) la funzione:

function z = thermoeq(n0, xi, nu, K, P)

• Definite una funzione anonima per esporre l’interfaccia richiesta da fzero
• Disegnate la funzione da azzerare
  • Attenzione: la nostra thermoeq assume che xi sia uno scalare
  • Per disegnare il suo valore occorrerà valutarla ripetutamente!
  • Cercate di capire quanti zeri vi siano
• Utilizzate fzero per individuare il punto di equilibrio
  • Scegliete il punto di partenza in modo da convergere allo zero corretto

Esempio: Lunghezza di Curve Parametriche

Supponiamo di voler calcolare la lunghezza di una curva ellittica

Per farlo, ci serve una descrizione formale della traiettoria

• Di solito, un traiettoria si descrive mediante una curva parametrica…
• …Cioè una funzione con input scalare ed output vettoriale:

\[F: \mathbb{R} \mapsto \mathbb{R}^n\]

In particolare, una ellissi è descritta da:

\[F(t) = \left(\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right) = \left(\begin{array}{c} a \cos t \\ b \sin t \end{array}\right)\]

• \(a\) e \(b\) sono le lunghezze dei due semi-assi
• L’unica variabile che compare è in questo caso \(t\)

Esempio: Lunghezza di Curve Parametriche

Il risultato può essere qualcosa di questo genere:

Esempio: Lunghezza di Curve Parametriche

Il risultato può essere qualcosa di questo genere:

• Per calcolare la lunghezza di una curva parametrica…
• …Possiamo immaginare di dividerla in segmenti infinitesimi

Esempio: Lunghezza di Curve Parametriche

Il risultato può essere qualcosa di questo genere:

• Ogni segmento infinitesimo corrisponde ad un vettore tangente
• L’equazione si ottiene derivando ogni componente di \(F(t)\)

\[F'(t) = \left(\begin{array}{c} f'_1(t) \\ f'_2(t) \end{array}\right) = \left(\begin{array}{c} -a \sin t \\ b \cos t \end{array}\right)\]

Esempio: Lunghezza di Curve Parametriche

Il risultato può essere qualcosa di questo genere:

La lunghezza di un vettore tangente è quindi data da:

\[\|F'(t)\| = \sqrt{{f'_1(t)}^2 + {f'_2(t)}^2}\]

Esempio: Lunghezza di Curve Parametriche

Il risultato può essere qualcosa di questo genere:

La lunghezza della curva si ottiene integrando quella del vett. tangente:

\[\int_{t_0}^{t_1} \|F'(t)\| dt = \int_{t_0}^{t_1} \sqrt{{f'_1(t)}^2 + {f'_2(t)^2}} dt\]

• Spesso si ottiene un integrale difficile da calcolare per via simbolica!

Metodi di Quadratura in Matlab

Matlab offre due funzioni principali per effettuare integrazione:

function Q = integral(F,XMIN,XMAX)
function Q = trapz(X,Y)

Funzionano in modo radicalmente diverso:

integral richiede una funzione F che abbia un singolo parametro:

function Y = F(X)

• L’intervallo di integrazione XMIN..XMAX viene diviso in sotto-intervalli
• Per ogni sotto-intervallo, F viene invocata per ottenere campioni
• L’integrale sui sotto-intervalli viene approssimato in base ai campioni
• Eventualmente, si ripete la suddivisione per aumentare la precisione

Metodi di Quadratura in Matlab

Matlab offre due funzioni principali per effetturare integrazioni:

function Q = integral(F,XMIN,XMAX)
function Q = trapz(X,Y)

Funzionano in modo radicalmente diverso:

trapz utilizza il metodo dei trapezi

• Si assume che la funzione da integrare sia stata già campionata
• I vettori X e Y contengono le coordinate \(x\) e \(y\) dei campioni
• Viene calcolata l’area dell’interpolazione lineare a tratti

La funzione trapz è particolarmente utile per dati sperimentali

• Non c’è una vera funzione da integrare, ma solo delle misurazioni!

Esempio: Lunghezza di Curve Parametriche

Nel caso della nostra ellissi, abbiamo:

\[L = \int_{0}^{2 \pi} \sqrt{(- a \sin t)^2 + (b \cos t)^2} dt\]

Per calcolare l’integrale, innanzitutto definiamo la funzione da integrare

• Se l’espressione è semplice, possiamo usare una funzione anonima:

dl = @(t) sqrt((-a.*sin(t)).^2 + (b.*cos(t)).^2)

• Altrimenti, definiamo un nuova funzione con function
• In entrambi i casi, usiamo gli operatori elemento per elemento…
• …Perché integral e trapz funzionano manipolando vettori

Esempio: Lunghezza di Curve Parametriche

Il prossimo passo dipende dal metodo di integrazione scelto

Se vogliamo usare integral, possiamo scrivere:

L = integral(dl, 0, 2*pi)

• La funzione Matlab si occupa del campionamento

Se vogliamo usare trapz, possiamo scrivere:

X = linspace(0, 2*pi) % Valori di t
Y = dl(X) % Lunghezze del vettore tangente
Q = trapz(X, Y)

• Il campionamento va fatto prima di invocare la funzione

Esempio: Dimensionamento di una Pista

Vogliamo progettare una pista ellittica

• Assumiamo che il semi-asse \(b\) sia fissato
• Il semi-asse \(a\) va invece deciso in modo che…
• …La pista abbia la stessa lunghezza di quella di Indianapolis

Stima di Parametri

Sappiamo che la lunghezza della pista è data da:

\[L(a) = \int_{0}^{2 \pi} \sqrt{(a \sin t)^2 + (b \cos t)^2} dt\]

• Solo \(a\) è variabile \(\Rightarrow\) la lunghezza è una funzione di \(a\), i.e. \(L(a)\)

Se \(L^*\) è la lunghezza desiderata, deve valere:

\[L(a) = L^*\]

Si tratta di una equazione non lineare in \(a\)!

• La cosa strana è che \(L(a)\) è calcolata via integrazione numerica
• …Ma se risolviamo l’eq. con metodi numerici, questo non importa
• …Perché ci basta poter calcolare la funzione da azzerare

Stima di Parametri

Dobbiamo risolvere:

\[L(a) - L^* = 0\]

Quindi potremmo scrivere:

b = ...;  % Valore di b
Lt = ...; % Lunghezza di Indianapolis
fz = @(a) curve_length(a, b) - Lt;
[asol, fval, flag] = fzero(fz, a0) % a0: stima iniziale

function L = curve_length(a, b)
    dl = @(t) sqrt((a.*sin(t)).^2 + (b.*cos(t)).^2);
    L = integral(dl, 0, 2*pi);
end

Stima di Parametri

Molti problemi di progettazione si possono affrontare così

• Supponiamo che il parametro da determinare si chiami \(x\)
• E di avere un vincolo su una grandezza \(y\) che dipende da \(x\)

A questo punto:

• Prima si trova il modo di calcolare \(y\) assumendo che \(x\) sia noto
• Poi si incapsula il metodo di calcolo in una funzione \(F(x)\)
• Quindi si risolve una equazione del tipo:

\[F(x) = y^*\]

• Dove \(y^*\) è il valore desiderato per \(y\)

Vedremo diversi esempi di qui alla fine del corso!

Esercizio: Bacino Idrico

Un piccolo bacino idrico è riempito artificialmente

Esercizio: Bacino Idrico

Un piccolo bacino idrico è riempito artificialmente

La portata d’acqua in ingresso (in \(m^3/h\)) è data da:

\[q(t) = a + b \sin\left(2\pi \frac{t}{24}\right) + c \sin\left(2\pi \frac{t}{13}\right)\]

• I coefficienti \(a, b, c\) sono noti

I dati del problema sono nel file es_flow.m

• La formula per la portata \(q(t)\) è già definita nella funzione:

function Q = flow_rate(t, a, b, c)

Esercizio: Bacino Idrico

Q1: Si definisca una funzione:

function W = intake(t0, t1, a, b, c)

• Che calcoli la quantità totale d’acqua che entra nel bacino…
• …Tra due estremi di tempo t0 e t1 (misurati in ore)
• Occorrerà calcolare (per via numerica) un integrale!

Si determini quanta acqua entra nel bacino in \(72\) ore

Q2: Quanto tempo ci vuole perché entrino \(200\, m^3\) d’acqua?

• Si tratta di un problema di stima di parametri
• Occorre determinare un valore di \(t_1\) che soddisfi la condizione

Esercizio: Bacino Idrico (2)

Un bacino idrico artificiale è alimentato naturalmente

La portata in ingresso (in \(m^3/h\)) è misurata ad intervalli regolari

• Un certo numero di misurazioni sono nel file flow.xlsx
• Il codice di lettura è disponibile nel file es_flow2.m

Esercizio: Bacino Idrico (2)

Q1: Si stimi la quantità d’acqua entrata nel periodo considerato

• Si effettui una integrazione a partire dai dati sperimentali

• Si utilizzi il metodo dei trapezi (funzione trapz)

Esercizio: Bacino Idrico (2)

Si assuma poi che parte della portata in ingresso sia dirottabile

In particolare, dirottiamo tutta la portata sopra un certo limite

• In un istante di tempo \(t\), la portata così limitata è data da:

\[\min(L, \mathit{pwl}(T, Q, t))\]

• \(\mathit{pwl}(T, Q, t)\) è l’approssimazione lineare a tratti della portata
  • \(T\) e \(Q\) sono i valori di tempi e portata nei dati sperimentali
  • È necessario usare una qualche forma di interpolazione perché…
  • …la vera portata è nota solo per gli istanti di tempo con misurazioni date
  • \(\mathit{pwl}(T, Q, t)\) può essere calcolata in Matlab con interp1
• \(L\) è il limite al di sopra del quale si incanala l’acqua altrove

Esercizio: Bacino Idrico (2)

Si assuma poi che parte della portata in ingresso sia dirottabile

Con \(L = 13\), la portata limitata appare così:

Esercizio: Bacino Idrico (2)

Q2: Si definisca la funzione:

function Wtot = flow_with_limit(T, Q, limit)

• Che dati i vettori T e Q con i tempi e le portate misurate…
• …E dato il valore limit del limite \(L\)…
• …Calcoli la quantità totale d’acqua arrivata nel periodo considerato

In pratica, la funzione deve calcolare:

\[\int_{\min(T)}^{\max(T)} \min(L, \mathit{pwl}(T, Q, t))\, dt\]

Si calcoli la quantità d’acqua arrivata, assumendo \(L = 13\)

Q3: Si determini il valore di \(L\) perché arrivino \(8,000\, m^3\) d’acqua

Esercizio: Posizionamento di una Pompa

Si deve posizionare una pompa su una condotta:

• La pompa deve servire due utenze
• La prima utenza si trova a \(1322\,m\) dalla condotta
• La seconda utenza si trova a \(2131\,m\) dalla condotta
• La seconda utenza è \(1\,Km\) a valle della prima

Le due utenze verranno servite:

• Costruendo delle condotte rettilinee…
• …Che connettono la pompa alle utenze

Vogliamo determinare la posizione ottimale della pompa

Esercizio: Posizionamento di una Pompa

La distanza totale della pompa dalle due utenze è data da:

\[dist(a) = \sqrt{(a - x_0)^2 + y_0^2} + \sqrt{(x_1 - a)^2 + y_1^2}\]

• Dove \(a\) è la posizione orizzontale della pompa
• \(x_0\), \(x_1\), \(y_0\), \(y_1\) sono le posizioni orizzontali e verticali delle utenze
• Si assume che per la condotta ha \(y = 0\)

Esercizio: Posizionamento di una Pompa

La posizione ottimale è quella che minimizza la distanza

• Disegnando l’andamento di \(dist(a)\) in funzione di \(a\)…
• …Possiamo notare che c’è un solo minimo

• Quindi un solo punto in cui la derivata si annulla

Esercizio: Posizionamento di una Pompa

Possiamo così calcolare il minimo risolvendo:

\[dist(a)' = 0\]

La derivata può essere calcolata in forma analitica:

\[dist(a)' = -\frac{x_0-a}{\sqrt{(x_0-a)^2 + y_0^2}} -\frac{x_1-a}{\sqrt{(x_1-a)^2 + y_1^2}}\]

In alternativa, possiamo approssimare la derivata per via numerica:

\[dist(a)' \simeq \frac{dist(a+\delta) - dist(a)}{\delta}\]

• Con \(\delta = a \sqrt{eps}\)
• \(eps\) è l’epsilon di macchina, accessibile mediante eps

Esercizio: Posizionamento di una Pompa

Il file es_pump_location.m contiene i dati del problema

Si calcoli il valore di \(a\) che minimizza la distanza, risolvendo \(dist'(a) = 0\)

• Q1: si utilizzi l’espressione analitica della derivata
  • Per calcolarla, definite una funzione:

  function dd = ddist(a, x0, x1, y0, y1)

• Q2: utilizzate l’approssimazione numerica
  • Per calcolarla, definite una funzione (f è la funzione da derivare):

  function dd = ddist2(f, a)

Confrontate i due risultati

Esercizio: Planata

Un aeroplanino di carta viene lanciato in orizzontale

Esercizio: Planata

Un aeroplanino di carta viene lanciato in orizzontale

La traiettoria nel tempo è descritta da una curva parametrica:

\[F(t) = \left(\begin{array}{c} v_x t \\ y_0 - \frac{1}{2}cgt^2 \end{array}\right)\]

• \(v_x\) è la velocità con cui viene lanciato
• \(y_0\) è la quota iniziale (in \(m\))
• \(g\) è l’accelerazione di gravità
• \(c\) è un coefficiente noto che tiene conto della resistenza dell’aria

I dati del problema sono disponibili nel file es_glide.m

Esercizio: Planata

Q1: Si definisca la funzione:

function L = glide_length(vx, g, c, t0, tf)

• Che calcoli la strada percorsa dall’aeroplanino…
• …Tra due istanti di tempo t0 e tf

Si determini la strada percorsa tra i due istanti specificati nel file

NOTA: strada percorsa = lunghezza della traiettoria percorsa

Q2: Si determini con che velocità \(v_x^*\) deve avvenire il lancio…

• …Perché la strada percorsa sia pari a \(15 m\)

Esercizio: Punto Intermedio

Si vuole piazzare un misuratore di velocità su un tratto di pista

Il tratto è definito da una parabola con estremi e coefficienti noti

Esercizio: Punto Intermedio

La parabola può essere vista come una curva parametrica

\[F(x) = \left(\begin{array}{c} x \\ a x^2 + b x + c \end{array}\right)\]

• Il parametro in questo caso è \(x\), che coincide con la prima coordinata
• Tutti i dati sono disponibili nel file es_halfway.m

Q1: si definisca una funzione:

function L = poly_length(p, x0, x1)

• Che, dati un polinomio p e due estremi per la coordinata \(x\)…
• …Calcoli lunghezza della curva polinomiale tra x0 e x1
• Ricordate che la derivata di un polinomio si può calcolare con polyder

Si utilizzi la funzione per calcolare la lunghezza del tratto di pista

Esercizio: Punto Intermedio

Q2: Il misuratore deve essere collocato a metà del tratto

• Si determinino le coordinate di un punto \((x', y')\)…
• …Che sia equidistante dai due estremi del tratto di pista

Il secondo quesito è un problema di stima di parametri

• Richiede di risolvere una equazione non lineare…
• …In cui la funzione da azzerare è calcolata via integrazione numerica

Le distanze sono uguali se la loro differenza è nulla, i.e.:

\[F((x, y)) = 0\]

Dove:

• \(F((x, y)) = D((x, y), (x_0, y_0)) - D((x, y), (x_1, y_1))\)
• Dove \(D((x, y), (x_0, y_0))\) è la distanza di \((x,y)\) da \((x_0, y_0)\), e così via