Calcolo numerico/Introduzione all'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 matematico 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. È 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. Su un calcolatore, algoritmi diversi portano però a risultati diversi. L'ordine con cui vengono svolte 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

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

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

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

Quando $a$ è positivo ma molto piccolo, il secondo procedimento è stabile mentre l'altro no, iniziano a esserci problemi per numeri dell'ordine di $10^{-16}$.

Esempio

Consideriamo la risoluzione del sistema lineare

$A𝒙=𝒃$

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

$\frac{det{A}_j}{detA}$

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

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

e

$A^j=[A^1,A^{j-1},𝒃,A^{j+1},A^n$

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

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

$C_2=2$

$C_n\ge n\ast c_{n-1}\ge n\ast (n-1)\ast c_{n-2}\ge \dots \ge n\ast (n-1)\ast \dots \ast 3\ast C_2\ge 2n!\ge n!$

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

Il metodo di Gauss trasforma il sistema in un sistema equivalente della forma $U(𝒙)=𝒄$ con $U$ matrice triangolare superiore.

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

Nell'algoritmo si combinano linearmente le due equazioni.

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

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

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

Il risultato di questa tecnica è che $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.

Esempio 1.1

Supponiamo di avere il sistema

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

e le due rette si intersecano nel punto $1x=y=1$.

Presa la prima equazione, se perturbiamo il coefficiente di $x$ ponendolo uguale a $1+1/100$ otteniamo:

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

e il punto di intersezione è

$\overline{x}=1/9,\overline{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. Consideriamo ad esempio una funzione di classe almeno $C^2$, e consideriamo

$f(x+h)=f(x)+h{f}^{\prime }(x)+h^2/2f^{(2)}(\xi ),x<\xi <x+h,\text{formula di Taylor}$

$f^{\prime }(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 )$, e va a 0 come $h$.

2. Metodo iterativo: ad esempio, vogliamo approssimare la soluzione dell'equazione $x=-\cos x$. Definisco una successione $x_{i+1}=\cos x_i$ con $i>0$, assegnando come valore iniziale $x_0=1$, e proseguendo così.