我正在使用以下函数来近似函数在某一点的导数:
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)
之前,我一直在期待相反的情况。
最佳答案
浮点运算(以及整数运算和定点运算)具有一定的粒度:值只能按一定的步长变化。对于 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的倍数,所以一定是round(h•d/s)•s,其中 round(y) 是 y 四舍五入到最接近的整数。显然,随着 h 变小,h•d/s 也变小,因此四舍五入的效果到一个整数相对较大。
另一种看待这个问题的方式是,对于给定的 f(x),在计算 f(x) 附近的值时存在一定量的误差。随着 h 变小,f(x+h)–f(x) 变小,但是错误保持不变。因此,误差相对于 h 增加。
关于python - 浮点运算错误,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/18995148/