ODE e Sistemi Dinamici

Una equazione differenziale ordinaria (ODE) è definita da:

$$\dot{x} = f(t, x)$$

• Dove $t$ è una variabile scalare
• $x$ è una funzione di $t$, potenzialmente con valori vettoriali, i.e.:

$$x : \mathbb{R} \mapsto \mathbb{R}^n$$

• $\dot{x}$ è la derivata di $x$ in $t$

Informalmente, una eq. diff. ord. definisce un sistema dinamico

• $x$ è lo stato del sistema, che evolve nel tempo $t$
• $f$ specifica come la derivata di $x$ dipenda dallo stato e dal tempo
• Il tempo è una variabile continua $\Rightarrow$ il sistema è tempo-continuo

Riduzione al Primo Ordine

È sufficiente limitarsi alla derivata del primo ordine

• Se ci sono derivate di ordine maggiore, e.g. $\ddot{x}$…
• …Possono essere sempre ricondotte al primo ordine…
• …Introducendo un elemento di stato per ogni derivata intermedia

Per esempio:

$$\ddot{x} = a t^2$$

Introducendo $y = \dot{x}$, diventa:

$$\left(\begin{array}{c} \dot{x} \\ \dot{y} \end{array}\right) = \left(\begin{array}{c} y \\ a t^2 \end{array}\right)$$

• Caso tipico: $x =$ posizione, $y =$ velocità

Problemi ai Valori Iniziali

Quando si risolve una ODE, l’obiettivo è determinare $x$

Ma $x$ è una funzione! Quindi vogliamo determinarne l’andamento

• Per farlo, è necessario fare alcune assunzioni addizionali…
• …Per esempio, possiamo specificare un valore iniziale

In questo modo otteniamo un problema ai valori inziali:

$$\begin{align} & \dot{x} = f(t, x) \\ & x(0) = x_0 \end{align}$$

Intuitivamente, si tratta di simulare il sistema dinamico

• Poiché il sistema è tempo-continuo, $f$ non determina lo stato futuro…
• …Ma il modo in cui lo stato varia (i.e. la sua derivata)

ODE in Matlab

Matlab mette a disposizione diversi risolutori di ODE

Quello che utilizzeremo principalmente è accessibile mediante:

function [T, X] = ode45(f, tspan, x0)

Per quanto riguarda i parametri di ingresso/uscita:

• f è la funzione che definisce la ODE. La sua interfaccia deve essere:

function dx = f(t, x)

• t è il “tempo”
• x è lo stato corrente (un vettore)
• dx è la derivata dello stato corrente (un vettore colonna)

ODE in Matlab

Matlab mette a disposizione diversi risolutori di ODE

Quello che utilizzeremo principalmente è accessibile mediante:

function [T, X] = ode45(f, tspan, x0)

Per quanto riguarda i parametri di ingresso/uscita:

• tspan è un vettore con i due estremi dell’intervallo di simulazione
  • E.g. tspan = [0, 3600]
• x0 è lo stato iniziale
• T contiene tutti gli istanti di tempo visitati dall’algoritmo
• Ogni riga della matrice X contiene uno stato visitato
  • Esattamente come per la nostra vecchia simulate

Esempio: L’Auguste Piccard PX-8

L’attrattiva principale dell’Expo 1964 fu un sottomarino

• Faceva immersioni con passeggeri nel Lago di Ginevra

Esempio: L’Auguste Piccard PX-8

Per muoversi in verticale un sottomarino carica/scarica acqua

In questo modo esso varia la sua densità

• Se la densità è maggiore di $\rho$, il sottomarino “cade” nel fluido
• Se è minore, il sottomarino “cade” verso l’alto
• Se è uguale, il sottomarino si muove per inerzia

Per precisione, il sottomarino è soggetto a tre forze principali:

• La forza di gravità, che lo accelera verso il basso
• La forza di galleggiamento, che lo accelera verso l’alto
• L’attrito aerodinamico dell’acqua (trascinamento)

Esempio: L’Auguste Piccard PX-8

