2024-04-18 04:53:39 发布
网友
我用以下代码计算一阶导数:
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)的值,并且期望值和{}之间的相对误差变得很大。在
x == 10
h
10E-14
Df
f(10)
为什么?这是怎么回事?在
你可能遇到了一些数值上的不稳定性,比如x=10和h=~1E-13np.exp.公司无论h是加上还是减去都非常接近10,所以数值的近似误差很小np.exp.公司以非常小的2*h划分显著的比例
f(x)的求值至多有一个|f(x)|*mu的舍入误差,其中mu是浮点类型的机器常数。因此,中心差分公式的总误差是近似的
f(x)
|f(x)|*mu
mu
2*|f(x)|*mu/(2*h) + |f'''(x)|/6 * h^2
在本例中,指数函数等于其所有导数,因此误差与
它在h = (3*mu)^(1/3)处有一个最小值,对于带有mu=1e-16的双格式,它在h=1e-5左右。在
h = (3*mu)^(1/3)
mu=1e-16
h=1e-5
如果在分母中使用求值点之间的实际差2*h,而不是{},则精度会提高。这可以在下面的loglog图中看到,它与精确导数的距离。在
2*h
除了@LutzL的答案之外,我还将从第5.7章的一本伟大的书Numerical Recipes 3rd Edition: The Art of Scientific Computing中添加一些关于数值导数的信息,特别是关于给定的x的最佳h值的选择:
x
1/3
epsilon * |f(x) * h|
1e-16
f(x+h) - f(x-h)
epsilon ** 1/3 * x
5e-6
symmetric
forward
backward
你可能遇到了一些数值上的不稳定性,比如x=10和h=~1E-13np.exp.公司无论h是加上还是减去都非常接近10,所以数值的近似误差很小np.exp.公司以非常小的2*h划分显著的比例
f(x)
的求值至多有一个|f(x)|*mu
的舍入误差,其中mu
是浮点类型的机器常数。因此,中心差分公式的总误差是近似的在本例中,指数函数等于其所有导数,因此误差与
^{pr2}$它在
h = (3*mu)^(1/3)
处有一个最小值,对于带有mu=1e-16
的双格式,它在h=1e-5
左右。在如果在分母中使用求值点之间的实际差},则精度会提高。这可以在下面的loglog图中看到,它与精确导数的距离。在
2*h
,而不是{除了@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
乘以{symmetric
公式除外。您可能需要使用forward
或backward
求值(如果该函数的求值代价很高,并且您已经预先计算了f(x)
。如果函数的求值成本较低,则可能需要使用高阶方法多次求值,以减小精度误差(请参见问题注释中提供的five-point stencil on wikipedia)。在相关问题 更多 >
编程相关推荐