Esempio: Argine di un Fiume

Supponiamo di dover progettare l’argine di un fiume

• Alle ascisse abbiamo una posizione orizzontale $x$
• Alle ordinate abbiamo l’altezza $y$

Esempio: Argine di un Fiume

L’argine deve:

• Essere definito da una curva parabolica
• Toccare il fiume in una posizione nota $(x_0, y_0)$
• Raggiungere il livello della strada nel punto $(x_1, y_1)$
• Raggiungere l’altezza massima $y_{max}$…
• In un punto di coordinata $x_2$ non nota

Di fatto, si tratta di progettare/tracciare una curva

• Allora possiamo provare a risolverlo…
• …Impostando un insieme di equazioni

Equazioni per il Problema

Proviamo a formulare equazioni per le condizioni da rispettare:

• La curva deve passare per il punto $(x_0, y_0)$:

$$\alpha_2 x_0^2 + \alpha_1 x_0 + \alpha_0 = y_0$$

• La curva deve passare per il punto $(x_1, y_1)$:

$$\alpha_2 x_1^2 + \alpha_1 x_1 + \alpha_0 = y_1$$

• La curva raggiungerà il max nel punto $x_2$ in cui la derivata si annulla:

$$2 \alpha_2 x_2 + \alpha_1 = 0$$

• In corrispondenza di $x_2$, l’altezza dovrà essere $y_{max}$:

$$\alpha_2 x_2^2 + \alpha_1 x_2 + \alpha_0 = y_{max}$$

Equazioni per il Problema

Nel complesso, abbiamo:

$$\begin{align} \alpha_2 x_0^2 + \alpha_1 x_0 + \alpha_0 &= y_0 \\ \alpha_2 x_1^2 + \alpha_1 x_1 + \alpha_0 &= y_1 \\ 2 \alpha_2 x_2 + \alpha_1 &= 0 \\ \alpha_2 x_2^2 + \alpha_1 x_2 + \alpha_0 &= y_{max} \end{align}$$

• Le variabili sono $\alpha_2, \alpha_1, \alpha_0$, ma anche $x_2$
• Ci sono delle espressioni non lineari, i.e. $\alpha_2 x_2, \alpha_2 x_2^2, \alpha_1 x_2$

In sostanza, abbiamo un sistema di equazioni non lineari:

• Quindi, non possiamo impostare il problema in forma matriciale!

Equazioni Non Lineari in Matlab

Per risolvere equazioni non lineari Matlab offre due funzioni:

function [X, FVAL, FLAG] = fzero(F, X0)
function [X, FVAL, FLAG] = fsolve(F, X0)

Entrambe risolvono equazioni nella forma:

$$f(x) = 0$$

• In altre parole, cercano di trovare un valore di $x$…
• …Che azzeri la funzione $f(x)$
• Procedono in moto iterativo
• Una di esse (fsolve) potrebbe non convergere

Una soluzione dell’equazione è detta uno zero della funzione

Equazioni Non Lineari in Matlab

Per risolvere equazioni non lineari Matlab offre due funzioni:

function [X, FVAL, FLAG] = fzero(F, X0)
function [X, FVAL, FLAG] = fsolve(F, X0)

Per quanto riguarda i parametri, abbiamo che:

• X0 è la stima iniziale di $x$ da cui partire
• F è la funzione da azzerare (passata come parametro)
• F deve avere un singolo parametro, i.e.:

function z = F(X)

• Il parametro (i.e X) corrisponde alla variabile $x$

Equazioni Non Lineari in Matlab

Per risolvere equazioni non lineari Matlab offre due funzioni:

function [X, FVAL, FLAG] = fzero(F, X0)
function [X, FVAL, FLAG] = fsolve(F, X0)

Per quanto riguarda i valori restituiti abbiamo che:

• X è l’ultimo valore di $x$ visitato
  • Se tutto è andato bene, è la soluzione
• FVAL è il valore di F in corrispondenza di X
  • Se tutto è andato bene, sarà molto vicino a 0
