浮点运算错误

我使用以下函数来近似函数在某一点的导数：

def prime_x(f, x, h):

    if not f(x+h) == f(x) and not h == 0.0: 
        return (f(x+h) - f(x)) / h
    else:
        raise PrecisionError

作为一个测试我正在通过f as fx and x如3.0。在哪里fx is:

def fx(x):

    import math
    return math.exp(x)*math.sin(x)

其中有exp(x)*(sin(x)+cos(x))作为衍生品。现在，根据谷歌和我的计算器

exp(3)*(sin(3)+cos(3)) = -17.050059.

到目前为止，一切都很好。但是当我决定用小值测试该函数时h我得到以下信息：

print prime_x(fx, 3.0, 10**-5)
-17.0502585578
print prime_x(fx, 3.0, 10**-10)
-17.0500591423
 print prime_x(fx, 3.0, 10**-12)
-17.0512493014
print prime_x(fx, 3.0, 10**-13)
-17.0352620898
print prime_x(fx, 3.0, 10**-16)
__main__.PrecisionError: Mantissa is 16 digits

为什么当 h 减小时（在某个点之后）误差会增大？我本来期待相反的情况，直到f(x+h)等于f(x).

Floating-point arithmetic (and integer arithmetic and fixed-point arithmetic) has a certain granularity: Values can only change by a certain step size. For IEEE-754 64-bit binary format, that step size is about 2^–52 times the value (about 2.22•10^–16). That is very small for physical measurements.

然而，当你做h非常小，f(x) and f(x+h）与步长相比并不是很大。差值只能是步长的整数倍。

当导数为d，f(x）是关于h•d。即使你计算 f(x) and f(x+h）以及可能的浮点格式，它们的差值的测量值必须是步长的倍数s，所以它必须是圆形的（h•d/s)•s，其中圆（y) is y四舍五入到最接近的整数。显然，正如你所做的那样h更小，h•d/s较小，因此四舍五入为整数的效果相对较大。

另一种看待这个问题的方法是，对于给定的 f(x)，计算 f( 周围的值时存在一定的误差x）。当你做的时候h更小，f(x+h)–f(x) 变小，但误差保持不变。因此，误差相对于h.