Partiamo da un Problema

Partiamo da un problema semplice, ma realistico:

Partiamo da un Problema

Partiamo da un problema semplice, ma realistico:

• Una palla da baseball viene lanciata in aria...
• ...In direzione verticale, alla velocità di 15 m/s

Alcune domande interessanti:

• Come varia la posizione della palla nel tempo?
• Quale è (approssimativamente) l'altezza massima?
• Quando (approssimativamente) viene raggiunta?

Modellare il Problema

Innanzitutto, occorre modellare il problema:

• La palla si muove in verticale: sia $x$ al sua posizione
• La posizione iniziale è $0$ e la velocità iniziale $15\, (m/s)$
• La palla è soggetta alla forza di gravità

Modellare il Problema

Possiamo descrivere il sistema con un'equazione differenziale

In particolare, la palla è soggetta ad una accelerazione costante:

$$\begin{align} & \ddot{x} = -g \end{align}$$

• La posizione $x$ è in realta una funzione del tempo, i.e. $x(t)$
• Per ricordarlo, indichiamo le sue derivate aggiungendo dei "punti"
• E.g. $\dot{x}$ è la velocità, $\ddot{x}$ l'accelerazione

Sappiamo inoltre che:

$$\begin{align} & x(0) = 0\\ & \dot{x}(0) = 15 \end{align}$$

Soluzione (per Via Analitica)

Il problema è semplice e può essere risolto per via analitica

A partire da:

$$\begin{align} & \frac{d \dot{x}}{dt} = -g \end{align}$$

Possiamo ottenere:

$$\begin{align} d \dot{x} & = -g\, dt\\ \int d \dot{x} & = \int -g\, dt\\ \dot{x} & = -gt + C_0 \end{align}$$

Soluzione (per Via Analitica)

Combinando le due equazioni:

$$\begin{align} & \dot{x} = -gt + C_0 & & \dot{x}(0) = 15 \end{align}$$

• Determiniamo che $C_0 = 15$

Con la stessa tecnica, possiamo derivare:

$$\begin{align} \frac{dx}{dt} & = -gt + 15\\ \int dx & = \int -gt + 15\, dt\\ x & = - \frac{1}{2} gt^2 + 15 t + C_1 \end{align}$$

• Poiché sappiamo che $x(0) = 0$, abbiamo $C_1 = 0$

Soluzione di una Equazione Differenziale

Alla fine del procedimento, otteniamo quindi:

$$\begin{align} x(t) = - \frac{1}{2} gt^2 + 15 t \end{align}$$

• È la definizione di $x$, in funzione del tempo
• Rappresenta la soluzione dell'equazione originale

Intuitivamente: definisce come alcune grandezze variano nel tempo

Posizione e Velocità della Palla

Possiamo disegnare la funzione su un piano cartesiano!

Posizione e Velocità della Palla

Guardando il grafico, possiamo osservare che:

• La posizione della palla ha un andamento parabolico (nel tempo)
• È un caso classico di moto uniformemente accelerato!

Inoltre:

• La massima altezza è di circa 12 m
• Viene raggiunta dopo circa 1 secondo e mezzo

Si possono ottenere i valori esatti, ma per ora non lo facciamo

Invece, ci concentriamo sulla formalizzazione del problema

Equazioni Differenziali Ordinarie

Noi ci focalizziamo su equazioni differenziali ordinarie

Si tratta di equazioni differenziali nella forma:

$$x^{(n)} = f(t, x, \dot{x}, \ddot{x}, \ldots, x^{(n-1)})$$

Dove:

• $x$ è una funzione della variabile $t$
• $\dot{x}, \ddot{x}\ldots x^{(i)}$ sono le derivate di $x$
• $n$ definisce l'ordine dell'equazione

Comprendono molte equazioni di interesse pratico

• Caso più generale: equazioni differenziali algebriche

Soluzione Numerica di una EDO

Vedremo come risolvere un EDO utilizzando metodi numerici

Conseguenza:

• La soluzione di una equazione differenziale è una funzione $x(t)$
• Quindi una soluzione numerica sarà una sequenza di valori
• Ogni elemento corrisponde al valore di $x(t)$ per un certo $t$

Altra conseguenza:

• Per le EDO a volte viene impiegata una forma implicita

$$f(t, x, \dot{x}, \ddot{x}, \ldots, x^{(n)}) = 0$$

• Quando usiamo metodi numerici, è equivalente alla forma esplicita

