Funzione fprintf

La funzione fprintf stampa con formattazione controllata:

L’interfaccia della funzione è:

function fprintf(FORMAT, E1, E2, ...)

• FORMAT è una stringa da stampare
• E1, E2 sono espressioni (e.g. nomi di variabile)
• FORMAT può contenere dei “segnaposto”, con il carattere %
  • E.g. %f è un segnaposto per un valore reale
  • E.g. %d è un segnaposto per un valore intero
• Il primo segnaposto è sostituito con il valore di E1
• Il secondo segnaposto è sostituito con il valore di E2, etc.

Funzione fprintf

Vediamo un esempio:

A = 10.5;
B = 2;
fprintf('A = %f, B = %f, A*B = %f\n', A, B, A*B)

Stampa: A = 10.500000, B = 2.000000, A*B = 21.000000

• “\n” è un carattere speciale e serve ad andare a capo

I segnaposto cono configurabili

Ci interessa un caso solo: “%.Nf” stampa un reale con N valori decimali

• E.g. %.3f stampa con 3 cifre decimali
• E.g. %.1f stampa con 1 cifra decimale

Ciclo while

Il ciclo for non è sempre adeguato ad implementare iterazioni:

• E.g. quando non è desiderabile limitare il numero di iterazioni a priori

Per questi casi, Matlab fornisce il ciclo while

La sintassi è:

while <espressione> % vera o falsa
  <corpo>
end

• Il corpo viene ripetuto…
• …Fintanto che <espressione> denota true (o \(\neq 0\))

Funzione Zeta di Riemann

Consideriamo la funzione Zeta di Riemann ed il nostro vecchio codice:

\[\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}\]

z = 0;                      % val. della somma
old_z = -Inf;               % vecchio z
for n = 1:10000               % 1e5 = iterazioni massime
  z = z + 1 ./ n.^s;
  if abs(z - old_z) < 1e-6  % 1e-6 è la tolleranza
    break                   % Interrompe il ciclo
  end
  old_z = z;                % rimpiazzo il vecchio z
end

• Cosa succede se 10000 iterazioni non bastano a convergere?

Funzione Zeta di Riemann

Consideriamo la funzione Zeta di Riemann ed il nostro vecchio codice:

\[\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}\]

z = 0;                      % val. della somma
n = 1;
old_z = -Inf;               % vecchio z
while abs(z - old_z) > 1e-6 % 1e-6 è la tolleranza
  old_z = z;                % memorizzo il vecchio z
  z = z + 1 ./ n.^s;
  n = n + 1;                % incremento n
end

• Fintanto che \(|z - old_z| > 1e-6\), il ciclo prosegue

Funzione Zeta di Riemann

Consideriamo la funzione Zeta di Riemann ed il nostro vecchio codice:

\[\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}\]

z = 0;                      % val. della somma
n = 1;
old_z = -Inf;               % vecchio z
while abs(z - old_z) > 1e-6 % 1e-6 è la tolleranza
  old_z = z;                % memorizzo il vecchio z
  z = z + 1 ./ n.^s;
  n = n + 1;                % incremento n
end

• Dobbiamo però gestire l’indice n esplicitamente

Cicli Infiniti

Il ciclo while non fornisce garanzie di terminazione

Per esempio:

n = 1;
s = 0;
while n < 10
  s = s + n;
end

• Questo ciclo non termina
• Perché n non viene incrementato!

Se vi capita, niente panico: basta premere [CTRL+C]

• Capita più spesso di quanto ci si possa aspettare :-)

Equilibrio di Sistemi Tempo-Discreti

Consideriamo un sistema dinamico tempo-discreto

In generale è definito da una equazione del tipo:

\[x^{(t+1)} = f(x^{(t)})\]

Per questo tipo di sistemi, abbiamo imparato a:

• Osservare l’andamento dello stato del tempo
• Identificare il tipo di comportamento (convergente, periodico…)
• Per i sistemi convergenti, individuare uno stato di equilibrio…
• …Simulando sufficientemente a lungo e misurando lo stato finale

Gli stati di equilibrio, però, si possono determinare a priori!

Determinazione degli Stati di Equilibrio

Uno stato è di equilibrio se viene “trasformato in se stesso”

Formalmente, uno stato di equilibrio \(x\) deve soddisfare:

\[x = f(x)\]

• Manca l’indice di tempo, perché lo stato a sx e dx è lo stesso

Se risolviamo l’equazione, determiniamo gli stati di equilibrio

