[EX] Differenze finite e Roundoff error

Ciao a tutti,

propongo un esericizio/tematica che viene data spesso per "random", ma che in realtà non lo è. Considerata l'approssimazione con differenze finite della derivata prima di una funzione $f$ nel punto $x$, si ha

$$\frac{f(x+h)-f(x-h)}{2h} = f'(x) + \frac{f^{(3)}(x)}{6} h^2 + \mathcal{O}(h^4)$$

dove $h$ è l'ampiezza del passo spaziale della discretizzazione. E' evidente che questa approssimazione è del secondo ordine, tuttavia, plottando in un grafico logaritmico il valore di $h$ contro l'approssimazione della derivata prima di $\cos(x)$ a $x=1$, si vede che per valori molto piccoli di $h$ le cose vanno male (il codice Python usato per produrre il grafico è in spoiler)

Testo nascosto, fai click qui per vederlo

Codice:

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return np.cos(x)

def df(x):
    return -np.sin(x)

def FD(f,xbar,h):
    return (f(xbar+h) - f(xbar-h))/(2.0*h)



hrange = np.zeros([15])
N = np.size(hrange)
h = 1

for i in range(N):
    hrange[i]= h/(10.0)**i

print(hrange)
err = []
xbar = 1.0

for k in range(N):
    err.append(np.abs( FD(f,xbar, hrange[k]) - df(xbar)) )
   
plt.figure(figsize=(10,5))
plt.title('Errore |FD(cos(x),1,h) - (-sin(1))| come funzione di $h$')
plt.loglog(hrange,err,'b-o',markerfacecolor='Orange',linewidth='2')
plt.xlabel('h')
plt.ylabel('error')

Come si vede, da $h<10^{-5}$ in poi si ha che l'accuratezza mano a mano viene persa : si dice, anche nei libri, che ciò è dovuto agli errori di arrotondamento, senza dilungarsi troppo.

La domanda dunque è la seguente: è possibile, invece, dimostrare che $h \approx 10^{-5}$ è il valore dal quale in poi si perde accuratezza, assumendo di lavorare in precisione doppia (cioè con $\mu =2^{-53} \approx 10^{-16}$ unità di roundoff)?


Si può assumere che l'approssimazione di $f(x)$ in aritmetica di macchina soddisfi $$\text{fl}(f(x)) = f(x)(1+\delta) \text{ con }|\delta| \leq \mu$$

Soluzione in spoiler

Testo nascosto, fai click qui per vederlo

Detta $g(x)$ l'approssimazione centrata con diff. finite della derivata prima descritta nel primo post, nei libri di solito si stima sempre $|f'(x)-g(x)|$. In realtà, considerando gli errori di arrotondamento serve introdurre $\tilde{g(x)}$, la rappresentazione floating point di tale rapporto incrementale. In particolare,

$$|g(x) - \tilde{g(x)}| = |\frac{f(x+h)-f(x-h) - f(x+h)(1+\delta) + f(x-h) (1+\delta)}{2h}|$$ dove ho usato il fatto che $\tilde{g(x)}$ è fatta solamente da valutazioni di funzioni, ricordando che $|\delta| < \mu$.

Ora, serve notare che $$|f(x+h)-f(x+h)(1+\delta)| \leq \mu$$ (stessa cosa per $f(x-h)$) e perciò:

$$|g(x) - \tilde{g(x)}| \leq \frac{2 \mu}{2h} = \frac{\mu}{h}$$

A questo punto, il **vero** errore commesso è $$|f'(x)-\tilde{g}(x)| = |f'(x)-g(x)+g(x)-\tilde{g(x)}| \leq |f'(x)-g(x)| + |g(x)-\tilde{g(x)}| \leq M h^2 + \frac{\mu}{h}$$

avendo indicato con $M=\max \frac{f^{(3)}(x)}{6}$, sull'intervallo di interesse, assumendo $f$ abb. regolare. In questo caso $|M| \leq 1$.

Già da qui si può notare il mantra "per $h$ troppo piccolo l'errore di arrotondamento è dominante". Di fatto, $\frac{1}{h}$ si impenna e la stima perde di significato, mentre nella pratica si osserva una perdita di accuratezza.

Si può però notare che $\partial_h ( \frac{1}{h^2} + \frac{\mu}{h} ) = 2h - \frac{\mu}{h^2}$ da cui segue che $$h \approx (\frac{\mu}{2})^{1/3} \approx 10^{-5}$$
($\mu = 2^{-53}$ in doppia precisione)

conferamando quanto osservato nel grafico sopra.