Riduzione al Primo Grado

Di fatto, possiamo limitarci ad EDO di primo ordine

Quindi alla forma (esplicita):

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

Il "trucco" è che $x$ può essere una funzione vettoriale

• Possiamo eliminare le derivate di ordine più alto...
• ...Introducendo nuove variabili (e.g. $v$ per $\dot{x}$)
• ...E nuove equazioni nella definizione di $f$ (e.g. $v = \dot{x}$)

Si capisce meglio con un esempio...

Riduzione al Primo Grado

Per il nostro problema, abbiamo una equazioni di secondo ordine:

$$\begin{align} & \ddot{x} = -g \end{align}$$

Introduciamo una variabile per $\dot{x}$ e l'equazione che la definisce:

$$\begin{align} & \dot{x} = v\\ & \dot{v} = -g \end{align}$$

Formalmente, abbiamo $f:\mathbb{R}^2 \rightarrow \mathbb{R}^2$

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

Le due componenti di $f$ sono:

$$\begin{align} & f_1((x, v), t) = v & & f_2((x, v), t) = -g \end{align}$$

Problema ai Valori Iniziali

Per risolvere una EDO non è sufficiente la definizione:

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

Occorre specificare delle informazioni addizionali!

• Molto spesso, si specifica il valore di $x$ ad un istante $t_0$
• Si chiama "valore iniziale"

Conseguenza: chiamiamo problema ai valori iniziali un sistema:

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

• Un problema ai valori iniziali ammette un'unica soluzione

EDO in Octave

Octave permette di risolvere EDO di primo ordine

In particolare, si utilizza la funzione:

x = lsode(f, x0, t)

Dove "ODE" sta per Ordinary Differential Equations e:

• f è la funzione che definisce l'equazione
• x0 è il valore di $x(t_0)$
• t è il vettore di valori di $t$ per cui si vuole conoscere $x(t)$
• x è il vettore con i valori di $x(t)$

Per altre informazioni, usate help!

Soluzione dell'Esempio in Octave

Il "cuore" della soluzione dell'esempio in Octave è:

function y = f(x, t)
    % x(1) = posizione, x(2) = velocità
    dx = x(2);
    dv = -g;
    y = [dx, dv];
end

x0 = [0, 15]; % Valori iniziali
t = linspace(0, 4); % Primi 4 secondi
sol = lsode(@f, x0, t); % soluzione (matrice [x, v])
plot(t, sol(:, 1)) % stampo i valori di x

La soluzione sarà anche disponibile sul sito del corso

Elementi di Informatica e
Applicazioni Numeriche T

Motivazioni per i Metodi Numerici

Un Secondo Problema

Supponiamo di dover verificare la progettazione di montagne russe

Un Secondo Problema

Ci interessa il tratto in ingresso ad un "giro della morte"

Un Secondo Problema

Ci interessa il tratto in ingresso ad un "giro della morte"

• Il tratto ha una forma parabolica, data dall'equazione:

$$y = \frac{b}{a^2} x^2 - \frac{2b}{a} x + b$$

• $b = 20\, m$ è la quota di partenza del carrello
• $a = 15\, m$ è il punto di ingresso al giro della morte

Vogliamo determinare la velocità del carrello nel punto $a$

• Se non è abbastanza alta sono guai!
• HP: il carrello ha velocità nulla all'inizio della discesa

Un Secondo Problema

Il carrello è soggetto alla forza di gravità. In sintesi, abbiamo:

Modellare il Problema

Come sempre, la prima cosa da fare è modellare il problema

• In particolare, ci interessa la velocità tangenziale
• Dobbiamo quindi capire come sia influenzata da $g$
• La relazione dipenderà dalla posizione sulla curva

Intuitivamente:

• Se riusciamo a modellare queste relazioni come EDO
• Possiamo calcolare come velocità e posizione cambino nel tempo

A questo punto diventa possibile risolvere il problema

Accelerazione Tangenziale

Consideriamo un punto sulla discesa di ingresso

• Per semplicità, centriamo gli assi sul punto considerato

Accelerazione Tangenziale

Ci interessa l'accelerazione tangenziale

• Questa è data dalla proiezione del vettore sulla tangente

Accelerazione Tangenziale

Sia ${\bf g}$ il vettore accelerazione e ${\bf d}$ il vettore direzione della tangente

