Comportamento di un Fluido in Condotta

Comportamento di un Fluido in Condotta

Per studiare il moto di un fluido in condotta:

• Possiamo ragionare sul bilancio di potenza meccanica...
• Per ottenere una equazione che caratterizza il moto

Nella forma più semplice possibile, abbiamo che:

$$H_1 - H_2 = W_{12}$$

Dove:

• $H_1$ ed $H_2$ sono i valori del carico totale ai due estremi
  • La differenza $H_1 - H_2$ rappresenta la forza motrice
• $W_{12}$ rappresenta le perdite di carico

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

Il carico totale è dato da:

$$H_i = \frac{p_i}{\rho g} + z + \frac{v_i^2}{2 g}$$

Dove:

• $p$ è la pressione
• $\rho$ è la densità del fluido
• $g$ è l'accelerazione di gravità
• $z$ è la quota
• $v$ è la velocità

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

Per quanto riguarda le perdite di carico:

$$W_{12} = \sum_{i = 1}^{n} \left(W_{i}^{(dist)} + W_{i}^{(conc)}\right)$$

Dove:

• $n$ è il numero di tratte della condotta
• $W_{i}^{(dist)}$ sono le perdite distribuite per la tratta $i$
• $W_{i}^{(conc)}$ sono le perdite concentrate per la tratta $i$

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

Sviluppando i due termini $W_{i}^{(dist)}$ e $W_{i}^{(conc)}$, otteniamo:

$$W_{i}^{(conc)} + W_{i}^{(dist)} = \frac{4f_i}{D_i} \frac{v_i^2}{2g}L_i + \sum_{j = 1}^{m_i} k_{i,j} \frac{v_i^2}{2g}$$

Dove:

• $D_i$ è il diametro idraulico della tratta $i$-ma
• $L_i$ è la lunghezza della tratta $i$-ma
• $v_i$ è la velocità nella tratta $i$-ma
• $k_{i,j}$ sono i coefficienti per le $m$ perdite concentrate
• $f$ è il fattore di attrito e dipende dalla velocità

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

In particolare, per il fattore di attrito sappiamo che:

$$f_i = f\left(\text{Re}_i, \frac{\varepsilon_i}{D_i}\right)$$

Dove:

• $\varepsilon_i$ è la scabrezza nella tratta $i$-ma
• $\text{Re}_i$ è il numero di Reynolds per la tratta $i$-ma

$$\text{Re}_i = \frac{v_i D_i}{\nu}$$

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

Per calcolare $f_i$ possiamo usare per esempio l'equazione di Churchill

$$f_i = \left( \left(\frac{8}{\text{Re}_i}\right)^{12} + \frac{1}{(A_i + B_i)^{^3/_2}} \right)^{^1/_{12}}$$

Con:

$$\begin{align} A_i &= \left( 2.457 \log \left( \frac{1}{(7/\text{Re}_i)^{0.9} + 0.27 (\varepsilon_i / D_i)} \right) \right)^{16} \\ B_i &= \left( \frac{37530}{\text{Re}_i} \right)^{16} \end{align}$$

Comportamento di un Fluido in Condotta

Entrando un po' più nei dettagli, abbiamo che:

Infine sappiamo che la portata $Q$ è data da:

$$Q = \pi \frac{D_i^2}{4} v_i$$

• La portata è costante sulla condotta
• Essa definisce la relazione tra $v_i$ e $D_i^2$ per ogni tratta

Diametro Costante, Determinare $H_1 - H_2$

Partiamo dal caso più semplice:

• Abbiamo una singola tratta (quindi un singolo diametro)
• Vogliamo determinare la forza motrice $H_1 - H_2$
• Tutte le altre informazioni sono note

Diametro Costante, Determinare $H_1 - H_2$

Partiamo dal caso più semplice:

• Abbiamo una singola tratta (quindi un singolo diametro)
• Vogliamo determinare la forza motrice $H_1 - H_2$
• Tutte le altre informazioni sono note

Questo tipo di problema è molto semplice!

• Abbiamo una soluzione in forma chiusa:

$$H_1 - H_2 = \frac{4f}{D} \frac{v^2}{2g}L + \sum_{j = 1}^{m} k_{j} \frac{v^2}{2g}$$

• Si tratta solo di calcolare il valore dell'espressione
• Possiamo usare il calcolatore, ma basta anche una calcolatrice!

Diametro Costante, Determinare $v$ e $Q$

Consideriamo una situazione più interessante:

• Abbiamo una singola tratta (quindi un singolo diametro)
• Vogliamo determinare la portata/la velocità
• Tutte le altre informazioni sono note

Diametro Costante, Determinare $v$ e $Q$

Consideriamo una situazione più interessante:

• Abbiamo una singola tratta (quindi un singolo diametro)
• Vogliamo determinare la portata/la velocità
• Tutte le altre informazioni sono note

L'equazione principale è la stessa di prima:

$$H_1 - H_2 = \frac{4f}{D} \frac{v^2}{2g}L + \sum_{j = 1}^{m} k_{j} \frac{v^2}{2g}$$

• Stavolta però $v$ non è nota e compare nell'espressione di $f$
• Infatti abbiamo che $f = f(\text{Re}, ^\varepsilon/_D)$...
• ...che è non lineare e difficile da trattare