• In generale \(x\) è un vettore, quindi \(x = f(x)\) è un sistema di equazioni
• Se \(f(x)\) è lineare, possiamo usare la forma matriciale:

\[A x = b\]

E possiamo risolverlo con i metodi visti in analisi numerica

Un Esempio: Pagerank

Per l’algoritmo pagerank, visto la scorsa lezione

L’equazione fondamentale è data da:

\[x^{(t+1)} = p \underbrace{\frac{1}{n}}_{\text{vettore}} + (1-p) P x^{(t)}\]

Uguagliando \(x^{(t+1)}\) e \(x^{(t)}\) otteniamo:

\[I x = p \frac{1}{n} - (1-p) P x\]

\[\underbrace{(I - (1-p) P)}_{A} x = \underbrace{p \frac{1}{n}}_{b}\]

• \(A\) è quadrata per costruzione (deriva da una transizione di stato)
• Quindi, in casi normali, la soluzione è data da \(x = A^{-1} b\)

Un Esempio: Pagerank

Supponendo di disporre delle variabili:

• p, per la probabilità di stancarsi
• P, per la matrice delle probabilità di click
• n, per il numero delle pagine

Possiamo prima costruire la matrice \(A\) e la colonna \(b\):

A = (eye(n) - (1-p)*P) % eye e' la matrice identita'
b = (ones(n,1) * p/n) % ones(n,1) per avere una colonna

E quindi possiamo calcolare la soluzione con una divisione sinistra:

xeq = A \ b

Equilibrio di Sistemi Tempo-Discreti Non Lineari

L’approccio è valido anche per sistemi dinamici non lineari

Uno stato, per essere di equilibrio, deve soddisfare:

\[x = f(x)\]

• I.e. la relazione che abbiamo già visto!
• …Che corrisponderà però \(x\) ad un sistema di equazioni non lineari

Ce ne occuperemo più avanti nel corso

Simulare o Non Simulare?

Quando risolvere un sistema di equazioni e quando simulare?

Simuliamo se:

• Ci interessa il comportamento nel tempo
• Il sistema non ha uno stato di equilibrio (e.g. periodico, caotico)
• Vogliamo determinare (empiricamente) la stabilità dell’equilibrio
• …

Risolviamo le equazioni se:

• Non ci interessa il comportamento nel tempo
• Non ci interessa la stabilità (basta uno stato di equilibrio)
• Se dobbiamo calcolare lo stato di equilibrio con alta precisione
• …

Stati di Equilibrio Multipli

Un sistema dinamico tempo discreto lineare:

• Ha sempre un solo stato di equilibrio…
• …A meno che la matrice dei coefficienti non sia singolare

Per esempio, per le previsioni del tempo (scorsa lezione) avevamo:

\[\begin{align} x^{(t+1)} = P x^{(t)} \quad\text{ con }\quad P = \left(\begin{array}{cc} 0.9 & 0.5 \\ 0.1 & 0.5 \\ \end{array}\right) \end{align}\]

Da cui si ottiene il sistema:

\[(I - P) x = 0 \quad\text{ con }\quad \begin{aligned} (I-P) = \left(\begin{array}{cc} 0.1 & -0.5 \\ -0.1 & 0.5 \\ \end{array}\right) \end{aligned}\]

• Possiamo usare la funzione det per calcolare il determinante:

det(eye(2) - A) % -1.1102e-17 (sarebbe 0)

Stati di Equilibrio Multipli

Cosa vuol dire in pratica?

• Un sistema sotto-determinato ha infinite soluzioni…
• …Quindi ci sono infiniti stati di equilibrio

Per esempio, supponiamo di avere due vasi comunicanti

• Il livello finale dipende da quanta acqua c’è nel sistema!
• In questo specifico caso, lo stato raggiunto dipende dallo stato iniziale…
• …Altre volte, ci sono dei vincoli sottointesi

Esercizio: Previsioni del Tempo

Nel caso delle previsioni del tempo:

Il tempo è bello o brutto, quindi la somma delle probabilità deve essere 1

• Partiamo dalle equazioni originali per l’equilibrio:

\[\begin{align} 0.1 x_g - 0.5 x_b = 0\\ -0.1 x_g + 0.5 x_b = 0 \end{align}\]

• Una delle due equazioni può essere rimossa…
• …Ma possiamo anche aggiungere \(x_g + x_b = 1\)

\[\begin{align} 0.1 x_g - 0.5 x_b = 0\\ x_g + x_b = 1 \end{align}\]

