假设我要在x = pi / 2处以数值计算f(x)= cos(x)的一阶导数。df / dx = -sin(x),因此为df / dx = -1。为此,我使用最简单的公式:
df / dx =(f(x + h)-f(x))/ h + O(h)。
O(h)是与h成正比的误差,因此,数学上O(h)随h变为零而变为零。我了解在计算机中情况有所不同。
我是否应该期望如果我使用的是双精度,这应该给我15到17个十进制有效数字,我应该能够将精确结果df / dx = -1.0乘以〜10 ^ -15?即我应该找到(df / dx)_numerical + 1.0〜10 ^ -15?
这是我发现的h的不同值的结果:
h值:
[0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10]
(df / dx)_数字:
[-0.9999999983332231,
-0.9999999999898844,
-0.9999999999175667,
-1.0000000005838656,
-0.999999993922529,
-1.000000082740371,
-1.000000082740371]
这是预期的吗?为什么?为什么对于h = 10 ^ -5可获得最佳结果?
使用有限差分计算导数很容易失去意义。问题是要减去的上下文,而不是要进行的除法。基本上,cos(x)-cos(x+delta)
当cos(x)
和cos(x+delta)
都接近零时,就可以正常工作...但是随着它们远离零(但仍然彼此接近),结果的精度将急剧下降。在某一点上,由于使用较小的增量而导致的精度提高将被重要性降低导致的精度下降所抵消。对您来说,这似乎发生在1e-5左右,但这不是根本(在该区域中,错误往往会反弹很多)。
关于如何进行数值稳定的有限差分的文章很多,但是第一法则是“不要对差分的微小性过于贪婪”。有关更深入(而不模糊)的信息,我可以推荐Forman Acton的(通常)有效的数值方法。
抱歉,我忘了提到最重要的部分。分子中的重要性丧失很重要,但是分子中的重要性丧失x+h
更为重要。而且该部分实际上很容易修复。
随着x
变大,近似值x+h
变差(再次失去重要性)。基本上,您用于评估的步长与分母中的步长不匹配。但!由于您实际上并不关心的确切值h
,因此只需找出四舍五入h
后最终使用的值x+h
,然后在分母中使用该值即可。从本质上讲,你算算x2=x+h
,然后h'=x2-x
。(Sterbenz定理)的计算h'
是精确x >> h
的,从而消除了特定的重要性损失。
示例代码:
import math
def calcDerivAt_orig(f, x, h):
x1 = x
x2 = x+h
y1 = f(x1)
y2 = f(x2)
return (y2-y1)/h
def calcDerivAt_fixed(f, x, h):
x1 = x
x2 = x+h
y1 = f(x1)
y2 = f(x2)
return (y2-y1)/(x2-x1)
for h in [0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10]:
dOrig = calcDerivAt_orig(math.cos, math.pi/2, h)
origErr = abs(-1 - dOrig)
dFixed = calcDerivAt_fixed(math.cos, math.pi/2, h)
fixedErr = abs(-1 - dFixed)
print("h = {}: origErr = {}, fixedErr = {}".format(h, origErr, fixedErr))
产生:
h = 0.0001: origErr = 1.66677693869e-09, fixedErr = 1.66666680457e-09
h = 1e-05: origErr = 1.01155750443e-11, fixedErr = 1.6666779068e-11
h = 1e-06: origErr = 8.24332824223e-11, fixedErr = 1.66644475996e-13
h = 1e-07: origErr = 5.83865622517e-10, fixedErr = 1.55431223448e-15
h = 1e-08: origErr = 6.07747097092e-09, fixedErr = 0.0
h = 1e-09: origErr = 8.27403709991e-08, fixedErr = 0.0
h = 1e-10: origErr = 8.27403709991e-08, fixedErr = 0.0
一点也不差。当然,我们在选择时作弊x=pi/2
,因为cos(x)
那里的值接近于零;在x=1.5
某程度上,由于h
我最初描述的重要性下降,您仍然会看到错误随着缩小而迅速膨胀:
h = 0.0001: origErr = 3.53519750529e-06, fixedErr = 3.5351976152e-06
h = 1e-05: origErr = 3.53675653653e-07, fixedErr = 3.53669118991e-07
h = 1e-06: origErr = 3.52831192041e-08, fixedErr = 3.53651797846e-08
h = 1e-07: origErr = 4.03034106089e-09, fixedErr = 3.44793649187e-09
h = 1e-08: origErr = 5.5453325265e-09, fixedErr = 5.1691428915e-10
h = 1e-09: origErr = 7.21702790862e-08, fixedErr = 1.03628251535e-08
h = 1e-10: origErr = 1.6659127966e-08, fixedErr = 6.58739718329e-08
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句