Diametro Costante, Determinare $v$ e $Q$

Possiamo risolverlo utilizzando la IPF

Ci serve una espressione per $v$, che chiameremo $g(f)$:

$$v = g(f) = \sqrt{ \frac{2g(H_1 - H_2)} {\left( 4f\frac{L}{D} + \sum_{j = 1}^m k_j \right)} }$$

• La dipendenza da $v$ è nascosta: $f$ dipende da $\text{Re}$ e quindi da $v$

Diametro Costante, Determinare $v$ e $Q$

Possiamo risolverlo utilizzando la IPF

Ci serve una espressione per $v$, che chiameremo $g(f)$:

$$v = g(f) = \sqrt{ \frac{2g(H_1 - H_2)} {\left( 4f\frac{L}{D} + \sum_{j = 1}^m k_j \right)} }$$

• La dipendenza da $v$ è nascosta: $f$ dipende da $\text{Re}$ e quindi da $v$

Quindi per ogni iterazione di IPF dobbiamo calcolare:

$$v^{(k)} \underbrace{ \xrightarrow{\text{Re} = \frac{vD}{\nu}} \text{Re}^{(k)} \xrightarrow{f = f(\text{Re}, ^\varepsilon/_D)} f^{(k)} \xrightarrow{v = g(f)} }_\text{funzione di transizione} v^{(k+1)}$$

• C'è qualche passo in più, ma è il solito algoritmo!

Diametro Costante, Determinare $v$ e $Q$

Questa volta usare il calcolatore è una buona idea

• Si può anche eseguire l'algoritmo a mano...
• ...ma i calcoli sono complessi...
• ...e può essere necessario effettuare diverse iterazioni

Ci rimane un problema da risolvere:

Da quale valore di $v$ possiamo partire?

• In molti casi partire con $v^{(0)} = 2 m/s$ è una buona idea
• In alternativa, possiamo dedurre $v^{(0)}$ a partire da altri valori

Vediamo qualche esempio...

Diametro Costante, Determinare $v$ e $Q$

Per esempio, possiamo partire da $f$:

• Chiamiamo questo valore di partenza $f^{(-1)}$
• Una buona stima per $f$ in molti casi è $0.005$
• Da $f^{(-1)}$ possiamo calcolare $v^{(0)}$:

$$f^{(-1)} \xrightarrow{v = g(f)} v^{(0)}$$

• Poi continuiamo con lo schema precedente:

$$v^{(k)} \xrightarrow{\text{Re} = \frac{vD}{\nu}} \text{Re}^{(k)} \xrightarrow{f = f(\text{Re}, ^\varepsilon/_D)} f^{(k)} \xrightarrow{v = g(f)} v^{(k+1)}$$

Si può usare lo stesso metodo per altre grandezze

• E.g. partiamo da $\text{Re}$ se abbiamo una idea del tipo di moto

Tratte Multiple, Determinare $v$ e $Q$

Se abbiamo più di una tratta:

La velocità in ogni tratta dipende dal $D_i$, ma la portata è costante

Quindi, conviene concentrarsi sull'ottenere la portata

Tratte Multiple, Determinare $v$ e $Q$

Se abbiamo più di una tratta:

La velocità in ogni tratta dipende dal $D_i$, ma la portata è costante

Quindi, conviene concentrarsi sull'ottenere la portata

Possiamo ottenere una espressione per la portata:

$$Q = h(f_1, f_2, \ldots) = \sqrt{ \frac{ \frac{(H_1 - H_2)\pi^2g}{8} } { \sum_{i = 1}^n \frac{1}{D_i^4} \left( \frac{4 f_i L_i}{D_i} + \sum_{j = 1}^{m_i} k_{i,j} \right) } }$$

• In queste slides la chiameremo $h(f_1, f_2, \ldots)$
• Dove $f_1, f_2, ...$ dipendono a loro volta da $Q$

Tratte Multiple, Determinare $v$ e $Q$

A questo punto possiamo applicare la IPF

Lo schema è il seguente:

$$Q^{(k)} \xrightarrow{v_i = \frac{4Q}{\pi D_i^2}} v_i^{(k)} \xrightarrow{\text{Re}_i = \frac{v_i D_i}{\nu}} \text{Re}_i^{(k)} \xrightarrow{f_i = f(\text{Re}_i, ^\varepsilon/_{D_i})} f_i^{(k)} \xrightarrow{Q = h(f_1, f_2, \ldots)} Q^{(k+1)}$$

Come nel caso precedente:

• Va fatto qualche calcolo in più del solito...
• ...ma l'algoritmo è sempre lo stesso

I calcoli sono tanti e conviene automatizzarli

• Una volta implementato il metodo...
• ...per risolvere un esercizio basterà cambiare i dati noti

Tratte Multiple, Determinare $v$ e $Q$

Per quanto riguarda il valore di partenza:

• Si può andare per tentativi
• Oppure possiamo ottenere $Q^{(0)}$ a partire da altri valori

Per esempio:

• Possiamo partire con $f_0 = f_1 = f_2 = 0.005$
• Quindi calcoliamo $Q^{(0)}$ e procediamo come al solito