Sia $L$ il volume dell’acqua a bordo, allora abbiamo che:

La forza di gravità è data da (asse cartesiano orientato verso l’alto):

$$F_g = -g (m + \rho L)$$

• $m$ è la massa del sottomarino
• $g$ è l’accelerazione di gravità
• $\rho$ è la densità dell’acqua

La forza di galleggiamento è data da (principio di Archimede):

$$F_b = g \rho V$$

• $V$ è il volume del sottomarino

Esempio: L’Auguste Piccard PX-8

Sia $L$ il volume dell’acqua a bordo, allora abbiamo che:

• L’attrito aerodinamico dell’acqua (trascinamento) è dato da:

$$F_t = -\frac{1}{2} \rho A C_D\, v |v|$$

• $A$ è l’area della sezione
• $C_D$ è un coefficiente che dipende dalla forma
• $v$ è la velocità
  • Il prodotto $v |v|$ ha lo stesso segno di $v$…
  • …ed il valore assoluto di $v^2$

Esempio: L’Auguste Piccard PX-8

Supponiamo che l’acqua a bordo vari nel modo seguente:

• Il valore $L \simeq 52.2$ bilancia le forze di gravità e galleggiamento

Esempio: L’Auguste Piccard PX-8

In pratica, $L$ varia secondo una funzione lineare a tratti

• In Matlab la si può valutare con interp1
• I punti Tp e Lp che definiscono la funzione sono in es_submarine.m
• Il file contiene anche tutte le informazioni sul PX-8

Si desidera calcolare la quota del PX-8 nel tempo

• Questa è regolata dall’equazione differenziale:

$$\ddot{x} = \frac{1}{m + \rho L} (F_g + F_b + F_t)$$

• La quota iniziale $x_0$ è $-5\, m$
• In particolare, quale è la quota dopo 10, 20 e 30 minuti?

Esempio: L’Auguste Piccard PX-8

Dobbiamo risolvere un problema ai valori iniziali

Per farlo, dovremo invocare:

function [T, X] = ode45(f, tspan, x0)

Quindi dobbiamo determinare:

• Come descrivere lo stato $X$ del sistema
  • Una buona scelta: $X = (\text{posizione}, \text{velocità})$
• L’intervallo di simulazione
  • Una buona scelta: $[0, 1800]$ (30 minuti)
• La funzione f che ne calcola la derivata

Esempio: L’Auguste Piccard PX-8

Per definire f è utile procedere nel solito modo:

• Prima definiamo una funzione che calcoli la derivata dello stato…
• …A partire da tutti i dati del problema

function dx = dstate(t,x,g,V,M,rhoW,Cd,S,Tp,Lp)
    % x(1) = altitudine, x(2) = velocita'
    Fb = g * rhoW * V; % Galleggiamento
    L = interp1(Tp, Lp, t); % Carico d'acqua
    Fg = -(g*M + g*rhoW * L); % Gravità
    Ft = -0.5 * x(2) * abs(x(2)) * Cd * S * rhoW;
    % Calcolo le derivate
    dx(1) = x(2);
    dx(2) = (Fb+Fg+Ft) / (M + rhoW * L);
end

Esempio: L’Auguste Piccard PX-8

Per definire f è utile procedere nel solito modo:

• Quindi definiamo una funzione anonima…
• …Che esponga solo due parametri (i.e. $t$ ed $X$)

La funzione anonima potrà essere utilizzata per invocare ode45

dX = @(t, x) dstate(t,x,g,V,M,rhoW,Cd,S,Tp,Lp)';
x0 = [-5, 0]; % Stato iniziale
tspan = [0, 1800];
[t, X] = ode45(dX, tspan, x0);

• Notate che dX deve restituire un vettore colonna…
• …Per questa ragione trasponiamo il risultato di dstate
• Ricordate che la quota iniziale è $-5\, m$

Esempio: L’Auguste Piccard PX-8

A questo punto possiamo recuperare quota e velocità

x = X(:, 1); % Quota (prima colonna)
v = X(:, 2); % Velocità (seconda colonna)

Possiamo disegnare l’andamento della quota:

plot(t, x)

Esempio: L’Auguste Piccard PX-8

A questo punto possiamo recuperare quota e velocità

x = X(:, 1); % Quota (prima colonna)
v = X(:, 2); % Velocità (seconda colonna)

Ed ottenere la quota dopo 10, 20, 30 minuti:

x10 = interp1(t, x, 10 * 60)
x20 = interp1(t, x, 20 * 60)
x30 = interp1(t, x, 30 * 60)

• Usiamo interp1 perché la “funzione quota” x è definita per punti

Esercizio: BMW-i8

Una BMW i8 accelera a tavoletta su un rettilineo

Esercizio: BMW-i8

Una BMW i8 accelera a tavoletta su un rettilineo

Supponiamo che il motore eroghi una forza costante $F$

• L’auto ha un motore elettrico, così l’assunzione non è così irrealistica

Si oppone alla direzione del moto la forza di trascinamento:

$$F_t = -\frac{1}{2} \rho C_D S v |v|$$

• $\rho$ è la densità dell’aria, $v$ è la velocità
• $S$ è la superficie della sezione dell’auto
• $C_D$ è un coefficiente di trascinamento

Esercizio: BMW-i8

Quindi il sistema è definito dall’ODE:

$$\ddot{x} = \frac{1}{m} (F + F_t)$$

• Dove $m$ è la massa dell’auto

Tutti i dati sono nel file es_bmw.m

• Si determini l’andamento di velocità e posizione per 60 sec
• Si disegni l’andamento della velocità
• Si determini il tempo necessario per raggiungere i 27.8 m/s
• Si determini la strada percorsa in tale tempo
• Si determini la massima velocità raggiunta

Esercizio: Riscaldamento di una Stanza (2)

Vogliamo riscaldare una stanza con un convettore

• Il convettore riscalda l’aria, che sua volta riscalda i muri
• …Che disperdono parte del calore verso l’esterno
• La temperature del convettore e dell’esterno sono costanti
• L’aria della stanza ed i muri hanno capacità termiche non trascurabili

Possiamo modellare il sistema utilizzando un circuito RC equivalente:

• Cerchi = generatori di tensione (temperature)
• Zig-zag = resistenze (resistenze termiche)
• Linee parallele = capacitori (capacità termiche)

Esercizio: Riscaldamento di una Stanza (2)

Vogliamo ottenere un modello tempo-continuo

Ogni capacità si può modellare con una equazione differenziale:

$$\begin{align} & \dot{T}_a = \frac{1}{C_a} (w_{ca} - w_{aw}) && \dot{T}_w = \frac{1}{C_w} (w_{aw} - w_{wo}) \end{align}$$

Esercizio: Riscaldamento di una Stanza (2)

Vogliamo ottenere un modello tempo-continuo

I flussi invece si calcolano mediante equazioni lineari:

$$\begin{align} & w_{ca} = \frac{1}{R_{ca}} (T_c - T_a) \\ & w_{aw} = \frac{1}{R_{aw}} (T_a - T_w) \\ & w_{wo} = \frac{1}{R_{wo}} (T_w - T_o) \\ \end{align}$$

Esercizio: Riscaldamento di una Stanza (2)

Quindi nel complesso abbiamo le equazioni differenziali:

$$\begin{align} & \dot{T}_a = \frac{1}{C_a} (w_{ca} - w_{aw}) && \dot{T}_w = \frac{1}{C_w} (w_{aw} - w_{wo}) \\ \end{align}$$

Dove:

$$\begin{align} & w_{ca} = \frac{1}{R_{ca}} (T_c - T_a) \\ & w_{aw} = \frac{1}{R_{aw}} (T_a - T_w) \\ & w_{wo} = \frac{1}{R_{wo}} (T_w - T_o) \\ \end{align}$$

Esercizio: Riscaldamento di una Stanza (2)

Tutti dati del problema sono disponibile ne file es_heating2.m

