Introduzione sull'aritmetica floating-point

Quattro parole chiave dell'aritmetica floating-point

Stabilità dell'algoritmo

Un algoritmo è l'implementazione del metodo che vogliamo usare per risolvere il problema matmatico che vogliamo trattare. La stabilità è in corrispondenza con l'errore algoritmico.

Condizionamento di un problema

A cui viene associato l'errore inerente (tipico del problema).

Complessità computazionale

Riguarda il tempo impiegato da un algoritmo per ottenere la soluzione. E' relazionata con il tempo di calcolo e con l'errore algoritmico.

Convergenza

associata a errore analitico o di approssimazione.

Stabilità di un algoritmo

Un algoritmo è una sequenza finita di operazioni. Per uno stesso problema possono esistere diversi algoritmi che portano alla soluzione. Algoritmi diversi su un calcolatore portano però a risultati diversi. L'ordine con cui vengono fatte le operazioni può portare a risultati leggermente diversi.

Si definiscono algoritmi numericamente instabili quelli in cui si ha un'elevata propagazione dell'errore, mentre si definiscono algoritmi numericamente stabili quelli in cui la propagazione dell'errore è contenuta.16 è inoltre il numero di cifre decimali rappresentabili con matlab.

Esempio

Voglio calcolare una quantità $b=1+a-1$$b=1+a-1$ (in aritmetica esatta $b=a$$b=a$).

Algoritmo 1: $(a+1)-1$$(a+1)-1$

Algoritmo 2: $(1-1)+a$$(1-1)+a$ Quando $a$$a$ è positivo ma molto piccolo, il secondo procedimento è stabile mentre l'altro no (cfr stabilità.txt), iniziano a esserci problemi per numeri dell'ordine di $10^{-16}$$10^{-16}$.

 

Esempio

Consideriamo la risoluzione del sistema lineare

$$A{\boldsymbol {x}}={\boldsymbol {b}}$$

$$A{\boldsymbol {x}}={\boldsymbol {b}}$$

con $A$$A$ matrice quadrata $n\times n$$n\times n$. Il sistema si può risolvere con il Metodo di Kramer: la soluzione $x_{j}$$x_j$ del sistema può essere calcolata come il rapporto tra due determinanti

$${\frac {\det A_{j}}{\det A}}$$

$${\frac {\det A_{j}}{\det A}}$$

dove, indicando con $A^{i}$$A^{i}$ le colonne della matrice,

$$A=[A^{1},A^{2},\dots ,A^{n}]$$

$$A=[A^{1},A^{2},\dots ,A^{n}]$$

e