Possiamo usare lo stesso metodo con altre grandezze:

• Per esempio $v_i$ oppure $\text{Re}_i$

Tratta Singola, Determinare $D$

Consideriamo un ultimo tipo di problema:

• Abbiamo una singola tratta
• Vogliamo determinare il diametro
• Tutte le altre informazioni sono note

È un esempio di problema di progetto

Tratta Singola, Determinare $D$

Consideriamo un ultimo tipo di problema:

• Abbiamo una singola tratta
• Vogliamo determinare il diametro
• Tutte le altre informazioni sono note

È un esempio di problema di progetto

Partiamo da una espressione per $D$

$$D = d(f) = \left( \frac{32 L f}{\pi^2 g (H_1 - H_2)} Q^2 \right)^{^1/_5}$$

• In queste slides, chiameremo questa funzione $d(f)$
• Come nei casi precedenti $f$ dipende a sua volta da $D$

Tratta Singola, Determinare $D$

A questo punto possiamo usare IPF

Lo schema è il seguente:

$$D^{(k)} \xrightarrow{v = \frac{4Q}{\pi D_i^2}} v^{(k)} \xrightarrow{\text{Re} = \frac{v D}{\nu}} \text{Re}^{(k)} \xrightarrow{f = f(\text{Re}, ^\varepsilon/_{D})} f^{(k)} \xrightarrow{D = d(f)} D^{(k+1)}$$

Per ottenere una stima iniziale possiamo:

• Andare per tentativi
• Partire da $f$ (e.g. $f = 0.005$) ed ottenere $D^{(0)}$
• Partire da $v$ (e.g. $v = 2 ^m/_s$)
• Partire da $\text{Re}$ (se abbiamo una idea sul tipo di moto)

Qualche Osservazione

Possiamo usare IPF per risolvere l'equazione di Bernoulli

• Possiamo eseguirla a mano o implementarla in Octave...
• ...Ma vi consiglio di implementarla in Octave

Per una buona ragione:

• I calcoli sono complessi e ripetitivi: è facile fare errori!
• Implementare costa tempo e un po' di sforzo...
• ...ma una volta che avete un programma funzionante...
• ...è molto facile eseguirlo su dati diversi

Elementi di Informatica e
Applicazioni Numeriche T

Zeri di Equazioni Non-Lineari

Limiti di IPF

• Con IPF possiamo risolvere equazioni non lineari...
• ... E sistemi di equazioni non lineari

Purtroppo IPF ha anche diversi limitazioni:

Una Formulazione Alternativa

Possiamo però formulare lo stesso problema in modo diverso

Data una equazione non lineare, è molto facile portarla nella forma:

$$f(x) = 0$$

• Per esempio, se partiamo da:

$$g(x) = h(x)$$

• Possiamo ottenere facilmente:

$$g(x) - h(x) = 0$$

Una Formulazione Alternativa

Possiamo però formulare lo stesso problema in modo diverso

Se abbiamo un sistema di $m$ equazioni in $n$ incognite

$$\begin{align} g_1(x_1, x_2, \dots x_n) = 0\\ g_2(x_1, x_2, \dots x_n) = 0\\ \ldots\\ g_m(x_1, x_2, \dots x_n) = 0 \end{align}$$

Possiamo vederlo come:

$$f(x) = 0$$

• Con $f : \mathbb{R}^n \rightarrow \mathbb{R}^m$ (da vettore a vettore)

Zero di una Funzione

In sostanza per risolvere un sistema/equazione non lineare...

...Possiamo concentrarci nel trovare una soluzione per:

$$f(x) = 0$$

Questa formulazione offre due grossi vantaggi:

• L'equazione originale viene manipolata in modo molto semplice
• La forma della funzione non dipende dall'incognita

Esempio: Equazione di Bernoulli

Consideriamo l'equazione di Bernoulli (singola tratta)

Possiamo portarla nella forma:

$$(H_1 - H_2) - \frac{4f}{D} \frac{v^2}{2g}L - \sum_{j = 1}^{m} k_{j} \frac{v^2}{2g} = 0$$

• Abbiamo solo dovuto portate tutti i termini su un lato
• Se cerchiamo il valore di $v$ (e quindi $Q$)...
• ...O il valore di $D$...
• ...La forma della funzione resta la stessa

Ora ci serve solo un metodo per trovare lo zero di una funzione!

Elementi di Informatica e
Applicazioni Numeriche T

Metodo della Bisezione

Partiamo da un Caso Semplice

Cominciamo facendo qualche ipotesi semplificativa:

• HP1: $f(x)$ è univariata e scalare:

$$f : \mathbb{R} \rightarrow \mathbb{R}$$

• HP2: $f(x)$ è continua

È un caso semplificato, ma comunque molto utile

• Include tutti le equazioni definite usando polinomi
• Include la maggior parte delle equazioni utili in fisica

Un Esempio

Come primo esempio, supponiamo di voler risolvere:

$$f(x) = x^3 - \frac{3}{2} x^2 + 1 = 0$$

Un Esempio

Il plot della Il plot della funzione, con l'unico zero in rosso

Un Esempio

Come primo esempio, supponiamo di voler risolvere:

$$f(x) = x^3 - \frac{3}{2} x^2 + 1 = 0$$

