Sistemi non Lineari (Approfondimento)

Consideriamo il sistema non lineare:

$$\begin{align} & 2 x^2 = e^y \\ & \frac{1}{2} x + \frac{1}{3} y + z^2 = 0 \\ & z = \sin(\pi x) + \cos(\pi y) \end{align}$$

Come possiamo risolverlo?

Sistemi non Lineari (Approfondimento)

Consideriamo il sistema non lineare:

$$\begin{align} & 2 x^2 = e^y \\ & \frac{1}{2} x + \frac{1}{3} y + z^2 = 0 \\ & z = \sin(\pi x) + \cos(\pi y) \end{align}$$

• Prima portiamo tutte le equazioni nella forma $f_i(x, y, z) = 0$

$$\begin{align} f_1(x, y, z) &= 2 x^2 - e^y = 0 \\ f_2(x, y, z) &= \frac{1}{2} x + \frac{1}{3} y + z^2 = 0 \\ f_3(x, y, z) &= z - \sin(\pi x) - \cos(\pi y) = 0 \end{align}$$

Sistemi non Lineari (Approfondimento)

Consideriamo il sistema non lineare:

$$\begin{align} & 2 x^2 = e^y \\ & \frac{1}{2} x + \frac{1}{3} y + z^2 = 0 \\ & z = \sin(\pi x) + \cos(\pi y) \end{align}$$

• Poi vediamo il sistema come una equazione vettoriale:

$$F((x, y, z)) = \left(\begin{array}{c} f_1(x, y, z)\\ f_2(x, y, z)\\ f_3(x, y, z) \end{array}\right) = (0, 0, 0)$$

• E la risolviamo con il metodo di Newton!

Sistemi non Lineari (Approfondimento)

Una possibile soluzione (disponibile sul sito del corso)

clear all

function y = f(x)
    y(1) = 2.*x(1).^2 - e.^x(2);
    y(2) = 1./2 .* x(1) + 1./3 .* x(2) + x(3).^2;
    y(3) = sin(pi .* x(1)) + cos(pi .* x(2)) - x(3);
end

x0 = [0, 0, 0]
[sol, fval, status] = fsolve(@f, x0)

• Per dettagli si veda la registrazione della lezione

Attenzione alla convergenza!

Consideriamo una leggera variziona del sistema:

$$\begin{align} & 2 x^2 = e^y \\ & \frac{1}{2} x + \frac{1}{3} y + z = 0 \quad \text{era: } z^2\\ & z = \sin(\pi x) + \cos(\pi y) \end{align}$$

• Con la variazione indicata, il sistema è impossibile

In questo caso il metodo di Newton non converge

• In senso stretto, $f(x^{(k)})$ non si annulla mai...
• ...Però le implementazioni al PC usano delle tolleranze

Attenzione alla convergenza!

Per via delle tolleranze, possono succedere due cose:

• I valori di $x^{(k+1)}$ e $x^{(k)}$ sono vicini, anche se $f(x^{(k+1)}) \neq 0$
• I valori di $f(x^{(k+1)})$ e $f(x^{(k)})$ sono vicini, anche se $f(x^{(k+1)}) \neq 0$

In entrambi i casi la funzione fsolve di Octave termina:

• Il valore del parametro di uscita INFO è rispettivamente 2 o 3
• Ma il problema non è risolto, infatti FVAL non è nullo

Come gestire il problema?

• Controllate sempre il valore di FVAL!
• Se non è sufficientemente piccolo, non fidatevi del risultato

Elementi di Informatica e
Applicazioni Numeriche T

Minimini Quadrati
e Leggi di Potenze

Determinazione del Tipo di un Fluido

In fluidodinamica si può distinguere la classe di un fluido

La distinzione viene fatta in base alla relazione tra:

• La velocità di deformazione $\dot{\gamma}$
• E lo sforzo tangenziale $\tau$

In particolare, consideriamo tre casi:

• Fluido Newtoniano: relazione lineare $\tau = \alpha \dot{\gamma}$
• Fluido di Bingham: relazione affine $\tau = \alpha \dot{\gamma} + \beta$
• Fluido pseudoplastico: legge di potenze: $\tau = \alpha \dot{\gamma}^{\beta}$

Determinazione del Tipo di un Fluido

Su un piano cartesiano, i tre tipi di fluido si presentano così:

Determinazione del Tipo di un Fluido

La determinazione viene effettuata in modo empirico:

• Si effettua una serie di misurazioni di $\tau, \dot{\gamma}$
• Quindi si fa una ipotesi sul tipo di fluido e si stimano i parametri

