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