Nel file es_weather.m:

• Determinate lo stato di equilibrio risolvendo il sistema lineare
• Verificate che coincide con lo stato raggiunto dalla simulazione

Equilibrio di Sistemi Tempo-Continui Lineari

Una stanza è ventilata mediante una sola apertura:

• Sia \(T_o\) temperatura esterna, \(T_a\) quella dell’aria interna
• Sia \(T_w\) la temperatura dei muri

Il flusso di calore tra l’esterno e l’aria è dato da:

\[i_1 = \frac{1}{R_1}(T_o - T_a)\]

• Dove \(R_1\) è la resistenza termica dell’apertura

Il flusso di calore tra l’aria e i muri è dato da:

\[i_2 = \frac{1}{R_2}(T_a - T_w)\]

• Dove \(R_1\) è la resistenza termica tra l’aria ed i muri

Equilibrio di Sistemi Tempo-Continui Lineari

Per quanto riguarda le temperature:

• La temperatura dell’esterno e dei muri si suppone costante
• La temperatura dell’aria varia con i flussi di calore

In particolare vale la relazione:

\[\frac{d T_a}{d t} = \frac{1}{C_a} (i_1 - i_2)\]

• Il flusso di calore \(i_1\) va dall’esterno all’aria, mentre \(i_2\) va dall’aria ai muri
• \(C_a\) è la capacità termica dell’aria

Per completezza, ricordiamo che:

\[\begin{align} i_1 &= \frac{1}{R_1}(T_o - T_a) & i_2 &= \frac{1}{R_2}(T_a - T_w) \end{align}\]

Equilibrio di Sistemi Tempo-Continui Lineari

Questo è un primo esempio di sistema dinamico tempo-continuo

I sistemi dinamici tempo continui sono descritti mediante equazioni differenziali

• Tipicamente, queste sono nella forma:

\[\dot{x} = f(x)\]

• Non sappiamo risolvere le equazioni differenziali…
• …Vedremo come farlo verso la fine del corso

Possiamo però già osservare che:

• Per definizione all’equilibrio lo stato non varia…
• …Il che vuol dire che le derivate si annullano

Il generale, un sistema tempo continuo all’equilibrio deve soddisfare:

\[f(x) = 0\]

Equilibrio di Sistemi Tempo-Continui Lineari

Nel nostro esempio abbiamo:

\[\begin{align} \frac{d T_a}{d t} &= \frac{1}{C_a} (i_1 - i_2) & i_1 &= \frac{1}{R_1}(T_o - T_a) & i_2 = \frac{1}{R_2}(T_a - T_w) \end{align}\]

Quindi, all’equilibrio avremo:

\[\begin{align} 0 &= \frac{1}{C_a} (i_1 - i_2) \\ i_1 &= \frac{1}{R_1}(T_o - T_a) \\ i_2 &= \frac{1}{R_2}(T_a - T_w) \end{align}\]

• È un sistema di equazioni lineari in \(i_1, i_2, T_a\)

E questo sappiamo come risolverlo!

Esercizio: Temperatura di una Stanza

Partiamo dal sistema originale:

\[\begin{align} 0 &= \frac{1}{C_a} (i_1 - i_2) \\ i_1 &= \frac{1}{R_1}(T_o - T_a) \\ i_2 &= \frac{1}{R_2}(T_a - T_w) \end{align}\]

• Le variabili sono \(i_1, i_2, T_a\)…
• …Ci conviene portarle a sx del segno =

Esercizio: Temperatura di una Stanza

Partiamo dal sistema originale:

\[\begin{align} \frac{1}{C_a} (i_1 - i_2) &=0 \\ i_1 + \frac{1}{R_1} T_a &= \frac{1}{R_1}T_o \\ i_2 - \frac{1}{R_2}T_a &= \frac{1}{R_2} T_w \end{align}\]

• Ora possiamo portarlo in forma matriciale

Esercizio: Temperatura di una Stanza

Partiamo dal sistema originale:

\[\left(\begin{array}{ccc} \frac{1}{C_a} & - \frac{1}{C_a} & \\ 1 & & \frac{1}{R_1} \\ & 1 & -\frac{1}{R_2} \\ \end{array}\right) \left(\begin{array}{c} T_a \\ i_1 \\ i_2 \end{array}\right) = \left(\begin{array}{c} 0 \\ \frac{1}{R_1}T_o \\ -\frac{1}{R_2}T_w \end{array}\right)\]