$$A^{j}=[A^{1},A^{j-1},{\boldsymbol {b}},A^{j+1},A^{n}$$

$$A^{j}=[A^{1},A^{j-1},{\boldsymbol {b}},A^{j+1},A^{n}$$

dove sostituisco alla j-esima colonna il termine noto. Il determinante si calcola con la formula dei cofattori. Chiamiamo $c_{n}$$c_n$ il numero di operazioni moltiplicative per calcolare il determinante di una matrice $A$$A$ di dimensione $n$$n$.

Il calcolo del determinante è ridotto al calcolo di n determinanti di matrici più piccole ($(n-1)\times (n-1)$$(n-1)\times (n-1)$) moltiplicate per coefficienti. Quindi $c_{n}=n*C_{n-1}+n$$c_{n}=n*C_{n-1}+n$ (operazione ricorsiva).

$$C_{2}=2$$

$$C_{2}=2$$

$$C_{n}\geq n*c_{n-1}\geq n*(n-1)*c_{n-2}\geq \dots \geq n*(n-1)*\dots *3*C_{2}\geq 2n!\geq n!$$

$$C_{n}\geq n*c_{n-1}\geq n*(n-1)*c_{n-2}\geq \dots \geq n*(n-1)*\dots *3*C_{2}\geq 2n!\geq n!$$

Per trovare tutte le soluzioni del sistema devo calcolare $n+1$$n+1$ determinanti di matrici di dimensione $n$$n$, quindi il costo del metodo di Kramer è maggiore uguale di $(n+1)n!=(n+1)!$$(n+1)n!=(n+1)!$ e richiede tempi molto lunghi. Il costo del metodo di Gauss è invece maggiore o uguale di $n^{3}/3$$n^{3}/3$.

Il metodo di Gauss trasforma il sistema in un sistema equivalente della forma $U({\boldsymbol {x}})={\boldsymbol {c}}$$U({\boldsymbol {x}})={\boldsymbol {c}}$ con U matrice triangolare superiore.

Calcolo un coefficiente $m_{k}$$m_{k}$ detto moltiplicatore, tale che $m_{k}={\frac {A_{ik}}{A_{kk}}}$$m_{k}={\frac {A_{ik}}{A_{kk}}}$ al passo $k-1$$k-1$.

Nell'algoritmo si combinano linearmente le due equazioni.

L'elemento $a_{kk}$$a_{kk}$ è detto elemento pivetale. Se $a_{kk}=0$$a_{kk}=0$ in aritmetica esatta si ha un problema, perché non si può dividere per 0. In aritmetica floating-point quando $a_{kk}$$a_{kk}$ è piccolo sorgono dei problemi.

Con il calcolatore si fa in modo che l'$a_{kk}$$a_{kk}$ che viene usato sia il meno piccolo possibile. Applico la tecnica del pivet parziale: cerco l'indice di equazione ${\bar {i}}$${\bar {i}}$ tale che $a_{{\bar {i}}k}$$a_{{\bar {i}}k}$ in modulo è il massimo per $j=k,\dots ,n$$j=k,\dots ,n$ dei valori assoluti di $a_{jk}$$a_{jk}$.

Prima di proseguire si scambia la riga $k$$k$-esima con la ${\bar {i}}$${\bar {i}}$-esima.

Il risultato di questa tecnica è che $m_{ik}={\frac {a_{ik}}{a_{kk}}}<1$$m_{ik}={\frac {a_{ik}}{a_{kk}}}<1$

 

Condizionamento

Esistono problemi in cui si ha un errore elevato indipendentemente dal metodo utilizzato per calcolare la soluzione. La difficoltà è una proprietà intrinseca del problema stesso.

Definizione 1.1

Si dice problema mal condizionato un problema in cui piccole variazioni dei dati inducono grandi variazioni nei risultati (dipendenza continua dai dati).

 

Esempio 1.1

Supponiamo di avere il sistema

$${\begin{cases}x+y=2\\1001x+1000y=2001\end{cases}}$$

$${\begin{cases}x+y=2\\1001x+1000y=2001\end{cases}}$$

e le due rette si intersecano nel punto $x=y=1$$x=y=1$.
Presa la prima equazione, se perturbiamo il coefficiente di x ponendolo uguale a $1+1/100$$1+1/100$ otteniamo:

$${\begin{cases}(1+1/100)x+y=2\\1001x+1000y=2001\end{cases}}$$

$${\begin{cases}(1+1/100)x+y=2\\1001x+1000y=2001\end{cases}}$$

e il punto di intersezione è

$${\bar {x}}=1/9,\,{\bar {y}}={\frac {1901}{900}},$$

$$\bar x = 1/9, \, \bar y = \frac{1901}{900},$$

cioè una piccola variazione dei dati da una grande variazione nella soluzione. (questo avviene perché le due rette sono quasi parallele).

 

Convergenza

Si ha a che fare con la convergenza in due casi:

1. processi di discretizzazione. Considero ad esempio una funzione di classe almeno $C^{2}$$C^2$, e consideriamo
   $$f(x+h)=f(x)+hf'(x)+h^{2}/2f^{(2)}(\xi ),\,x<\xi <x+h,\quad {\hbox{formula di Taylor}}$$

   $$f(x+h)=f(x)+hf'(x)+h^{2}/2f^{(2)}(\xi ),\,x<\xi <x+h,\quad {\hbox{formula di Taylor}}$$
   $$f'(x)={\frac {f(x+h)-f(x)}{h}}-h/2f^{(2)}(\xi ).$$

   $$f'(x)={\frac {f(x+h)-f(x)}{h}}-h/2f^{(2)}(\xi ).$$
   La discretizzazione ha un errore che si comporta come $h/2f^{(2)}(\xi )$$h/2f^{(2)}(\xi )$, e va a 0 come h.
2. metodo iterativo: ad esempio, voglio approssimare la soluzione dell'equazione $x=-\cos x$$x=-\cos x$. Definisco una successione $x_{i+1}=\cos x_{i}$$x_{i+1}=\cos x_{i}$ con $i>0$$i>0$, assegnando come valore iniziale $x_{0}=1$$x_{0}=1$, e proseguendo così.