Come possiamo trovare uno zero?

Partiamo da una osservazione:

• Nel disegno uno zero lo si vede
• Perché $f(x)$ attraversa l'asse delle $x$
• Ossia tra due punti $a$ e $b$, $f(x)$ cambia di segno

Vediamo di essere un po' più formali...

Teorema del Valore Intermedio

Teorema del valore intermedio:

Sia data una funzione $f : \mathbb{R} \rightarrow \mathbb{R}$:

• Se $f$ è continua su un intervallo $[a, b]$
• Allora nell'intervallo $[a, b]$ assumerà tutti i valori tra $f(a)$ e $f(b)$

Corollario importante (teorema di Bolzano):

Sia data una funzione $f : \mathbb{R} \rightarrow \mathbb{R}$ continua su $[a, b]$

• Se $f(a)$ e $f(b)$ hanno segno opposto
• Allora $f(x)$ ha uno zero nell'intervallo $[a, b]$

Intervallo Contenente uno Zero

Quindi, per trovare un intervallo contenente uno zero:

• Possiamo valutare $f(x)$ in vari punti
  • Lo facciamo già per disegnare una funzione!
• Se troviamo due valori $a$ e $b$ tali che $f(a)\, f(b) <0$...
• ...Allora sappiamo che c'è almeno uno zero tra $a$ e $b$

A questo punto, per trovare uno zero:

• Possiamo stringere l'intervallo $[a, b]$...
• ...in modo che continui a contenere uno zero
• Ci fermiamo quando l'intervallo diventa sufficientemente piccolo

Possiamo farlo utilizzando il metodo della bisezione

Metodo della Bisezione

Partiamo da un intervallo $[a, b]$ tale che $f(a) \, f(b) < 0$

Metodo della Bisezione

Otteniamo un nuovo punto $c$ a metà tra $a$ e $b$

Metodo della Bisezione

Aggiorniamo l'intervallo $[a, b]$

Metodo della Bisezione

Ripetiamo il processo

Metodo della Bisezione

Ripetiamo il processo

Metodo della Bisezione

Ripetiamo il processo

Metodo della Bisezione

Ripetiamo il processo

Metodo della Bisezione

Ripetiamo il processo

Metodo della Bisezione

Ci fermiamo quando la distanza tra $a$ (o $b$) e $c$ è abbastanza piccola

Metodo della Bisezione

Le operazioni critiche sono due:

Op 1: Il calcolo del valore di $c$, che è dato da:

$$c = \frac{1}{2}(a+b)$$

Metodo della Bisezione

Le operazioni critiche sono due:

Op 1: Il calcolo del valore di $c$, che è dato da:

$$c = \frac{1}{2}(a+b)$$

Op 2: L'aggiornamento dell'intervallo $[a, b]$

• Se $f(c)$ ha lo stesso segno di $f(a)$, viene rimpiazzato $b$
• Altrimenti viene rimpiazzato $b$
• In questo modo si mantiene la prorietà $f(a) \, f(b) < 0$

Metodo della Bisezione

Possiamo descrivere l'algoritmo mediante pseudo-codice:

for <numero massimo di iterazioni>
    c = (a + b)/2
    if $|a - c|$ < $\varepsilon$
        break
    end
    if <$f(a)$ e $f(c)$ hanno lo stesso segno>
        a = c
    else
        b = c
    end
end

Elementi di Informatica e
Applicazioni Numeriche T

Metodo della Bisezione:
Analisi di Convergenza

Convergenza del Metodo della Bisezione

Il metodo della bisezione ha ottime proprietà di convergenza:

• Converge sempre (entro i limiti di precisione della macchina)
• Non è particolarmente veloce a convergere...
• ...però la velocità di convergenza è molto stabile

Proviamo a caratterizzare meglio la convergenza

Per farlo, dobbiamo:

• Ottenere una espressione per l'errore (i.e. $e = x - x^*$)
• Analizzare il limite:

$$\lim_{k \rightarrow \infty} \frac{\left|e^{(k+1)}\right|}{\left|e^{(k)}\right|}$$

Convergenza del Metodo della Bisezione

• Alla prima iterazione, l'errore è al massimo:

$$e^{(0)} = \frac{1}{2}\left|a - b \right|$$

• Alla seconda iterazione l'intervallo $[a, b]$ si dimezza:

$$e^{(1)} = \frac{1}{2}\frac{\left|a - b \right|}{2}$$

Generalizzando, abbiamo che:

$$e^{(k)} = \frac{1}{2^{k+1}}\left|a - b \right|$$

Velocità di Convergenza

È ora di dare qualche definizione in più sulla convergenza:

Un metodo iterativo converge se:

$$\lim_{k \rightarrow \infty} \left|e^{(k)}\right| = 0$$

Consideriamo poi il rapporto tra errori consecutivi:

$$\lim_{k \rightarrow \infty} \frac{\left|e^{(k+1)}\right|}{\left|e^{(k)}\right|} = \mu$$

• Perché via sia convergenza, deve valere $\mu < 1$

A seconda del valore di $\mu$, distinguiamo diversi tipi di convergenza

Velocità di Convergenza

In particolare abbiamo che:

Caso 1: $0 < \mu < 1$, la convergenza è lineare

• L'errore si riduce di un fattore costante

Caso 2: Se $\mu = 1^-$, la convergenza a sub-lineare

• La convergenza rallenta quando ci si avvicina alla soluzione

Caso 2: Se $\mu = 0$, la convergenza è super-lineare

• La convergenza accelera avvicinandosi alla soluzione

Vediamo come si comporta il metodo della bisezione...

Convergenza del Metodo della Bisezione

Innanzitutto possiamo osservare che il metodo converge:

$$\lim_{k \rightarrow \infty} e^{(k)} = \lim_{k \rightarrow \infty} \frac{1}{2^{k+1}}\left|a - b \right| = 0$$

• Non occorrono ipotesi particolari: il metodo converge sempre

La convergenza è di tipo lineare, perché:

$$\lim_{k \rightarrow \infty} = \frac{\left|e^{(k+1)}\right|}{\left|e^{(k)}\right|} = \frac{1}{2}$$

• L'errore (nel caso peggiore) si dimezza ad ogni iterazione
• Quindi si riduce di un fattore 10 in $\log_2 10 \simeq 3.3$ iterazioni
• Quindi guadagniamo una cifra di precisione ogni circa 3 iterazioni

Convergenza e Zeri Multipli

Se l'intervallo $[a, b]$ contiene più di uno zero

• Il metodo delle bisezione convergerà sempre ad uno zero...
• ...ma non è facile determinare quale

Dipende tutto dai valori scelti per $[a, b]$

Consideriamo una versione modificata del nostro problema:

$$x^3 - \frac{5}{2} x^2 + 1 = 0$$

• La funzione ha tre zeri
• Lo zero di convergenza dipende dalla scelta di $a$ e $b$

Convergenza e Zeri Multipli

In questo plot vediamo i tre zeri:

Convergenza e Zeri Multipli

Con questa scelta di $a$ e $b$, convergiamo allo zero di sx

Convergenza e Zeri Multipli

Con questa scelta di $a$ e $b$, convergiamo allo zero di dx

Convergenza e Zeri Multipli

Se l'intervallo $[a, b]$ contiene più di uno zero

• Il metodo delle bisezione convergerà sempre ad uno zero...
• ...ma non è facile determinare quale

Dipende tutto dai valori scelti per $[a, b]$

In generale "scegliere" verso quale zero convergere è difficile:

• Vale di fatto per tutti i metodi di soluzione
• Ci sono alcune eccezioni, ma non le approfondiremo

Elementi di Informatica e
Applicazioni Numeriche T

Metodo di Newton-Raphson

Limiti del Metodo di Bisezione

Il metodo di bisezione da due grandi vantaggi:

• È molto semplice da implementare
• È molto robusto dal punto di vista della convergenza

Limiti del Metodo di Bisezione

Il metodo di bisezione da due grandi vantaggi:

• È molto semplice da implementare
• È molto robusto dal punto di vista della convergenza

...Un difetto tollerabile (nel nostro caso):

• La convergenza è relativamente lenta

Limiti del Metodo di Bisezione

Il metodo di bisezione da due grandi vantaggi:

• È molto semplice da implementare
• È molto robusto dal punto di vista della convergenza

...Un difetto tollerabile (nel nostro caso):

• La convergenza è relativamente lenta

...E due grossi problemi:

• Richiede di conoscere un intervallo $[a, b]$ che racchiude uno zero
• Funziona solo per funzioni scalari e univariate

Vediamo un secondo metodo di soluzione...

Metodo di Newton-Raphson

Il nostro secondo approccio è il metodo di Newton-Raphson:

Ha caratteristiche complementari rispetto al metodo della bisezione

• Richiede un singolo punto di partenza
• In molti casi converge molto rapidamente...
• ...Ma può anche non convergere

È uno dei metodi più noti per la soluzione di equazioni non lineari

• È applicabile anche a sistemi di equazioni
• Inizialmente però ci focalizziamo su funzioni $f : \mathbb{R} \rightarrow \mathbb{R}$

Metodo di Newton-Raphson

Partiamo da un singolo valore di $x$

Metodo di Newton-Raphson

Approssimiamo la funzione con la tangente in $x^{(0)}$

Metodo di Newton-Raphson

Troviamo uno zero per la funzione approssimata

Metodo di Newton-Raphson

Ripetiamo il processo

Metodo di Newton-Raphson

Ripetiamo il processo

Metodo di Newton-Raphson

Ripetiamo il processo

Metodo di Newton-Raphson

Finché non arriviamo a convergenza

Metodo di Newton-Raphson

L'idea principale del metodo è:

• Approssimare $f(x)$ con il suo sviluppo in serie di Taylor di ordine $1$
• L'operazione richiede che $f(x)$ sia derivabile

All'iterazione $k$-ma abbiamo quindi:

$$f(x) \simeq f(x^{(k)}) + f'(x^{(k)}) (x - x^{(k)})$$

Da cui possiamo ottenere lo zero della funzione approssimata:

$$x = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})}$$

Metodo di Newton-Raphson

L'iterazione fondamentale del metodo è quindi:

$$x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})}$$

Lo pseudo-codice dell'algoritmo è:

$x$ = <stima iniziale>
for <numero massimo di iterazioni>
    $x_p$ = $x$
    $x$ = $x$ - $f(x)$ / $f^\prime(x)$
    if $|x - x_p|$ < $\varepsilon$
        break
    end
end

Un Esempio Numerico

Vediamo un esempio numerico, con:

• $f(x) = x^3 - \frac{5}{2} x^2 + 1$
• E quindi $f'(x) = 3x^2 - 5 x$

Elementi di Informatica e
Applicazioni Numeriche T

Metodo di Newton-Raphson:
Analisi di Convergenza

Convergenza del Metodo di Newton-Raphson

• Il metodo di Newton-Raphson non garantisce la convergenza...
• ...Ma quando converge, lo fa molto velocemente

Proviamo a caratterizzare la convergenza

Partiamo dall'iterazione fondamentale

$$x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})}$$

• È la descrizione di un punto fisso!

$$x = g(x) = x - \frac{f(x)}{f'(x)}$$

Convergenza del Metodo di Newton-Raphson

Possiamo applicare un risultato noto per IPF

• In particolare, in un intorno di $x^*$, abbiamo:

$$\frac{\left|e^{(k+1)}\right|}{\left|e^{(k)}\right|} = \left|g'(x)\right|$$

• Quindi la convergenza dipende dal valore di $\left|g'(x)\right|$:

$$g'(x) = \left(x - \frac{f(x)}{f'(x)}\right)' = 1 - \frac{f'(x)^2 - f(x) f''(x)}{f'(x)^2} = \frac{f(x) f''(x)}{f'(x)^2}$$

Convergenza del Metodo di Newton-Raphson

Per avere convergenza, in un intorno di $x^*$ deve valere:

$$\left|g'(x)\right| = \left|\frac{f(x) f''(x)}{f'(x)^2}\right| = \frac{\left|f(x)\right| \left|f''(x)\right|}{f'(x)^2} < 1$$

A questo punto ci limitiamo a delle considerazioni informali

Considerazione 1:

• Poiché $f(x^*)$ vale $0$...
• ...Se $f$ è continua il valore di $\left|f(x)\right|$ sarà molto piccolo
• Quindi se $f$ è continua è probabile che il metodo converga

Convergenza del Metodo di Newton-Raphson

Per avere convergenza, in un intorno di $x^*$ deve valere:

$$\left|g'(x)\right| = \left|\frac{f(x) f''(x)}{f'(x)^2}\right| = \frac{\left|f(x)\right| \left|f''(x)\right|}{f'(x)^2} < 1$$

A questo punto ci limitiamo a delle considerazioni informali

Considerazione 2:

• Se $f'(x^*) \neq 0$, $f''(x^*)$ è limitata e $f$ è continua, allora:

$$\lim_{k \rightarrow \infty} \frac{\left|e^{(k+1)}\right|}{\left|e^{(k)}\right|} = \lim_{k \rightarrow \infty} \left|g'(x^{(k)})\right| = 0$$

• La convergenza è super-lineare!
• In particolare: il numero di cifre corrette raddoppia ad ogni iterazione

Convergenza del Metodo di Newton-Raphson

Per avere convergenza, in un intorno di $x^*$ deve valere:

$$\left|g'(x)\right| = \left|\frac{f(x) f''(x)}{f'(x)^2}\right| = \frac{\left|f(x)\right| \left|f''(x)\right|}{f'(x)^2} < 1$$

A questo punto ci limitiamo a delle considerazioni informali

Considerazione 3:

• Se $f'(x^*) = 0$ (lo zero è su un estremo), il denominatore è piccolo
• In linea di principio il metodo potrebbe divergere
• In pratica, se $f'$ è continua la convergenza c'è, ma è solo lineare

Convergenza del Metodo di Newton-Raphson

Per avere convergenza, in un intorno di $x^*$ deve valere:

$$\left|g'(x)\right| = \left|\frac{f(x) f''(x)}{f'(x)^2}\right| = \frac{\left|f(x)\right| \left|f''(x)\right|}{f'(x)^2} < 1$$

A questo punto ci limitiamo a delle considerazioni informali

Considerazione 4:

• Se $f''(x^*) =$ (lo zero è su un flesso), il numeratore sarà piccolo
• La convergenza sarà ancora più veloce

Convergenza del Metodo di Newton-Raphson

Per avere convergenza, in un intorno di $x^*$ deve valere:

$$\left|g'(x)\right| = \left|\frac{f(x) f''(x)}{f'(x)^2}\right| = \frac{\left|f(x)\right| \left|f''(x)\right|}{f'(x)^2} < 1$$

A questo punto ci limitiamo a delle considerazioni informali

In generale:

• Caratterizzare i casi di convergenza del metodo è difficile
• Il più delle volte il metodo converge
• I casi di divergenza capitano perlopiù se $f'$ non è continua
• Se per un dato $x^{(k)}$ abbiamo $f'(x^{(k)}) = 0$, il metodo non funziona

Elementi di Informatica e
Applicazioni Numeriche T

Metodo di Newton-Raphson:
Derivata Non Nota

Difficoltà nell'Ottenre $f'(x)$

Il metodo di Newton-Raphson richiede di conoscere la derivata di $f$

A volte però $f'(x)$ è difficile da ottenere per via analitica

• E.g. la derivata in $v$ dell'equazioni di Bernoulli:

$$(H_1 - H_2) - \frac{4f}{D} \frac{v^2}{2g}L - \sum_{j = 1}^{m} k_{j} \frac{v^2}{2g} = 0$$

• Ricordate che il fattore di attrito dipende da $v$!
• Ottenere una espressione per la derivata in questo caso è difficile

Per fortuna, si può evitare di farlo

Approssimazione del Valore di $f'(x)$

Partiamo da una osservazione:

• Non ci interessa avere une espressione per $f'(x)$...
• ...Ma solo poterla calcolare in punti specifici

Per definizione, abbiamo che:

$$f'(x) = \lim_{\delta \rightarrow 0} \frac{f(x+\delta) - f(x)}{\delta}$$

Quindi possiamo calcolare una approssimazione di $f'(x)$ con:

$$f'(x) \simeq \frac{f(x+\varepsilon) - f(x)}{\varepsilon}$$

• Con $\varepsilon$ sufficientemente piccolo

Approssimazione del Valore di $f'(x)$

Si tratta di un calcolo approssimato

• Ma in fondo il calcolatore ha una precisione finita...
• ...quindi anche se avessimo una espressione per $f'(x)$...
• ...al momento di calcolarla avremmo una approssimazione!

Inoltre, possiamo sempre migliorare l'approssimazione:

$$f'(x) \simeq \frac{f(x+\varepsilon) - f(x)}{\varepsilon}$$

• Basta prendere un $\varepsilon$ molto piccolo, no?
• In teoria, sì. In pratica, è una brutta idea

Vediamo perché...

Errori di Cancellazione

Il problema è legato all'utilizzo dei numeri in virgola mobile

Supponiamo di avere:

• $\varepsilon = 0.001$
• $f(x) = 1,000,000,000$
• $f(x + \varepsilon) = 1,000,000,001$

E di rappresentare i numeri in virgola mobile, con:

• Tre cifre decimali (più il segno) per la mantissa
• Una cifra (più il segno) per l'esponente

Per semplicità, lavoriamo con cifre decimali (anziché con i bit)

Errori di Cancellazione

Con le ipotesi specificate, abbiamo:

• $\varepsilon = 0.001 = 1.000 \times 10^{-3}$
• $f(x) = 1,000,000,000 = 1.000 \times 10^{9}$
• $f(x + \varepsilon) = 1,000,000,001 = \underline{1.000 \times 10^{9}}$ [e qui butta male]

E quindi:

$$\frac{f(x+\varepsilon) - f(x)}{\varepsilon} = \frac{1.000 \times 10^{9} - 1.000 \times 10^{9}}{1.000 \times 10^{-3}} = 0$$

Mentre il vero risultato è

$$\frac{f(x+\varepsilon) - f(x)}{\varepsilon} = \frac{1}{0.001} = 1,000$$

Errori di Cancellazione

Questo è un esempio di errore di cancellazione

• Una sottrazione tra due numeri in virgola mobile...
• Che dovrebbe essere non-zero...
• Restituisce invece un risultato nullo

Nel nostro caso:

• Usando i numeri in virgola mobile...
• $f(x+\varepsilon)$ e $f(x)$ hanno la stessa rappresentazione!
• Il problema è dato dal fatto che i due numeri sono molto grandi
• Ma differiscono per un valore molto piccolo

Errori di Cancellazione

Vediamo un altro esempio classico:

• $a = 1,000,000,000$ (molto grande)
• $b = 1,000$ (piccolo rispetto ad $a$)
• E quindi: $a + b = 1,000,001,000$

Invece, con le stesse assunzioni di prima, abbiamo:

• $a = 1.000 \times 10^9$
• $b = 1.000 \times 10^3$
• E quindi: $a + b = 1.000 \times 10^9 + \underline{0.000 \times 10^9} = 1.000 \times 10^9$

Perché per poter sommare i numeri, l'esponente deve essere identico

Errori di Cancellazione

In generale, bisogna stare in guardia ogni volta che:

• Vogliamo sommare/sottrarre numeri grandi con numeri piccoli
• Vogliamo distinguere differenze piccole tra numeri grandi

In generale:

• Ogni volta che vogliamo "mischiare" livelli di precisione diversa...
• È possibile che sorgano errori di cancellazione

Infatti, i numeri in virgola mobile:

• Sono accurati nel rappresentare numeri piccoli...
• Ma non accurati con i numeri grandi!

Calcolo "Sicuro" della Derivata Approssimata

Tornando alla nostra derivata approssimata:

$$\frac{f(x+\varepsilon) - f(x)}{\varepsilon}$$

Esiste una buona regola per ottenere...

• Una approssimazione ragionevolmente accurata
• Degli errori di cancellazione limitati

In particolare:

• Se $\epsilon$ è il più piccolo numero rappresentabile sul calcolatore
• Allora un buon valore per $\varepsilon$ è $\sqrt{\epsilon}\,x$...
• O semplicemente $\sqrt{\epsilon}$ se $x = 0$

Elementi di Informatica e
Applicazioni Numeriche T

Metodo di Newton-Raphson:
Caso Multidimensionale

Sistemi di Equazioni Non Lineari

• Ci siamo concentrati su equazioni in una singola variabile
• Equivalentemente: zeri di funzioni $f : \mathbb{R} \rightarrow \mathbb{R}$

Però, Il metodo di N-R funziona anche per sistemi non lineari

• È uno dei grandi vantaggi del metodo rispetto alla bisezione
• Consideriamo solo sistemi determinati ($n$ eq., $n$ incognite)

L'idea di base resta la stessa:

• Utilizziamo una approssimazione lineare di $f(x)$
• E applichiamo la IPF fino a convergenza

Questa volta però abbiamo $f : \mathbb{R}^n \rightarrow \mathbb{R}^n$

Newton-Raphson per Sistemi di Non Lineari

Formalmente, dobbiamo risolvere:

$$f(x) = 0$$

• Con $f : \mathbb{R}^n \rightarrow \mathbb{R}^n$

Possiamo rappresentare l'equazione vettoriale nella forma:

$$\begin{align} f_1(x_1, x_2, \ldots, x_n) &= 0\\ f_2(x_1, x_2, \ldots, x_n) &= 0\\ \ldots\\ f_n(x_1, x_2, \ldots, x_n) &= 0 \end{align}$$

• Dove $f_i(x)$ è la $i$-ma componente di $f$

Newton-Raphson per Sistemi di Non Lineari

Applichiamo una approssimazione lineare ad ogni componente:

$$\begin{align} f_1(x^{(k)}) + \nabla f_1(x^{(k)}) (x - x^{(k)}) &= 0\\ f_2(x^{(k)}) + \nabla f_2(x^{(k)}) (x - x^{(k)}) &= 0\\ \ldots\\ f_n(x^{(k)}) + \nabla f_n(x^{(k)}) (x - x^{(k)}) &= 0 \end{align}$$

Dove:

• $x$ è il vettore delle variabili $x_1, x_2, \ldots$
• $x^{(k)}$ è un vettore con i valori di $x_1^{(k)}, x_2^{(k)}, \ldots$
• $\nabla f_1(x^{(k)})$ è il gradiente di $f_i$ calcolato in $x^{(k)}$

Newton-Raphson per Sistemi di Non Lineari

In forma matriciale, abbiamo:

$$f(x^{(k)}) + J_f(x^{(k)}) (x - x^{(k)}) = 0$$

Dove:

• $f(x^{(k)})$ è il vettore con $f_1(x^{(k)}), f_2(x^{(k)}), \ldots$
• $J_f(x^{(k)})$ è lo Jacobiano di $f$ calcolato in $x^{(k)}$
• La riga $i$-ma dello Jacobiano coincide con il gradiente $\nabla f_i(x)$

Lo Jacobiano si può calcolare con le stesse tecniche viste per $f'(x)$:

• O otteniamo una espressione per tutti i termini dei i gradienti
• Oppure possiamo approssimarli con un rapporto incrementale

Newton-Raphson per Sistemi di Non Lineari

Di fatto, abbiamo un normale sistema lineare:

$$f(x^{(k)}) + J_f(x^{(k)}) (x - x^{(k)}) = 0$$

• Lo possiamo risolvere come al solito
• E.g. con la riduzione di Gauss, o la divisione sinistra in Octave

La soluzione è data da:

$$x^{(k+1)} = x^{(k)} - J_f^{-1}(x^{(k)}) f(x^{(k)})$$

• Dove l'inverso dello Jacobiano fa le veci di $^1/_{f'(x)}$ nel caso scalare
• Abbiamo così ottenuto l'iterazione fondamentale del metodo N-R...
• ...Nel caso multidimensionale

