Valutazione numericamente stabile di sqrt (x + a) - sqrt (x)

Question

Esiste un modo elegante di valutare numericamente la seguente espressione per l'intero intervallo di parametri x, a> = 0?

f(x,a) = sqrt(x+a) - sqrt(x) 

è Inoltre vi qualsiasi linguaggio di programmazione o una libreria che non forniscono questo tipo di funzione? Se sì, sotto quale nome? Non ho alcun problema specifico con l'espressione sopra, ma lo ho incontrato molte volte in passato e ho sempre pensato che questo problema doveva essere risolto prima!

Answer 1

Sì, c'è! A condizione che almeno uno di x e a è positivo, è possibile utilizzare:

f(x, a) = a/(sqrt(x + a) + sqrt(x)) 

che è perfettamente numericamente stabile, ma non vale la pena una funzione di libreria a sé stante. Naturalmente, quando x = a = 0, il risultato dovrebbe essere 0.

Spiegazione: sqrt(x + a) - sqrt(x) è uguale a (sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x))/(sqrt(x + a) + sqrt(x)). Ora moltiplica i primi due termini per ottenere sqrt(x+a)^2 - sqrt(x)^2, che si semplifica a a.

Ecco un esempio che illustra la stabilità: il caso fastidioso per l'espressione originale è dove x + a e x sono molto vicine in termini di valore (o equivalentemente quando a è molto più piccolo in grandezza rispetto x). Ad esempio, se x = 1 e a sono di piccole dimensioni, sappiamo da un'espansione di Taylor intorno a che sqrt(1 + a) deve essere , quindi sqrt(1 + a) - sqrt(1) dovrebbe essere vicino a a/2 - a^2/8. Proviamo questo per una particolare scelta di piccoli a. Ecco la funzione originale (scritto in Python, in questo caso, ma è possibile trattarlo come pseudocodice):

def f(x, a): 
    return sqrt(x + a) - sqrt(x) 

ed ecco la versione stabile:

def g(x, a): 
    if a == 0: 
     return 0.0 
    else: 
     return a/((sqrt(x + a) + sqrt(x)) 

Ora vediamo cosa otteniamo con x = 1 e a = 2e-10:

>>> a = 2e-10 
>>> f(1, a) 
1.000000082740371e-10 
>>> g(1, a) 
9.999999999500001e-11 

Il valore che dovremmo abbiamo ottenuto è (fino a precisione della macchina): a/2 - a^2/8 - per questo particolare a , i termini di ordine cubico e superiore sono insignificanti nel contesto dei galleggianti a precisione doppia IEEE 754, che forniscono solo circa 16 cifre decimali di precisione. Facciamo calcolare tale valore per il confronto:

>>> a/2 - a**2/8 
9.999999999500001e-11