數學穩定評估sqrt（x + a） - sqrt（x）

Question

是否存在一個優化的方法來在整個參數範圍x中進行數值穩定評估以下表達式：a> = 0？

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

也有任何編程語言或庫，它提供這種功能？如果是，以什麼名字？我現在沒有使用上述表達式的具體問題，但在過去遇到過很多次，並一直認爲這個問題一定是以前解決的！

Answer 1

是的，有！前提是x至少一個a是肯定的，你可以使用：

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

這是完全數值穩定，但不值得在自己的權利的庫函數。當然，當x = a = 0，結果應該是0。

說明：sqrt(x + a) - sqrt(x)等於(sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x))/(sqrt(x + a) + sqrt(x))。現在乘以前兩項得到sqrt(x+a)^2 - sqrt(x)^2，簡化爲a。

下面是一個例子展示了穩定性：對於原始表達式的麻煩的情況是x + a和x在數值上非常接近（或等效地當a在量值上比x小得多）。例如，如果x = 1和a很小，我們從泰勒擴展1左右知道sqrt(1 + a)應該是1 + a/2 - a^2/8 + O(a^3)，所以sqrt(1 + a) - sqrt(1)應該接近a/2 - a^2/8。讓我們嘗試一個小的a的特定選擇。這裏的原始功能（用Python編寫的，在這種情況下，但你可以把它當作僞代碼）：

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

和這裏的穩定版本：

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

現在讓我們看看我們得到與x = 1和a = 2e-10：

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

我們應該得到的值（最多機牀精度）：a/2 - a^2/8 - 這個特殊a在IEEE 754雙精度浮點數的上下文中，三次和更高階的條件是微不足道的，它只提供大約16個十進制數字的精度。讓我們來計算該值進行比較：

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