In pratica, possiamo utilizzare il metodo dei minimi quadrati!

• Possiamo eseguire il metodo per ognuno dei tre casi:

$$\begin{align} (a)\quad \tau &= \alpha \dot{\gamma} & (b)\quad \tau &= \alpha \dot{\gamma} + \beta & (c)\quad \tau &= \alpha \dot{\gamma}^\beta \end{align}$$

• Scegliamo il caso che minimizza la somma dei quadrati dei residui:

$$e = \sum_{i = 1}^{n} (\tau_i - f_{\alpha,\beta}(\dot{\gamma}_i))^2$$

Un Esempio Pratico

Siano date le misurazioni di $\tau, \dot{\gamma}$ della figura seguente:

Un Esempio Pratico

Si determini il tipo di fluido ed i parametri della relazione tra $\tau$ e $\dot{\gamma}$

Un Esempio Pratico

Partiamo dal caso più semplice, ossia il fluido di Bingham

$$f_{\alpha,\beta}(\dot{\gamma}) = \alpha \dot{\gamma} + \beta$$

• Ci basta fare una interpolazione con un polinomio di grado 1!

x = % misurazioni di gamma
y = % misurazioni di tau

% Minimi quadrati
p1 = polyfit(x, y, 1) % p1 = [1.07631, 0.82395]

% Funzione approssimante
f1 = @(x) (polyval(p1, x))

Un Esempio Pratico

Vediamo su un grafico la qualità dell'approssimazione:

Un Esempio Pratico

In termini numerici, guardiamo alla somma dei quadrati dei residui:

$$e = \sum_{i = 1}^{n} (\tau_i - f_{\alpha,\beta}(\dot{\gamma}_i))^2$$

In Octave:

E1 = sum((y - f1(x)).^2)

• Il risultato in questo caso è $1.2065$
• Una soluzione completa è disponibile sul sito del corso
• Per dettagli si veda la registrazione della lezione

Un Esempio Pratico

Proviamo ipotizzare che il fluido sia Newtoniano, quindi:

$$f_{\alpha,\beta}(\dot{\gamma}) = \alpha \dot{\gamma}$$

• Non possiamo usare polyfit: non garantirebbe di avere $\beta = 0$

Però possiamo eseguire il metodo risolvendo un sistema lineare!

• Dalla teoria classica abbiamo che i parametri $\beta$ sono dati da:

$$\beta = (\Phi^T\Phi) \Phi^T y$$

• Dove la funzione approssimante è nella forma:

$$f_\beta(x) = \beta_1 g_1(x) + \beta_2 g_2(x) + \ldots$$

Un Esempio Pratico

Nel nostro caso, la funzione approssimante ha un solo componente:

$$g_1(x) \leftrightarrow \dot{\gamma}$$

• Ed abbiamo un solo parametro, i.e. il coefficiente $\alpha$

Quindi la matrice $\Phi$ sara fatta come segue:

$$\Phi = \left(\begin{array}{c} g_1(x_1)\\ g_2(x_2)\\ \vdots \end{array}\right) = \left(\begin{array}{c} \dot{\gamma}_1\\ \dot{\gamma}_2\\ \vdots \end{array}\right)$$

• $\Phi$ è semplicemente il vettore delle nostre misurazioni $\dot{\gamma}$!

Un Esempio Pratico

In Octave, assumendo di avere un fluido Newtoniano, otteniamo:

x = % colonna con le misurazioni di gamma
y = % colonna con la misurazioni di tau

Phi = x
p2 = Phi \ y % p = alpha = 1.6901

f2 = @(x) (p2 .* x) % Funzione approssimante

E2 = sum((y - f2(x)).^2)

• E2 vale $15.990$
• Molto peggio di prima: il fluido decisamente non è Newtoniano

Un Esempio Pratico

Rimane da considerare solo il caso del fluido pseudoplastico, con:

$$f_{\alpha,\beta}(\dot{\gamma}) = \alpha \dot{\gamma}^\beta$$

Di nuovo non possiamo usare polyfit, ma per una ragione diversa:

• Il problema è che la funzione è non lineare in $\beta$...
• ...Ma noi abbiamo visto solo il metodo lineare dei minimi quadrati!

Il metodo non si applica se $f$ dipende non linearmente dai parametri

In realtà, c'è una soluzione semplice: possiamo "barare" ;-)

• Possiamo cambiare la scala dei dati...
• ...in modo che la funzione diventi lineare!

Un Esempio Pratico

