sqrt(x+a) - sqrt(x)的数值稳定计算

Question

对于全参数范围x，一个>= 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