我同意@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。
有关更多详细信息，请参阅notes中关于为微分方程选择步长的内容。
此外，如果你想要非常非常准确的结果，你可以使用特征库。首先尝试下面的代码，不做任何修改，然后用Eigen include取消第一行的注解：

// #include "Eigen/Dense" // uncomment this, and check the result

#include <functional>
#include <iostream>
#include <limits>
#include <cmath>

typedef std::function<double(double)> RealFunc;
typedef std::function<double(std::function<double(double)>, double)> RealFuncDerivative;

double FibonacciFunc(double x) {
    return pow(x, 3) + 2.0 * pow(x, 2) + 10.0 * x - 20.0;
}

double derivative(RealFunc f, double x) {
    double h = sqrt(std::numeric_limits<double>::epsilon());
    return (f(x + h) - f(x - h)) / (2.0 * h);
}

double NewtonsMethod(RealFunc f, RealFuncDerivative d, double x0, double precision) {
    double x = x0;
    for (size_t i = 0;; i++) {
        x = x - (f(x) / d(f, x));

        if (abs(f(x)) < precision) {
            return x;
        }
    }
}

int main() {
    RealFunc f{FibonacciFunc};
    RealFuncDerivative d{derivative};

    std::cout << NewtonsMethod(f, d, 1.0, 10e-4) << "\n";
}

第一次运行的结果应该是1.41176，使用Eigen库的结果应该是1.36881（均使用2次迭代计算）。第二个结果与使用5位小数的分析真实结果完全相同。
您观察到的结果差异是由于Eigen影响代码中的浮点精度设置和舍入模式，以优化其自身操作的数值稳定性和性能。这可能会导致更好的结果，特别是当你有敏感的数值算法，如牛顿-拉夫森方法。
请注意，使用Eigen，即使使用大的h（如h=10而不是h=SQRT(DBL_EPSILON)），算法仍然收敛得很好，尽管相当慢（1.36877在47次迭代下）。所以总而言之，这不是必要的，但仍然首选使用尽可能小的h。

Newton_Raphson假设你可以有两个函数f（x）和它的导数f '（x）。如果你没有可用的导数作为一个函数，必须估计从原始函数的导数，那么你应该使用另一个求根算法。
Wikipedia root finding给出了一些建议，就像任何数值分析文本一样。

1)第一例：

-相对舍入误差，双精度约为2^{-16}，浮点数约为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），请更正此值。

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

对于一些小DX。

关于f（x）你知道些什么？如果你只有f作为一个黑盒子，你唯一能做的就是用数字近似导数。但准确性通常不是那么好。
如果你能接触到计算f的代码，你可以做得更好。试试"automatic differentiation"。有一些很好的图书馆可以使用。使用一点库魔法，您可以轻松地将函数转换为自动计算导数的函数。有关简单的C++示例，请参阅source code in this德语讨论。

你肯定想考虑John Cook的建议来选择h，但你通常不想使用中心差来近似导数。主要原因是它需要额外的函数计算，如果你使用一个向前的差异，即，

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

然后你就可以免费得到f（x）的值，因为你需要在牛顿法中计算它。当你有一个标量方程时，这不是什么大问题，但是如果x是一个向量，那么f '（x）是一个矩阵（雅可比矩阵），你需要做n个额外的函数计算来使用中心差分方法近似它。

除了约翰·D。Cooks在上面回答说，重要的是不仅要考虑浮点精度，还要考虑函数f（x）的鲁棒性。例如，在金融中，常见的情况是f（x）实际上是蒙特卡罗模拟，并且f（x）的值有一些噪声。在这些情况下，使用非常小的步长可能会严重降低导数的精度。

通常，信号噪声对导数质量的影响比其他任何东西都大。如果f（x）中确实存在噪声，Savtizky-Golay是一种优秀的平滑算法，通常用于计算良好的导数。简而言之，SG将多项式局部拟合到您的数据，然后此多项式可用于计算导数。
保罗