Guardiamo cosa succede considerando $\log(f_{\alpha,\beta})$ anziché $f_{\alpha,\beta}$:

$$\log\, f_{\alpha,\beta}(\dot{\gamma}) = \log \left( \alpha \dot{\gamma}^\beta\right) = \log \alpha + \beta \log x$$

In altre parole:

• Il logaritmo di $f_{\alpha,\beta}(\dot{\gamma})$...
• ...dipende linearmente dal logaritmo di $x$

Quindi, utilizzare il metodo dei minimi quadrati, trasformiamo i dati:

$$\begin{align} & \tau_i \longrightarrow \log \tau_i \\ & \dot{\gamma}_i \longrightarrow \log \dot{\gamma}_i \\ \end{align}$$

Un Esempio Pratico

Questi sono i dati trasformati su un grafico cartesiano:

Un Esempio Pratico

Sembra proprio una retta => è probabile che il fluido si pseudoplastico

Un Esempio Pratico

A questo punto possiamo approssimare i dati trasformati:

x = % misurazioni di gamma
y = % misurazioni di tau

xl = log(x) % trasformo gamma
yl = log(y) % trasformo tau

p3 = polyfit(xl, yl) % risultato: [0.49591, 0.68855]

Poiché abbiamo $\log\, f_{\alpha,\beta}(\dot{\gamma}) = \log \alpha + \beta \log x$, allora:

• p3(1) contiene $\beta$
• p3(2) contiene $\log \alpha$

Un Esempio Pratico

A questo punto possiamo applicare ai parametri alla forma iniziale:

$$f_{\alpha,\beta}(\dot{\gamma}) = \alpha\, x^\beta$$

In Octave, abbiamo:

% p3(1) = alpha = 0.49591
% p3(2) = log(beta) = 0.68855
% Quindi la funzione approssimante è:
f3 = @(x) (e.^p3(2) .* x.^p3(1) )

E3 = sum((y - f3(x)).^2)

• Il valore di E3 è $0.58175$, il migliore finora
• Il fluido è pseudoplastico con $\alpha \simeq 0.49591, \beta \simeq e^{0.68855}$

Un Esempio Pratico

Vediamo la qualità dell'approssimazione su un grafico:

Elementi di Informatica e
Applicazioni Numeriche T

Lookup di Funzioni
Definite per Punti

Determinare lo Stato ad un Tempo Dato

Consideriamo il nostro esempio delle BMW

Lo stato è descritto da $(x, v)$ e si evolve secondo l'EDO:

$$\dot{v} = \frac{1}{m} \left( F - \frac{1}{2} \rho v^2 C_D A \right)$$

Dove abbiamo che:

• $F = 1000\, N$, $m = 161\, Kg$
• $\rho \simeq 1.25\, Kg/m^3$ (densità dell'aria)
• $A = 1.2\, m^2$ (superficie esposta)
• $C_D = 0.26$ (resistenza aerodinamica)

Domanda: quanto vale la velocità dopo $40$ secondi?

Determinare lo Stato ad un Tempo Dato

Il modo più semplice per rispondere è risolvere l'EDO fino a $t = 40$:

function dx = f(x, t)
    ... % Definizione delle variabili
    % x(1) = posizione, x(2) = velocità
    dx(1) = x(2);
    dx(2) = 1./m .* (F - 0.5.*rho.* x(2).^2.*C_D.*A);
end
t = linspace(0, 40);
x0 = [0, 0];
sol = lsode(@f, x0, t)
v = sol(:, 2)
v40 = v(end) % Il risultato!

• Alla fine il risultato è l'ultimo elemento nel vettore v

Determinare lo Stato ad un Tempo Dato

Proviamo a complicare la cosa:

• Determiniamo la velocità dopo 40 secondi...
• ...E la strada percorsa dopo 20 secondi

Soluzione banale:

• Risolvere la EDO fino a $t = 40$ per ottenere la velocità
• Risolvere la EDO (di nuovo) fino a $t = 20$ per ottenere la strada

Così facendo però ripetiamo i conti due volte

Vediamo un metodo alternativo

Determinare lo Stato ad un Tempo Dato

La soluzione della EDO è un vettore di valutazioni della funzione $v(t)$:

Determinare lo Stato ad un Tempo Dato

Se ci interessa il valore di $v(t_0)$ per un $t_0$ che è stato valutato:

• Ci basta recuperare l'indice i0 di $t_0$ nel vettore t...
• Il risultato poi è dato da v(i0)

Determinare lo Stato ad un Tempo Dato

Se ci interessa $v(t_0)$ per un $t_0$ che è non stato valutato:

• Cerchiamo un indice i per cui t(i) >= t0 && t(i+1) <= t0
• Oppure: t(i) <= t0 && t(i+1) >= t0 (se t non è un "tempo")

Determinare lo Stato ad un Tempo Dato

Se ci interessa $v(t_0)$ per un $t_0$ che è non stato valutato:

• Poi per esempio otteniamo v(i)

Determinare lo Stato ad un Tempo Dato

Se ci interessa $v(t_0)$ per un $t_0$ che è non stato valutato:

• Poi per esempio otteniamo v(i)
• Oppure possiamo usare usare una approssimazione lineare

Determinare lo Stato ad un Tempo Dato

In sintesi, dati due vettori v e t:

• Troviamo un indice i tale che t "attraversa" t0
• Poi otteniamo il valore di v (eventualmente interpolando)

Determinare lo Stato ad un Tempo Dato

Equivale ad approssimare $y = f(x)$ con una spezzata!

• Punto sinistro = approssimazione a gradini

Determinare lo Stato ad un Tempo Dato

Equivale ad approssimare $y = f(x)$ con una spezzata!

• Punto sinistro = approssimazione a gradini
• Interpolazione lineare = approssimazione lineare a tratti

Determinare lo Stato ad un Tempo Dato

Il codice in Octave:

function xx = lookup_interp(t, x, t0)
    for ii = 1:length(t)-1
        if (t(ii) <= t0 && t(ii+1) >= t0) || ...
           (t(ii+1) <= t0 && t(ii) >= t0)
            % Coefficiente angolare
            m = (x(ii+1)-x(ii)) ./ (t(ii+1)-t(ii));
            % Interpolazione lineare
            xx = x(ii) + m .* (t0 - t(ii));
            break
        end
    end
end

Utilizzi dell l'Interpolazione con Spezzate

Quanto tempo di vuole per arrivare a $100\, Km/h$?

• Ho una condizione su $v$ (anziché su $t$)
• Mi interessa il valore di $t$ (anziché di $v$)

lookup_interp(v, t, 27.78) % 100 Km/h = 27.78 m/s

Quanto vale la velocità dopo $100\, m$?

• Ho una condizione su $x$
• Mi interessa il valore di $v$

lookup_interp(x, v, 100)

Interpolazione con Spezzate: interp1

In alternativa, potete usare la funzione predefinita:

function Y1 = interp1(X, Y, XI)

...Che restituisce il valore interpolato di Y quando X vale XI

• Di default, interp1 utilizza una interpolazione lineare
• Si può modificarne il comportamento (vedere help)

Attenzione: interp1 richiede che X sia ordinato

• Se X rappresenta il tempo, allora l'ordinamento è naturale
• Se X rappresenta qualcos'altro (e.g. strada) allora non è detto

Buone notizie: la nostra lookup_interp non ha questo problema!

Un Altro Esempio

Consideriamo ancora le nostre montagne russe:

• Dove: $a = 15\, m, b = 20\, m$
• Domanda: quanto vale la velocità nel punto $a$?

Un Altro Esempio

Vediamo le parti importanti della soluzione in Octave:

function y = f(x, t)
    % x(1) = pos., x(2) = vel. (vars: da definire)
    fp = (2.*b)./a.^2 .* x(1) + - (2.*b)./a;
    y(1) = x(2) .* (1 ./ sqrt(1 + fp.^2));
    y(2) = -g .* (fp ./ sqrt(1 + fp.^2));
end

t = linspace(0, 4); % In 4 secondi oltrepasso "a"
x0 = [0, 0];
sol = lsode(@f, x0, t);
v = sol(:, 2)
VA = lookup_interp(x, v, a); % Risultato: 1.6633

Elementi di Informatica e
Applicazioni Numeriche T

EDO e Stato Stazionario

Determinazione dello Stato Stazionario

Torniamo all'esempio delle BMW

Lo stato è descritto da $(x, v)$ e si evolve secondo l'EDO:

$$\dot{v} = \frac{1}{m} \left( F - \frac{1}{2} \rho v^2 C_D A \right)$$

Abbiamo sempre che:

• $F = 1000\, N$, $m = 161\, Kg$
• $\rho \simeq 1.25\, Kg/m^3$ (densità dell'aria)
• $A = 1.2\, m^2$ (superficie esposta)
• $C_D = 0.26$ (resistenza aerodinamica)

Nuova Domanda: quanto vale la velocità massima?

Determinazione dello Stato Stazionario

In questo caso la velocità massima è la velocità "a regime"

Determinazione dello Stato Stazionario

Lo stato stazionario di $v$ viene raggiunto quando $\dot{v} = 0$

• Risolviamo quindi la EDO fino ad un tempo sufficientemente grande:

t = linspace(0, 60);
sol = lsode(@f, [0, 0], t);
v = sol(:, 2)

• Recuperiamo i valori di $\dot{v} = \dot{v} = \frac{1}{m} \left( F - \frac{1}{2} \rho v^2 C_D A \right)$

dv = 1./m .* (F - 0.5 .* rho .* v.^2 .* C_D .* A)

• Infine recuperiamo il valore di $v$ quando $\dot{v} = 0$

% Uso una tolleranza: $\dot{v} = 0$ a +$\infty$
vmax = lookup_interp(dv, v, 1e-5)

Un Metodo Alternativo

Consideriamo la formula fondamentale del metodo di Eulero:

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

È la descrizione di un sistema dinamico tempo-discreto

• Quidi possiamo cercare lo stato stazionario con IPF
• In pratica, si tratta di eseguire il metodo di Eulero...
• ...Finché $x(t+\delta)$ e $x(t)$ non sono sufficientemente vicini

Rispetto al metodo precedente:

• Più efficiente: ci fermiamo allo stato stazionario!
• Meno accurato: lsode usa delle approssimazioni molto migliori

Elementi di Informatica e
Applicazioni Numeriche T

Potenzialità del Problema
ai Valori Iniziali

Problema ai Valori Iniziali?

Consideriamo ancora l'esempio delle montagne russe:

Domanda: quanto deve valere $v$ nel punto $0$ perché in $a$ sia $20\, m/s$?

Problema ai Valori Iniziali?

Proviamo a modellare 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}$$

• Dove $\phi'(x) = 2\, {^b/_{a^2}} x - 2\, {^b/_a}$

Sappiamo che, per un qualche tempo $f_f$ deve valere:

$$\begin{align} & x(t_f) = a && v(t_f) = 20 \end{align}$$

Possiamo risolverlo con i metodi visti finora?

Problema ai Valori Iniziali!

Provate a pensare dal metodo di Eulero:

• Se potessimo eseguirlo "a ritroso", partendo da $x = a$ e $v = 20$...
• ...Dovremmo poter determinare la velocità quando $x$ vale $0$

Intuitivamente:

• Se potessimo "invertire" la direzione del tempo...
• ...Ci ritroveremmo con un nuovo problema ai valori iniziali

$$\begin{align} & \dot{x}'(\tau) =\ ? \\ & \dot{v}'(\tau) =\ ? \\ & x'(0) = a & v'(0) = 20\\ \end{align}$$

Dobbiamo solo determinare le funzioni per il calcolo di $\dot{x}'(\tau), \dot{v}'(\tau)$

Inversione Temporale di una EDO

Per come l'abbiamo definita, la nuova variabile $\tau$ corrisponde a:

$$\tau = g(t) = t_f - t$$

Quindi, se ci riferiamo con $x'(\tau)$ a $x(g(\tau))$, possiamo ottenere:

$$\dot{x}'(\tau) = \dot{x}(g(t))\, g'(t) = - f(x(g(t), g(t)) = - f(x', t_f - t)$$

• Dove $f(x, t)$ è la funzione che definiva $\dot{x}(t)$

Quindi possiamo ottenere la nuova EDO con:

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

• Ossia negando la funzione $f(x, t)$ originaria
• E sostituendo ogni occorrenza di $t$ in $f(x, t)$ con $t_f - t$

Un Esempio Pratico

Nel nostro caso, il problema con il tempo invertito è dato da:

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

Per semplicità le variabili non sono state rinominate

• Se la variabile $t$ fosse comparsa nella espressione di $\dot{x}, \dot{v}$...
• ...Avremmo dovuto sostituirla con $t_f - t$

In questo caso, avremmo dovuto conoscere il valore di $t_f$

Un Esempio Pratico

Vediamo la soluzione (disegnata con plotyy, se vi interessa)

Un Esempio Pratico

Per ottenere il valore di velocità richiesto:

xf = [a, 20]; % Condizioni al contorno
tau = linspace(0, 4); % Tempo invertito
sol = lsode(@f, xf, tau);

% Ottengo i due vettori che compongono la soluzione
x = sol(:, 1);
v = sol(:, 2);

v0 = lookup_interp(x, v, 0) % velocità per x = 0

• Il risultato è $2.7558\, m/s$
• Il codice è disponibile sul sito del corso
• Per dettagli si veda la registrazione della funzione