• Ogni riga è la trascrizione di una equazione
• Ogni colonna della matrice è associata ad una variabile…
• …Perché viene moltiplicata per tale variabile

Se chiamiamo la matrice \(A\) ed il termine noto \(b\), abbiamo:

\[A x = b\]

• Possiamo risolverlo con una divisione sinistra!

Esercizio: Temperatura di una Stanza

Partite dal file es_temperature.m nello start-kit

• Impostate la matrice dei coefficienti
• Impostate la colonna dei termini noti
• Risolvete il sistema per \(T_o = 17, T_w = 21\) e per \(T_o = 28, T_w = 21\)
• Stampate il valore di \(T_a\) all’equilibrio

Problemi di Interpolazione (Vincolata)

Si vuole progettare una arcata a ridosso di una parete verticale

La curva che descrive l’arcata:

• Deve essere ancorata ad un punto noto \((x_0, y_0)\) sulla parete
• Deve essere ancorata ad un punto noto \((x_2, y_2)\) a terra
• Deve raggiungere l’altezza massima per \(x = x_1\) (con \(x_1\) noto)

Problemi di Interpolazione (Vincolata)

Un approccio: trattiamo la curva come una funzione \(f(x)\)

In questo modo possiamo tradurre le condizioni in equazioni:

• Deve essere ancorata ad un punto noto \((x_0, y_0)\) sulla parete

\[f(x_0) = y_0\]

• Deve essere ancorata ad un punto noto \((x_2, y_2)\) a terra

\[f(x_2) = y_2\]

• Deve raggiungere l’altezza massima per \(x = x_1\) (con \(x_1\) noto)

