math - 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 .

这是一个证明稳定性的示例:原始表达式的麻烦情况是 where 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 double 浮点数的上下文中无关紧要，它仅提供大约 16 位十进制数字的精度。让我们计算该值进行比较:

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