• Allora il vettore proiezione ha direzione ${\bf d}$ ed intensità ${\bf g}^T {\bf d}$

Accelerazione Tangenziale

Per avere la proiezione, ci serve quindi il vettore direzione

Algebricamente, tutti i punti sulla tangente hanno la forma:

$$\left(\begin{array}{c} x \\ y \end{array}\right) = \left(\begin{array}{c} x \\ \phi' x \end{array}\right)$$

Dove:

• $\phi(x)$ è la funzione che definisce la curva
• $\phi'$ è la derivata di $\phi(x)$ nel punto considerato

Per il vettore direzione deve valere:

$$\sqrt{x^2 + y^2} = 1$$

Accelerazione Tangenziale

Combinando le due equazioni otteniamo:

$$\sqrt{x^2 + (\phi' x)^2} = 1$$

E con un po' di derivazioni (semplici) deduciamo che:

$${\bf d} = \left(\begin{array}{c} \frac{1}{\sqrt{1 + \phi'^2}} \\ \frac{\phi'}{\sqrt{1 + \phi'^2}} \end{array}\right)$$

Poiché ${\bf g}^T = (0, -g)$, possiamo ottenere l'accelerazione tangenziale:

$$\dot{v} = {\bf g}^T {\bf d} = -g \frac{\phi'}{\sqrt{1 + \phi'^2}}$$

Accelerazione e Posizione

In sintesi l'accelerazione tangenziale è data da:

$$\dot{v} = {\bf g}^T {\bf d} = -g \frac{\phi'}{\sqrt{1 + \phi'^2}}$$

• In realtà, $\dot{v}$ è solo l'intensità del vettore
• La direzione del vettore accelerazione è data da ${\bf d}$

L'EDO caratterizza il comportamento della velocità tangenziale

• Risolvendo l'EDO, potremmo ottenere la funzione $v$
• Il problema è che $\phi'$ dipende dalla posizione $x$

Ci occorre un'altra equazione per caratterizzare la funzione $x$

Velocità e Posizione

Supponiamo che $v$ sia nota (l'abbiamo già caratterizzata)

• Si può ottenere la componente $x$ della velocità con una proiezione
• Per farlo, ci serve la velocità tangenziale in forma vettoriale

Accelerazione e Posizione

La velocità tangenziale in forma di vettore è data da:

$${\bf v} = v\, {\bf d} = v \left(\begin{array}{c} \frac{1}{\sqrt{1 + \phi'^2}} \\ \frac{\phi'}{\sqrt{1 + \phi'^2}} \end{array}\right)$$

• La direzione della velocità tangenziale è ovviamente ${\bf d}$

La componente $x$ della velocità tangenziale è quindi:

$$\dot{x} = {\bf v}^T \left(\begin{array}{c} 1 \\ 0 \end{array}\right) = v \frac{1}{\sqrt{1 + \phi'^2}}$$

• È semplicemente il primo termine del vettore

Un Modello per il Secondo Problema

Possiamo così ottenere una modello completo per il problema:

$$\begin{align} & \dot{x} = v \frac{1}{\sqrt{1 + \phi'(x)^2}} \\ & \dot{v} = -g \frac{\phi'(x)}{\sqrt{1 + \phi'(x)^2}} \end{align}$$

Con i valori iniziali:

$$\begin{align} & x(0) = 0 \\ & v(0) = 0 \end{align}$$

È una EDO definita da una funzione vettoriale!

Un Modello per il Secondo Problema

Un altro punto di vista per il modello:

Definiamo lo stato del sistema con due scalari:

• La velocità tangenziale $v$
• La posizione $x$
• Questa è la decisione di modellazione più importante

Determiniamo le due equazioni che definiscono $\dot{v}$ e $\dot{x}$

• Il procedimento è lo stesso che abbiamo visto
• Le due equazioni corrispondono alle due componenti di $f$

Risolvere il Secondo Problema

Potremmo tentare di risolvere il problema per via analitica

Le due EDO sono però molto difficili da integrare:

$$\begin{align} & \dot{x} = v \frac{1}{\sqrt{1 + \phi'(x)^2}} & & \dot{v} = -g \frac{\phi'(x)}{\sqrt{1 + \phi'(x)^2}} \end{align}$$

Dove, nel nostro caso (traiettoria parabolica):

$$\begin{align} & \phi(x) = \frac{b}{a^2} x^2 - \frac{2b}{a} x + b \\ & \phi'(x) = \frac{2b}{a^2} x - \frac{2b}{a} \end{align}$$

Risolvere il Secondo Problema

Se utilizziamo un metodo numerico, invece:

• Dobbiamo solo poter calcolare le due derivate $\dot{v}, \dot{x}$
• Possiamo farlo con (soluzione completa sul sito del corso):

function y = f(x, t)
    % x(1) = posizione, x(2) = velocità tangenziale

    phip = @(x) ((2.*b)./a.^2 .* x + - (2.* b)./a);
    fp = phip(x(1)); % derivata di phi in x

    dx =x(2) .* (1 ./ sqrt(1 + fp.^2));
    dv = -g .* (fp ./ sqrt(1 + fp.^2));
    y = [dx, dv];
end

Risolvere il Secondo Problema

La soluzione per i primi 3 secondi è la seguente:

Elementi di Informatica e
Applicazioni Numeriche T

Metodo Numerici per
Equazioni Differenziali Ordinarie

Soluzione Numerica di una EDO

Data una EDO del primo ordine in forma esplicita:

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

Ed un valore iniziale:

$$x(t_0) = x_0$$

Come possiamo risolvere il problema ai valori iniziali?

• Soluzione = una sequenza di valori di $x$
• Iniziamo con un metodo molto semplice
• Idea principale: rimpiazzare $\dot{x}$ con il rapporto incrementale

Metodo di Eulero Diretto

Metodo di Eulero (Diretto)

Rimpiazzando $\dot{x}$ con il rapporto incrementale otteniamo:

$$\frac{x(t + \delta) - x(t)}{\delta} \simeq \dot{x}(t) = f(x(t), t)$$

• Dove le dipendenza di $\dot{x}$ e $x$ da $t$ sono esplicitate

Con qualche piccola manipolazione otteniamo:

$$x(t + \delta) = x(t) + \delta\, f(x(t), t)$$

• Ossia una formula per calcolare il prossimo valore di $x$
• A partire da $x(t_0)$ possiamo calcolare $x(t)$ nei punti desiderati

Metodo di Eulero Diretto

Metodo di Eulero (Diretto)

Rimpiazzando $\dot{x}$ con il rapporto incrementale otteniamo:

$$\frac{x(t + \delta) - x(t)}{\delta} \simeq \dot{x}(t) = f(x(t), t)$$

• Dove le dipendenza di $\dot{x}$ e $x$ da $t$ sono esplicitate

La formula somiglia un po' a quella dell'IPF

$$x(t + \delta) = x(t) + \delta\, f(x(t), t)$$

• In questo caso però il sistema dinamico è tempo-continuo
• E ci interessa l'andamento, invece che lo stato stazionario

Metodo di Eulero Diretto

Un possibile codice per il metodo:

function x = euler_direct(f, x0, t)
    x(1, :) = x0; % prima riga (x è un vettore riga)
    for ii = 1:length(t)-1
        xc = x(ii, :); % x corrente
        tc = t(ii); % tempo corrente
        x(ii+1, :) = xc + f(xc, tc) .* (t(ii+1) - tc);
    end
end

• Alla fine, ogni riga di x contiene i valori desiderati
• Nel nostro caso, i valori di $x(t)$ e $v(t)$

Risultati, per il Nostro Esempio

Linea nera: lsode, linea rossa: metodo di Eulero, 50 punti

Metodo di Eulero: Accuratezza e Robustezza

Il metodo di Eulero non è particolarmente robusto:

• Usa per una approssimazione lineare
• Se $f$ è fortemente non lineare, l'approssimazione è cattiva

Il metodo di Eulero non è particolarmente preciso:

• Ogni passo si basa sui risultati precedenti...
• ...Quindi gli errori si accumulano!

Come gestire i due inconvenienti?

• Per la precisione: consideriamo punti più ravvicinati
• Per la robustezza: possiamo modificare un po' il metodo

Metodo di Eulero Inverso

Idea: definiamo il rapporto incrementale sul prossimo istante

$$\frac{x(t) - x(t-\delta)}{\delta} \simeq \dot{x}(t) = f(x(t), t)$$

Con qualche piccola manipolazione otteniamo:

$$x(t) = x(t-\delta) + \delta\, f(x(t), t)$$

Traslando di $\delta$ unità di tempo, finalmente abbiamo:

$$x(t+\delta) = x(t) + \delta\, f(x(t+\delta), t+\delta)$$

• Non è più una formula, ma una equazione non lineare
• Può essere risolta con il metodo di Newton-Raphson

Metodo di Eulero Inverso

Intuitivamente, quando risolviamo:

$$x(t+\delta) = x(t) + \delta\, f(x(t+\delta), t+\delta)$$

• Stiamo cercando un valore futuro di $x$ tale che...
• ...L'approssimazione sia consistente con il valore di $x$ corrente

Questa tecnica rende il metodo più robusto

• Il prezzo da pagare è una complessità maggiore...
• ...Ed un tempo di calcolo più elevato

Attenzione: il metodo non è necessariamente più accurato

• In generale, le differenze sono più evidenti se $\delta$ non è molto piccolo

Metodo di Eulero Inverso: Risultati

Nero: lsode, Blu: Eulero, Rosso: Eulero Inverso, 50 punti

Metodo del Punto Intermedio

L'idea base del metodo di Eulero è utilizzare l'approssimazione:

$$\frac{x(t+\delta) - x(t)}{\delta} \simeq \dot{x}(t)$$

• Cioè approssimare $\dot{x}(t)$ con il rapporto incrementale

In realtà, però, è più accurata l'approssimazione seguente:

$$\frac{x(t+\delta) - x(t)}{\delta} \simeq \dot{x}\left(\frac{\delta}{2}\right)$$

• In altre parole, il rapporto incrementale...
• ...è una approssimazione della pendenza nel punto medio

Metodo del Punto Intermedio

Procedendo come al solito otteniamo:

$$x(t+\delta) \simeq x(t) + \delta\, f\left(t + \frac{\delta}{2}, x\left(t+\frac{\delta}{2}\right)\right)$$

• La formula richiede $x(t+^\delta/_2)$, che ancora non conosciamo
• Una soluzione: lo approssimiamo con il metodo di Eulero

$$x\left(t + \frac{\delta}{2}\right) \simeq x(t) + \frac{\delta}{2} f(t, x(t))$$

Combinando, otteniamo il metodo del punto intermedio:

$$x(t+\delta) \simeq x(t) + \delta\, f\left(t + \frac{\delta}{2}, x(t) + \frac{\delta}{2} f(t, x(t))\right)$$

Metodo del Punto Intermedio

La formula del metodo del punto intermedio:

$$x(t+\delta) \simeq x(t) + \delta\, f\left(t + \frac{\delta}{2}, x(t) + \frac{\delta}{2} f(t, x(t))\right)$$

...Non è equivalente a raddoppiare i passi nel metodo di Eulero:

$$x(t+\delta) \simeq \underbrace{x(t) + \frac{\delta}{2} f(t, x(t))}_{\text{Eulero per } t + ^\delta/_2} + \frac{\delta}{2}\, f\left(t + \frac{\delta}{2}, x(t) + \frac{\delta}{2} f(t, x(t))\right)$$

• L'istante di tempo $^\delta/_2$ viene usato solo per correggere la derivata
• In questo modo si evita che gli errori si accumulino
• Aumentando i punti intermedi, si ottengono i metodi di Runge-Kutta

Metodo di Eulero Inverso: Risultati

Nero: lsode, Blu: Eulero, Rosso: Punto Intermedio, 50 punti

Elementi di Informatica e
Applicazioni Numeriche T

Un Esercizio

Un Problema di Fluidodinamica

Un serbatoio atmosferico scarica acqua attraverso una condotta

• Il serbatoio è cilindrico, di raggio 40cm
• Il livello dell'acqua nel serbatoio è di 4 metri
• La condotta è lunga 15 m
• La condotta ha diametro 10cm e scabrezza 1mm
• Per le perdite di carico dovute all'imbocco, si consideri $k = 0.5$

Quanto tempo impiega il serbatoio a svuotarsi?

Per risolvere il problema occorre utilizzare una EDO

Un Problema di Fluidodinamica

I dati del problema in sintesi:

Soluzione

• Per i dettagli della soluzione si veda il codice sul sito del corso...
• ...E la registrazione delle lezione

Prima osservazione chiave:

• Se $H$ è la quantità d'acqua nel serbatoio, allora $-Q = \dot{H}$
• La portata è la derivata della quantità d'acqua!

Seconda osservazione:

• I metodi numerici per le EDO richiedono solo di valutare $f$...
• ...Che serve a calcolare il valore della derivata...
• ...Quindi possiamo usare fsolve per calcolare $Q$