• FLAG è un numero che indica come sono andate le cose
  • Se vale 1, l’algoritmo ha trovato uno zero
  • Altrimenti, qualcosa è andato storto (consultate help)

Equazioni Non Lineari in Matlab

Per risolvere equazioni non lineari Matlab offre due funzioni:

function [X, FVAL, FLAG] = fzero(F, X0)
function [X, FVAL, FLAG] = fsolve(F, X0)

Le due funzioni usano una combinazione di metodi numerici

• fzero è basata (grossomodo) sul metodo della bisezione
  • Cerca di trovare un secondo punto x1…
  • …Tale che F(X0) e F(X1) abbiamo segni opposti
  • Se ce la fa e se F è continua, allora converge sempre
• fsolve usa metodi completamente diversi (i.e. trust region)
  • Non è garantito che converga…
  • …Ma può essere più veloce di fzero

Equazioni Non Lineari in Matlab

Per risolvere equazioni non lineari Matlab offre due funzioni:

function [X, FVAL, FLAG] = fzero(F, X0)
function [X, FVAL, FLAG] = fsolve(F, X0)

Il valore di x0 è molto importante:

• Se è scelto male, può impedire la convergenza
  • Nel caso di fzero, perché non si riesce a trovare un x1 adeguato
  • Nel caso di fsolve, per i limiti del metodo stesso
• Se ci sono più soluzioni, può determinare quale sia restituita
  • Il più delle volte, sarà la soluzione più vicina ad x0…
  • …Ma in generale non è garantito

Di solito, il valore di x0 si sceglie per intuizione o per tentativi

Sistemi di Equazioni Non Lineari

Per risolvere un sistema di equazioni non lineari…

…Cerchiamo lo zero di una funzione vettoriale, ossia risolviamo:

$$\begin{align} & F(x) = 0 && \text{con } F: \mathbb{R}^n \mapsto \mathbb{R}^m \end{align}$$

• Dove $n$ ed $m$ sono il numero di variabili e di equazioni.

Per esempio:

$$\begin{align} x^3 &= e^y \\ \log x &= y+1 \end{align} \longleftrightarrow F((x,y)) = (0, 0)$$

• Con:

$$F((x,y)) = \left(\begin{array}{cc} x^3 - e^y, & \log x - y-1 \end{array}\right)$$

Sistemi di Equazioni Non Lineari

Per i sistemi di equazioni non lineari…

…Possiamo usare solo fsolve

function [X, FVAL, FLAG] = fsolve(F, X0)

In questo caso:

• La funzione F dovrà avere ancora un singolo parametro:

function Z = F(X)

• Questa volta, però, X sarà un vettore anziché uno scalare
• E lo stesso vale per il valore di ritorno Z

fsolve cercherà un vettore X per cui F(X) restituisca un vettore nullo

Argine di un Fiume: Equazioni

Torniamo al nostro problema

• Innanzitutto portiamo tutti i membri a sx del segno =:

$$\begin{align} \alpha_2 x_0^2 + \alpha_1 x_0 + \alpha_0 - y_0 &= 0\\ \alpha_2 x_1^2 + \alpha_1 x_1 + \alpha_0 - y_1 &= 0\\ 2 \alpha_2 x_2 + \alpha_1 &= 0\\ \alpha_2 x_2^2 + \alpha_1 x_2 + \alpha_0 - y_{max} &= 0 \end{align}$$

Le espressioni definiscono i termini di una funzione vettoriale:

• La funzione dovrà avere come input un singolo vettore (e.g. “$\bf x$”)
• Il vettore dovrà contenere tutte le variabili:

$${\bf x} = (\alpha_2, \alpha_1, \alpha_0, x_2)$$

• L’ordine può (come sempre) essere scelto liberamente

Argine di un Fiume: Implementazione

A questo punto possiamo codificare la funzione in Matlab:

function z = f(a2, a1, a0, x0, x1, y0, y1, x2, ymax)
    % La curva deve passare per (x0, y0)
    z(1) = a2*x0^2 + a1*x0 + a0 - y0;
    % La curva deve passare per (x1, y1)
    z(2) = a2*x1^2 + a1*x1 + a0 - y1;
    % La derivata si deve annullare in X(4)
    z(3) = 2*a2*x2 + a1;
    % in X(4) la curva deve valere ymax
    z(4) = a2*x2^2 + a1*x2 + a0 - ymax;
end

• La variabile restituita z è un vettore di 4 elementi

Argine di un Fiume: Implementazione

La nostra f è la funzione da azzerare

• Ma ha troppi argomenti!
• fsolve richiede una funzione con un singolo argomento

Risolviamo il problema con una funzione anonima:

x0 = 0;
x1 = 7;
y0 = 0;
y1 = 2;
ymax = 5;
% Espongo un solo parametro
fz = @(X) f(X(1),X(2),X(3),x0,x1,y0,y1,X(4),ymax);

• Il parametro da esporre è X, il vettore delle variabili

Implementazione

Per trovare uno zero a questo punto usiamo:

[psol, fval, flag] = fsolve(fz, x0)

• In questo caso x0 specificherà un valore per ogni variabile
• Per esempio, potremmo avere: x0 = [-1, 1, 0, 3]

Scegliere il valore $x_0$ per i sistemi di equazioni può essere complicato

• Ci sono molti valori da decidere…
• …E disegnare la funzione può essere complicato (troppe dimensioni)

Come linea guida:

• Usate l’intuizione! Specie se il problema ha senso fisico
• Se qualcosa non torna, fate diversi tentativi

Soluzione

Per il nostro problema, l’argine avrà questa forma:

• Il massimo è raggiunto per $x_2 = 3.9446$

Esercizio: Argine di un Fiume (2)

Consideriamo di nuovo l’esempio di costruzione di un argine

Come in precedenza:

• Alle ascisse abbiamo una posizione orizzontale $x$
• Alle ordinate abbiamo l’altezza $y$

Esercizio: Argine di un Fiume (2)

L’argine deve:

• Essere definito da una curva parabolica $f$
• Toccare il fiume in una posizione nota $(x_0, y_0)$
• Toccare la strada in un punto $(x_1, y_1)$…
• …Di cui non è nota la coordinata $x_1$
• Raggiungere l’altezza massima $y_{max}$
• Avere una superficie complessiva pari a $s$

La superficie sarà data da:

$$S = \int_{\underbrace{x_0}_{= 0}}^{x_1} f(x) \, dx$$

Esercizio: Argine di un Fiume (2)

Il file es_riverbank2.m contiene i dati del problema:

• Impostate il problema di definizione della curva…
• …Vedrete che otterrete un sistema non lineare
• Risolvetelo con gli strumenti messi a disposizione da Matlab…
• …Dopo aver introdotto una o più funzioni opportunamente definite
• Disegnate l’andamento della curva dell’argine

Suggerimento:

• Usate una funzione “normale” per definire le equazioni…
• …Ed una funzione anonima per esporre solo i parametri desiderati
• Non è obbligatorio farlo, ma verrà comodo in futuro

Esercizio: Equilibri per Shepherd

Si consideri il modello di Shepherd:

$$x^{(k+1)} = \frac{r x^{(k)}}{1 + \left(\frac{x^{(k)}}{N}\right)^2 }$$

È un modello tempo-discreto per l’evoluzione di una popolazione:

• $x^{(k)}$ è il numero di individui al passo $k$-mo
• $r$ è un tasso di crescita
• $N$ è il valore di popolazione per il cui $r$ dimezza

Esercizio: Equilibri per Shepherd

Si consideri il modello di Shepherd:

$$x^{(k+1)} = \frac{r x^{(k)}}{1 + \left(\frac{x^{(k)}}{N}\right)^2 }$$

Il file es_shepherd.m nello start-kit contiene un simulatore

• Vogliamo provare a determinare lo stato di equilibrio…
• …Risolvendo l’equazione non lineare:

$$x = \frac{r x}{1 + \left(\frac{x}{N}\right)^2 }$$

• In pratica, richiediamo che lo stato non vari (i.e $x^{(k+1)} = x^{(k)}$)

Esercizio: Equilibri per Shepherd

Estendete il codice in es_shepherd.m

Costruite opportunamente la funzione da azzerare

• Disegnate l’andamento della funzione da azzerare
• …Così da individuare visivamente la posizione degli zeri

Trovate lo zero con fzero o fsolve

• Verificate anche i valori restituiti di fval e flag

Provate a variare il valore di partenza $x_0$

• Cercate di ottenere due punti di equilibrio diversi

Esercizio: Equilibri per Crescita Logistica

Si consideri il modello di Crescita Logistica:

$$x^{(k+1)} = r x^{(k)} \left(1 - \frac{x^{(k)}}{N}\right)$$

Che può essere utilizzato per l’evoluzione di una popolazione

• $x^{(k)}$ è il numero individui al passo $k$-mo
• $r$ è un tasso di crescita
• $N$ è il massimo valore della popolazione sostenibile

Il file es_logi_eq.m nello start-kit contiene un simulatore

• Vogliamo provare a determinare lo stato di equilibrio…
• …Risolvendo una equazione non lineare

Esercizio: Equilibri per Crescita Logistica

Estendete il codice in es_logi_eq.m

Costruite opportunamente la funzione da azzerare

• Disegnate l’andamento della funzione da azzerare…
• …Così da individuare visivamente la posizione degli zeri

Trovate lo zero con fzero o fsolve

• Verificate anche i valori restituiti di fval e flag

Provate a variare il valore di partenza $x_0$

• Osservate se viene trovato uno zero diverso

Esercizio: Crescita Logistica, Preda-Predatore

Si consideri questo modello preda-predatore visto a lezione:

$$\begin{align} & H^{(t+1)} = \overbrace{r \left(1 - \frac{H^{(t)}}{k} \right) H^{(t)}}^{\text{crescita logistica}} - \overbrace{s\, H^{(t)} P^{(t)}}^{\text{prede eliminate}} \\ & P^{(t+1)} = \underbrace{u\, P^{(t)}}_{\text{calo in assenza di prede}} + v \underbrace{\left( s\, H^{(t)} P^{(t)} \right)}_{\text{prede eliminate}} \end{align}$$

• $H$ è il numero di prede e $P$ quello di predatori
• $k$ è la massima popolazione di prede sostenibile
• $s$ è la frazione di $H^{(t)}$ che un predatore può “mangiare”
• $u < 1$ è il ritmo di scomparsa dei predatori in assenza di prede
• $v$ è il “bonus riproduttivo” per ogni preda “mangiata”

Esercizio: Crescita Logistica, Preda-Predatore

Il file es_logi_pp_eq.m contiene un simulatore:

• Si estenda il codice così da determinare un punto di equilibrio…
• …Risolvendo un sistema di equazioni non lineari

Il punto di partenza $x_0$ sarà in questo caso un vettore:

• Conterrà il numero iniziale di prede e predatori
• Per tentativi, trovate un $x_0$ tale che la soluzione trovata…
• …coincida con lo stato stabile raggiunto dalla simulazione

Esercizio: Gas Isoentropico

In un fluido comprimibile isoentropico all’equilibrio…

…Pressione e temperature si distribuiscono secondo le equazioni:

$$\begin{align} & 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} } \\ & T = T_0 \left[ 1 + \left(1 - \frac{1}{\gamma}\right) \frac{M}{R T_0} g (z_0 - z) \right] \end{align}$$

• $P, T$ sono la pressione e temperatura a quota $z$
• $P_0, T_0$ sono la pressione e temperatura a quota $z_0$
• $R$ è la costante dei gas perfetti
• $M$ è la massa del gas
• $\gamma$ è un parametro che caratterizza il gas

Esercizio: Gas Isoentropico

