浮点算术错误

Question

我正在使用以下函数来近似函数的导数:

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的fx和x为3.0.在哪里fx:

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

Answer 1

浮点运算(以及整数运算和定点运算)具有一定的粒度:值只能改变一定的步长.对于IEEE-754 64位二进制格式,该步长约为该值的2 ^-52倍(约2.22•10 ^-16).这对于物理测量来说非常小.

但是,当你使h非常小时,f(x)和f(x + h)之间的差异与步长相比并不是很大.差异只能是步长的整数倍.

当导数为d时,f(x)的变化约为h · d.即使您以浮点格式计算f(x)和f(x + h),它们的差值的测量值也必须是步长s的倍数,因此它必须是圆的(h • d/小号)• 小号,其中round(Ý)是Ý取整到最接近的整数.显然,当你使h变小时,h • d/s 较小,因此将其舍入为整数的效果相对较大.

另一种观察方式是,对于给定的f(x),在计算f(x)周围的值时存在一定量的误差.当你使h变小时,f(x + h)-f(x)变小,但误差保持不变.因此,误差相对于h增加.

Answer 2

当您减去两个几乎相同的数字时,结果的精度远低于任何一个输入.这降低了整体结果的精度.

假设您有以下两个数字,小数点后15位:

  1.000000000000001
- 1.000000000000000
= 0.000000000000001

看看发生了什么？结果只有一个好数字.