Newton-Raphson per Sistemi di Non Lineari

Il resto dell'algoritmo rimane pressoché invariato:

$x$ = <stima iniziale>
for <numero massimo di iterazioni>
    $x_p$ = $x$
    $x$ = $x$ - $J_f(x)$ / $f(x)$ % divisione sinistra
    if $\|x - x_p\|$ < $\varepsilon$
        break
    end
end

• L'unica accortezza è usare la norma Euclidea $\|\cdot\|$ invece di $|\cdot|$...
• Perché $x$ e $x_p$ sono vettori

Elementi di Informatica e
Applicazioni Numeriche T

Considerazioni Finali

Considerazioni Finali

Abbiamo visto due metodi per risolvere equazioni non lineari

Entrambi trovano uno zero di una funzione

• È più facile formulare un problema in questi termini...
• ...Piuttosto che usare IPF

Detto questo IPF rimane un po' più flessibile:

• Per esempio, il metodo di Newton-Raphson...
• ...Non è altro che un caso particolare di IPF!

Considerazioni Finali

I due metodi visti hanno caratteristiche complementari:

Il metodo della bisezione ha una convergenza molto robusta...

• ...Ma richiede più informazioni, è più lento...
• ...E non è applicabile a sistemi di equazioni

Il metodo di Newton-Raphson:

• Converge molto rapidamente ed è applicabile a sistemi...
• ...Ma la convergenza non è garantita

In pratica, si utilizza spesso una combinazione di metodi

• Octave fornisce delle funzioni di questo tipo!