Se sono noti $P_0, T_0, \gamma, M, z_0-z$, determinare $P$ e $T$ è facile:

• Perché le due equazioni sono esplicite in $P$ e $T$
• Vale a dire: $P$ e $T$ compaiono da sole a sx dell’uguale

Supponiamo invece di disporre di $P, T, P_0, T_0$ e $z_0 - z$

• A partire dai dati noti (e.g. misurati con qualche strumento)…
• …Vogliamo determinare le caratteristiche del gas (i.e $\gamma$ e $M$)
• Le due equazioni date sono implicite in $\gamma$ e $M$…
• …Perché $\gamma$ e $M$ non possono essere isolate a sx del segno “=”

Quindi dobbiamo risolvere un sistema di due equazioni non lineari

Esercizio: Gas Isoentropico

Il file es_isentropic.m contiene i dati del problema

Definite una funzione:

function z = f(p, T, p0, T0, gm, M, z0z, R, g)

• Che corrisponda al sistema di equazioni da risolvere…
• …Ma esponga tutti i parametri

Definite quindi una funzione anonima che esponga solo $\gamma$ e $M$

• Usatela per determinare $\gamma$ e $M$
• NOTA: Usando due funzioni (normale + anonima)
• …È molto veloce definire quali siano le variabili ed i parametri
• Le variabili sono quelle “esposte” nella funzione anonima!

Esercizio: Dodge-Metzner

Consideriamo un fluido di potenza (dilatante o pseudoplastico)

• Se il fluido è in moto turbolento, il fattore di attrito $f$ in una condotta
• …è definito dalla seguente relazione (Dodge-Metzner):

$$\frac{1}{\sqrt{f}} = \frac{4}{n^{0.75}} \log\left(Re_{pl} f^{1 - n/2}\right) - \frac{0.4}{n^{1.2}}$$

• Dove $n$ è un parametro che caratterizza il fluido
• $Re_{pl}$ è un numero di Reynolds modificato

Vogliamo utilizzare la relazione per calcolare il valore di $f$

Si supponga che i valori di $Re_{pl}$ e $n$ siano noti

Esercizio: Dodge-Metzner

Il file es_dodge_metzner.m contiene i dati del problema

• Definite nel file una funzione ausiliaria:

function z = dodge_metzner(f, n, Re_pl)

• Dovremo rendere la funzione utilizzabile in fzero o fsolve…
• …usando una funzione anonima per esporre solo alcuni parametri
• Si disegni l’andamento della funzione da azzerare
• Si determini il valore di $f$ con fzero o fsolve

In un secondo momento:

• Si assuma che $f$ valga $0.0020$…
• …E si determini il valore corrispondente di $n$

Esercizio: Buckingham-Reiner

Consideriamo un fluido di Bingham

• Se il fluido è in moto laminare, il fattore di attrito $f$ in una condotta
• …è definito dalla seguente relazione (Buckingham-Reiner):

$$f = \frac{64}{Re} \left(1 + \frac{1}{6}\frac{He}{Re} - \frac{64}{3} \frac{He^4}{f^3 Re^7} \right)$$

• Dove $He$ è il numero di Hedstrom
• E $Re$ è il numero di Reynolds

Vogliamo utilizzare la relazione per calcolare il valore di $f$

Si supponga che i valori di $Re$ e $He$ siano noti

Esercizio: Buckingham-Reiner

Il file es_buckingham_reiner.m contiene i dati del problema

• Definite nel file una funzione ausiliaria:

function z = buckingham_reiner(f, Re, He)

• Dovremo rendere la funzione utilizzabile in fzero o fsolve…
• …usando una funzione anonima per esporre solo alcuni parametri
• Si disegni l’andamento della funzione da azzerare
• Si determini il valore di $f$ con fzero o fsolve

Attenzione: Per $f = 0$ l’equazione perde senso fisico!

• Scegliete il range per il disegno di conseguenza
• Conviene usare più dei 100 punti che linspace usa di default