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)的值，并且期望值和{}之间的相对误差变得很大。在

为什么？这是怎么回事？在

Answer 1

你可能遇到了一些数值上的不稳定性，比如x=10和h=~1E-13np.exp.公司无论h是加上还是减去都非常接近10，所以数值的近似误差很小np.exp.公司以非常小的2*h划分显著的比例

Answer 2

f(x)的求值至多有一个|f(x)|*mu的舍入误差，其中mu是浮点类型的机器常数。因此，中心差分公式的总误差是近似的

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

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

^{pr2}$

它在h = (3*mu)^(1/3)处有一个最小值，对于带有mu=1e-16的双格式，它在h=1e-5左右。在

如果在分母中使用求值点之间的实际差2*h，而不是{}，则精度会提高。这可以在下面的loglog图中看到，它与精确导数的距离。在

Answer 3

除了@LutzL的答案之外，我还将从第5.7章的一本伟大的书Numerical Recipes 3rd Edition: The Art of Scientific Computing中添加一些关于数值导数的信息，特别是关于给定的x的最佳h值的选择：

• 始终选择h，以便h和{}之间存在一个完全可表示的数字。像1/3这样有趣的东西应该被避免，除非x等于{}的东西。在
• 舍入误差大约是epsilon * |f(x) * h|，其中epsilon是浮点精度，Python用双精度表示浮点数，所以它是1e-16。对于更复杂的函数（精度误差会进一步增加），它可能会有所不同，尽管这不是您的情况。在
• 最佳h的选择：对于简单的转发情况，不需要详细说明，除非您的x接近于零（您将在书中找到更多信息），这就是您的情况。您可能希望使用更高的x值。在这种情况下，已经提供了补充答案。在f(x+h) - f(x-h)的例子中，它相当于epsilon ** 1/3 * x，因此大约5e-6乘以{}，在像您这样的小值的情况下，哪种选择可能有点困难。不过，与@LutzL发布的实际结果非常接近（如果有人可以这么说的话，请记住浮点运算…）。在
• 您可以使用其他导数公式，但您使用的symmetric公式除外。您可能需要使用forward或backward求值（如果该函数的求值代价很高，并且您已经预先计算了f(x)。如果函数的求值成本较低，则可能需要使用高阶方法多次求值，以减小精度误差（请参见问题注释中提供的five-point stencil on wikipedia）。在