为什么不同的步长在逼近导数时会产生很大的数值误差？

编程问答 2022-05-31

问题描述

我正在尝试一些代码来帮助我逼近任意函数的导数。我看到另一个 post 上有四个选项：

  1. 有限差分
  2. 自动导数
  3. 符号微分
  4. 手动计算导数

我发现我的方法在第一个选项中效果最好，该选项带有“容易出现数值错误”的注释。所以我知道这种方法并不准确，这很好。

话虽如此，我对不同数据类型可以存储的大小数字进行了一些研究，并在此 post 中发现它可以非常小（大约 10–308 ) 和“在正常范围内，基本操作的结果将在格式的正常精度内准确”。

话虽如此，对于以下我探索不同大小间隔的代码片段，我似乎得到了非常糟糕的结果；最小的差异不应小于 10–27（10–9，立方），远大于极限值。我希望得到更具体的答复？

epsilon = 0.01 # is "small" w.r.t. to 3

def approx_derivative(func): # rough derivative factory function
  return lambda x : (func(x + epsilon) - func(x)) / epsilon

while epsilon > 10**-9:
  nth_deriv = lambda x : x ** 3 # 0th derivative
  for i in range(5): # should read about 27,27,18,6,0
    print(nth_deriv(3),end=',')
    nth_deriv = approx_derivative(nth_deriv) # take derivative
  print('\n')
  epsilon *= 0.1

输出为：

27,27.090099999999495,18.0599999999842,6.000000002615025,-3.552713678800501e-07,27.009000999996147,18.00600000123609,6.000000496442226,-0.007105427357601002,27.00090001006572,18.000599766310188,6.004086117172847,-71.05427357601002,27.000090000228735,18.000072543600254,3.5527136788005005,355271.36788005003,27.000009005462285,17.998047496803334,0.0,3552713678.8005,27.000000848431675,18.11883976188255,-35527136788004.99,27.0000001023618,3.552713678800497e+17,27.000002233990003,

正如我们在前几个示例中看到的，结果并不准确，但相当不错。但是，对于某些间隔大小，某些值会被放大；其他归0；有些是完全错误的，比如给出一半的值，尽管直觉认为它们应该对较小的 epsilon 变得更准确。我可以将哪些主要因素归因于此错误？我应该注意/注意什么？是否有我应该担心的错误（例如除以 0）？

epsilon 是否有一个通常被认为是使用浮点数进行计算的“最佳”值？或者是否有根据您的输入选择合适大小的 epsilon 的“经验法则”？是否有比我实现的导数更适合使用的导数定义？

解决方法

首先，浮点类型中可表示的最小值在这里无关紧要。精度是一个问题。 Python 规范没有具体说明 Python 实现使用什么浮点格式，但许多使用 IEEE-754 binary64 格式，其有效位有 53 位，这意味着值的最小可表示变化约为 2−52 乘以表示的数字的大小。我们将研究这会如何影响您的计算。

计算导数 f' 的近似值，计算 (x+e)3−x3.结果是 3ex2 + 3e2x + e3，大约 3ex2，然后除以 e 得到大约 3 x2。请注意，涉及到在幅度上 x3 左右的两个项，它们相差约 3ex2。由于您在 x = 3 处取导数，因此 3ex2 恰好等于 ex 3，因此它比 x3 低 e 倍。只要 e 是适度的，x3 附近的两个 binary64 数字很容易相差 ex3 不失准确性。

计算二阶导数 f'' 的近似值，计算 ((x+2e)3−( x+e)3)−((x+e)3−x3) = (x+2e)3 −2(x+e)3+x3 = 6e2x + 6e3。在此过程中，它被 e2 划分，但我们现在可以忽略它以进行错误分析。原始项的量级约为 x3，但我们的减法结果约为 6e2x。因此，这取决于能够相差大约 e2 的数字。由于 2−52 大约是 2•10−16，当 e 大约是 10−8 时，这会相当失败，在此之前错误较少。当 e 大约是这个大小时，binary64 格式根本无法准确地记录数字之间的差异，以在以后减去它们时保持准确性。

对于三阶导数，您需要表示 e3 附近相对于主要星等的差异，因此它在 e = 附近相当失败10−5，从那时起e3就接近2•10−16的极限。 >

derivativeieee-754numericnumericprecisionpython