Python 数值微分和 h 的最小值

Question

我使用以下代码计算一阶导数：

def f(x):
   f = np.exp(x)
   return f

def dfdx(x):
   Df = (f(x+h)-f(x-h)) / (2*h)
   return Df

例如，对于 x == 10，这很好用。但是当我将 h 设置为 10E-14 左右或以下时， Df 开始得到真正远离预期值 f(10) 的值，并且预期值和 Df 之间的相对误差变得很大。

这是为什么？这里发生了什么？

Answer 1

你可能遇到了一些数值不稳定的问题，对于x = 10和h =~ 1E-13，np.exp的自变量非常接近10，无论是加还是减h，所以近似误差很小np.exp 的值被非常小的 2 * h.

的除法显着缩放

Answer 2

这 Python tutorial 解释了精度有限的原因。综上所述，小数最终是用二进制表示的，精度在17位有效数字左右。所以，你是对的，它在 10E-14 之后变得模糊。

Answer 3

f(x) 的计算最多有 |f(x)|*mu 的舍入误差，其中 mu 是浮点类型的机器常量。因此中心差分公式的总误差约为

2*|f(x)|*mu/(2*h)  +  |f'''(x)|/6 * h^2

在本例中，指数函数等于它的所有导数，因此误差与

成正比

mu/h + h^2/6

的最小值在 h = (3*mu)^(1/3)，对于 mu=1e-16 的双精度格式约为 h=1e-5。

如果在分母中使用评估点之间的实际差异 (x+h)-(x-h) 而不是 2*h，则精度会提高。这可以在下面的精确导数距离的对数对数图中看出。

Answer 4

除了 I will add some info from a great book Numerical Recipes 3rd Edition: The Art of Scientific Computing第5.7章关于数值导数的回答，特别是最佳h值的选择对于给定的 x:

始终选择 h，以便 h 和 x 相差一个可精确表示的数字。应该避免像 1/3 这样有趣的东西，除非 x 等于 14.3333333.
舍入误差约为epsilon * |f(x) * h|，其中epsilon是浮点精度，Python表示双精度浮点数，所以是1e-16。对于更复杂的函数（精度误差会进一步出现），它可能会有所不同，尽管这不是你的情况。
最佳 h 的选择：没有详细说明，对于简单的前向情况，它将是 sqrt(epsilon) * x，除非您的 x 接近于零（您将在书），这是你的情况。在这种情况下，您可能希望使用更高的 x 值，已提供补充答案。在您的示例中 f(x+h) - f(x-h) 的情况下，它将达到 epsilon ** 1/3 * x，因此大约 5e-6 乘以 x，对于像这样的小值，选择可能有点困难你的。不过，与发布的实际结果非常接近（如果可以这么说，请记住浮点运算...）。
您可以使用其他导数公式，但您正在使用的 symmetric 除外。您可能想要使用 forward 或 backward 评估（如果函数的评估成本很高并且您事先计算了 f(x)。如果您的函数评估成本很低，您可能需要评估它多次使用高阶方法使精度误差更小（请参阅问题评论中提供的five-point stencil on wikipedia）。

Python Numerical Differentiation and the minimum value for h