• Assumendo che i valori iniziali di $T_a$ e $T_w$ siano noti…
• …vogliamo determinare l’andamento delle temperature per 2 ore
• Si disegni su un unico grafico l’andamento di $T_a$ e $T_w$
• Si determini la temperatura finale dell’aria e dei muri
• Si determini la temperature dell’aria e dei muri dopo un’ora
• Si determini il tempo necessario perché l’aria raggiunga i 20°C

Esercizio: Dalla Terra alla Luna

Nel libro di Jule Verne “Dalle Terra alla Luna”…

• …Il veicolo con i protagonisti viene “sparato” verso la Luna

Esercizio: Dalla Terra alla Luna

Durante il viaggio, la navicella è soggetta a forze gravitazionali

Esse sono regolate dalla legge di gravitazione di Newton:

$$F_{12} = -G \frac{m_1 m_2}{r_{12} \|r_{12}\|}$$

• $F_{12}$ è la forza esercitata dal corpo 2 sul corpo 1
• $G$ è la costante di gravitazione
• $m_1$ ed $m_2$ sono le masse del corpo 1 e 2
• $r_{12}$ è la distanza dal corpo 1 al corpo 2, i.e. $$r_{12} = x_1 - x_2$$
  • $x_1$ e $x_2$ sono le posizioni di 1 e 2
  • Funziona per scalari (la forma vettoriale è leggermente diversa)

Esercizio: Dalla Terra alla Luna

Si desidera modellare il moto della navicella

• Assumiamo per semplicità che la Terra e la Luna siano in fisse
• Quindi la nave viaggerà lungo una traiettoria verticale
• Il moto sarà regolato dell’equazione differenziale:

$$\ddot{x} = \frac{1}{m_s} (F_{se} + F_{sm})$$

• $m_s$ è la massa della navicella
• $F_{se}$ è l’attrazione esercitata dalla Terra sulla navicella
• $F_{sm}$ è l’attrazione esercitata dalla Luna sulla navicella

Esercizio: Dalla Terra alla Luna

Inoltre, assumiamo che:

• Il centro della Terra abbia quota 0
• La navicella parta da una quota $r_E$ (i.e. il raggio della Terra)

La navicella deve così raggiungere la quota $D-r_M$

• Con $D =$ distanza Terra-Luna e $r_M =$ raggio lunare

I dati del problema sono nel file es_moonshot.m

• Si disegni l’andamento della quota della navicella
  • Come velocità iniziale, si usino sia $11,000\, m/s$ che $11,100\, m/s$
  • Si verificano problemi di calcolo? Perché?
• Si determini la quota massima raggiunta: è superiore a $D-r_M$?

Esercizio: Serbatoi Comunicanti

Tre serbatoi comunicano attraverso condotte

Esercizio: Serbatoi Comunicanti

Tre serbatoi comunicano attraverso condotte

• Il problema può essere modellato con un circuito RC
• Serbatoi = capacitori, resistenze = condotte
• Flussi = correnti, Differenza di pressione = differenza di tensione
• Il modello è approssimativo, ma a noi basterà

Esercizio: Serbatoi Comunicanti

Il sistema è descritto dalle equazioni differenziali:

$$\begin{align} & \dot{P}_1 = \frac{1}{C_1} (q_{31} - q_{12}) && \dot{P}_2 = \frac{1}{C_2} (q_{12} - q_{23}) \\ & \dot{P}_2 = \frac{1}{C_3} (q_{23} - q_{31}) \end{align}$$

Con:

$$\begin{align} & q_{12} = \frac{1}{R_{12}} (P_1 - P_2) && q_{23} = \frac{1}{R_{23}} (P_2 - P_3) \\ & q_{31} = \frac{1}{R_{31}} (P_3 - P_1) \end{align}$$

Esercizio: Serbatoi Comunicanti

Tutti i dati del problema sono nel file es_tubes.m

• Assumendo che tutte le pressioni iniziali siano note…
• …Si determini e disegni l’andamento delle tre pressioni nel tempo
• Si determini la pressione nel serbatoio 1 dopo 600 minuti
• Dopo quanto tempo la differenza tra $P_1$ e $P_2$ è minore di $10^{-1}$?