Python数值微分和h的最小值答案

【问题标题】：Python Numerical Differentiation and the minimum value for hPython数值微分和h的最小值
【发布时间】：2019-03-16 18:33:37
【问题描述】：

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

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 之间的相对误差变得很大。

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

【问题讨论】：

【解决方案1】：

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，则精度会提高。这可以在以下到精确导数的距离的对数图中看到。

【讨论】：

【解决方案2】：

您可能会遇到一些数值不稳定，对于 x = 10 和 h =~ 1E-13，np.exp 的参数无论是加还是减 h 都非常接近 10，因此值的近似误差很小np.exp 被非常小的 2 * h 除法显着缩放。

【讨论】：

【解决方案3】：

除了@LutzL 的回答之外，我将添加一些来自第 5.7 章关于数值导数的好书Numerical Recipes 3rd Edition: The Art of Scientific Computing 的信息，尤其是关于选择最优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，在像您这样的小值的情况下，选择可能有点困难。不过，与@LutzL 发布的实际结果非常接近（如果有人可以这么说，请记住浮点运算......）。
您可以使用其他派生公式，除了您正在使用的symmetric。您可能想要使用forward 或backward 评估（如果该函数的评估成本很高，并且您已经预先计算了f(x)。如果您的函数评估成本低，您可能希望使用更高阶多次评估它使精度误差更小的方法（请参阅您问题的评论中提供的five-point stencil on wikipedia）。

【讨论】：

【解决方案4】：

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

【讨论】：