\[f'(x_1) = y_1\]

Così come sono ci dicono ben poco…

Tracciamento di Curve

Ci serve una assunzione sulla classe della funzione \(f(x)\)

Per esempio: \(f(x)\) è polinomiale. Formalmente:

\[f(x) = \sum_{i=0}^n \alpha_i x^i\]

Le nostre condizioni allora diventano:

\[\begin{align} & \text{passaggio per } (x_0, y_0) && \sum_{i=0}^n \alpha_i x_0^i = y_0 \\ & \text{passaggio per } (x_1, y_1) && \sum_{i=0}^n \alpha_i x_2^i = y_2 \\ & \text{annullamento di } f'(x_1) && \sum_{i=1}^n i \alpha_i x_1^{i-1} = 0 \end{align}\]

Tracciamento di Curve

Ci siamo quasi! Guardiamole meglio:

\[\begin{align} & \text{passaggio per } (x_0, y_0) && \sum_{i=0}^n \alpha_i x_0^i = y_0 \\ & \text{passaggio per } (x_1, y_1) && \sum_{i=0}^n \alpha_i x_2^i = y_2 \\ & \text{annullamento di } f'(x_1) && \sum_{i=1}^n i \alpha_i x_1^{i-1} = 0 \end{align}\]

Quali sono le incognite?

• Sono i parametri della funzione \(\alpha_i\) e non le \(x\)!

Che grado di polinomio ci serve?

• Tre condizioni \(\Rightarrow\) tre variabili, i.e. \(\alpha_0, \alpha_1, \alpha_2\) \(\Rightarrow\) secondo grado

Tracciamento di Curve

In questo modo otteniamo il sistema:

\[\begin{align} \alpha_2 x_0^2 + \alpha_1 x_0 + \alpha_0 &= y_0 \\ \alpha_2 x_2^2 + \alpha_1 x_2 + \alpha_0 &= y_2 \\ 2 \alpha_2 x_1 + \alpha_1 &= 0 \end{align}\]

• Che è lineare nelle incognite \(\alpha_2, \alpha_1, \alpha_0\)

La tecnica vista è un metodo generale per progettare curve:

• Si ipotizza una struttura per la curva da costruire (e.g. polinomio)
• Si traducono i vincoli del problema in equazioni
• Si risolvono le equazioni per determinare i parametri
  • Per ora consideriamo il caso in cui le equazioni sono lineari

Esercizio: Progettazione di un’Arcata

Consideriamo il sistema per il problema di progettazione dell’arcata:

\[\begin{align} \alpha_2 x_0^2 + \alpha_1 x_0 + \alpha_0 &= y_0 \\ \alpha_2 x_2^2 + \alpha_1 x_2 + \alpha_0 &= y_2 \\ 2 \alpha_2 x_1 + \alpha_1 &= 0 \end{align}\]

A partire dal file es_arc.m nello start-kit:

• Si determinino i coefficienti della curva risolvendo il sistema
• Si disegni la forma dell’arcata
• Si stampi a video il valore dell’altezza massima

Esercizio: Letto di un Fiume

Si vuole progettare lo scavo per il letto di un fiume

La sezione dello scavo deve presentarsi come segue:

• La coordinata \(x\) rappresenta una posizione orizzontale
• La coordinata \(d\) rappresenta la profondità dello scavo
  • Per questa ragione la sezione si presenta “al contrario”

Esercizio: Letto di un Fiume

Si vuole progettare lo scavo per il letto di un fiume

• La sezione deve essere descritta da una curva parabolica
• Deve passare per i punti noti \((x_0, y_0)\) e \((x_1, y_1)\)
• L’area della sezione determina la portata massima…
• …E deve essere pari ad un valore prestabilito \(s_1\)

Se \(f(x)\) è la funzione che descrive la curva, l’area della sezione è:

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

A partire dal file es_riverbed.m nello start-kit:

• Si determinino i coefficienti della curva
• Si disegni la forma della sezione

Esercizio: Controllo di Accelerazione

Si vuole controllare l’accelerazione di un carrello automatico

Il profilo di velocità in accelerazione deve presentarsi come segue:

• La coordinata \(t\) rappresenta il tempo
• La coordinata \(v\) rappresenta la velocità

Curve di questo tipo si utilizzano nelle centraline di controllo di auto e moto

Esercizio: Controllo di Accelerazione

Si vuole controllare l’accelerazione di un carrello automatico

• Il profilo deve seguire un andamento polinomiale
  • Il grado del polinomio è da determinare
  • Servirà un coefficiente per ogni condizione specificata
• La velocità iniziale e finale ed il tempo di accelerazione \(t_1\) sono noti
• Perché le variazioni non siano troppo brusca…
• …Si richiede che la derivata della velocità in \(t_0\) e \(t_1\) sia nulla

A partire dal file es_acceleration.m nello start-kit:

• Si determinino i coefficienti della curva
• Si disegni il profilo della velocità in accelerazione

Esercizio: Controllo di Frenata

Si vuole controllare l’arresto di un carrello automatico

Il profilo di velocità in frenata deve presentarsi come segue:

• La coordinata \(t\) rappresenta il tempo
• La coordinata \(v\) rappresenta la velocità

Possiamo usare la curva per programmare una centralina di controllo

Esercizio: Controllo di Frenata

Si vuole controllare l’arresto di un carrello automatico

• Il profilo è dato da un polinomio, di grado da determinare
• La velocità iniziale e finale ed il tempo di frenata \(t_1\) sono noti
• Si richiede che la derivata della velocità in \(t_1\) sia nulla
• Lo spazio di frenata \(S\) deve essere pari ad un valore \(s_1\), dove:

\[S = \int_{\underbrace{t_0}_{=0}}^{t_1} f(t) \, dt\]

A partire dal file es_brake.m nello start-kit:

• Si determinino i coefficienti della curva
• Si disegni il profilo della velocità in frenata

Esercizio: Progettazione di un Telaio

Si vuole progettare un telaio per una bicicletta

La forma del telaio deve apparire come segue:

Esercizio: Progettazione di un Telaio

Si vuole progettare un telaio per una bicicletta

• La forma del telaio è descritta da due curve paraboliche \(f_1\) ed \(f_2\)
• Le due curve originano in un punto comune \((x_0 y_0)\)
• La curva \(f_1\) deve passare per l’ancoraggio della sella in \((x_1, y_1)\)
• La curva \(f_2\) deve passare per l’ancoraggio dei pedali in \((x_2, y_2)\)
• Le due curve devono congiungersi in un punto \((x_3, y_3)\)…
• …Di cui è nota solo la coordinata \(x_3\)
• Le derivate di \(f_1\) in \(x_0\) ed \(f_2\) in \(x_3\) devono essere uguali

Esercizio: Progettazione di un Telaio

A partire dal file es_frame.m nello start-kit:

• Si determinino i coefficienti delle due curve
• Si disegni il profilo del telaio

Attenzione:

Ci sono due condizioni che coinvolgono entrambe le curve:

• Le due curve devono congiungersi in un punto \((x_3, ???)\)
• Le derivate di \(f_1\) in \(x_0\) ed \(f_2\) in \(x_3\) devono essere uguali

Non possono essere formulate separatamente!

• Occorrerà definire un’unico sistema di equazioni…
• …in cui compaiono sia \(f_1(x)\) che \(f_2(x)\)