在C/C++中实现派生

Question

f(x)通常以编程方式计算的导数如何确保最大精度？

我正在实现Newton-Raphson方法,它需要获取函数的导数.

Answer 1

我同意@erikkallen这(f(x + h) - f(x - h)) / 2 * h是数值近似衍生物的常用方法.但是,获得正确的步长h有点微妙.

近似误差(f(x + h) - f(x - h)) / 2 * h随着h变小而减小,表示你应该h尽可能小.但随着h变小,浮点减法的误差增加,因为分子需要减去几乎相等的数字.如果h太小,你可以松一个减法中的精度很高.所以在实践中你必须选择一个不太小的值,h以最小化近似误差和数值误差的组合.

根据经验,您可以尝试h = SQRT(DBL_EPSILON)哪里DBL_EPSILON是最小的双精度数e,以便1 + e != 1在机器精度方面. DBL_EPSILON是10^-15这样你可以使用h = 10^-7或10^-8.

有关更多详细信息,请参阅有关选择微分方程的步长的这些说明.

Answer 2

Newton_Raphson假设你可以有两个函数f(x)及其导数f'(x).如果您没有作为函数可用的导数并且必须估计原始函数的导数,那么您应该使用另一个根查找算法.

维基百科根发现提供了几个建议,就像任何数值分析文本一样.

Answer 3

1)第一种情况:

- 相对舍入误差,对于double为2 ^ { - 16},对于float为2 ^ { - 7}.

我们可以计算总误差:

假设您正在使用双浮动操作.因此,h的最佳值是2sqrt(DBL_EPSILON/f''(x)).你不知道f''(x).但你必须估计这个值.例如,如果f''(x)约为1,则h的最佳值为2 ^ { - 7},但如果f''(x)约为10 ^ 6,则h的最佳值为2 ^ { - 10}!

2)第二种情况:

注意,第二近似误差比第一近似误差快0.但如果f'''(x)非常滞后,那么第一个选项更为可取:

注意,在第一种情况下,h与e成比例,但在第二种情况下,h与e ^ {1/3}成比例.对于双浮动操作,e ^ {1/3}是2 ^ { - 5}或2 ^ { - 6}.(我想f'''(x)约为1).

哪种方式更好？如果您不知道f''(x)和f'''(x),或者您无法估计这些值,那么这是未知的.据信第二种选择是优选的.但如果你知道f'''(x)非常大,请使用第一个.

h的最佳值是多少？假设f''(x)和f'''(x)约为1.还假设我们使用双浮动操作.然后在第一种情况下,h约为2 ^ { - 8},在第一种情况下,h约为2 ^ { - 5}.如果您知道f''(x)或f'''(x),请更正此值.

Answer 4

fprime(x) = (f(x+dx) - f(x-dx)) / (2*dx)

对于一些小dx.

Answer 5

你对f(x)有什么了解？如果你只有f作为黑盒子,你唯一能做的就是数字逼近导数.但准确性通常不是那么好.

你可以做很多更好,如果你可以触摸,计算F中的代码.尝试"自动区分".有一些很好的库可用.通过一些库魔法,您可以轻松地将您的功能转换为自动计算衍生物的功能.有关简单的C++示例,请参阅此德语讨论中的源代码.

Answer 6

您肯定要考虑John Cook关于选择h的建议，但是您通常不希望使用居中的差异来近似导数。主要原因是，如果使用前向差异，则需要额外的功能评估，即

f'(x) = (f(x+h) - f(x))/h

然后，您将免费获得f（x）的值，因为您需要已经为牛顿方法计算了它。当您有一个标量方程时，这没什么大不了的，但是如果x是一个向量，则f'（x）是一个矩阵（雅可比行列式），您将需要进行n个额外的函数求值来近似它使用居中差异法。