C++ 数值导数

如何解决C++ 数值导数

这个问题在不同平台上出现了一千次。但是，我仍然需要了解一些东西。 这是一个完整的例子：

#include <iostream>
#include <iomanip>
#include <cmath>

// function to derivative
double f (double var,void * params){
  (void)(params); 
  return pow (var,1.5);
}

// Naive central step derivative
double derivative1(double var,double f(double,void*),double h){ 
return (f(var+h,NULL) - f(var-h,NULL) )/2.0/h; 
}

// Richardson 5-point rule
double gderivative(double var,double h0){
 return(4.0*derivative1(var,f,h0) - derivative1(var,2.0*h0))/3.0;
}


int main (void){


  for (int i=10;i>0;i--){
  double h0=pow(10,-i);
  double x=2.0;
  double exact = 1.5 * sqrt(x);
  double test1=derivative1(x,h0);
  double gtest=gderivative(x,h0);
  std::cout << "h0 = " << h0 << std::endl;
  std::cout << "Exact      = "  << std::scientific<<std::setprecision(15) << exact << std::endl;
  std::cout << "Naive step = "  << std::setprecision(15) << test1 <<",diff = " << std::setprecision(15)<< exact-test1 << ",percent error = " << (exact-test1)/exact*100.0 << std::endl;
  std::cout << "Richardson = "  << std::setprecision(15) << gtest <<",diff = " << std::setprecision(15)<< exact-gtest << ",percent error = " << (exact-gtest)/exact*100.0 << std::endl;



}


 return 0;
}

输出是

h0 = 1e-10
Exact      = 2.121320343559643e+00
Naive step = 2.121318676273631e+00,diff = 1.667286011475255e-06,percent error = 7.859661632610832e-05
Richardson = 2.121318306199290e+00,diff = 2.037360352868944e-06,percent error = 9.604208808228318e-05
h0 = 1.000000000000000e-09
Exact      = 2.121320343559643e+00
Naive step = 2.121320674675076e+00,diff = -3.311154328500265e-07,percent error = -1.560893119491818e-05
Richardson = 2.121320748689944e+00,diff = -4.051303013064000e-07,percent error = -1.909802555452698e-05
h0 = 1.000000000000000e-08
Exact      = 2.121320343559643e+00
Naive step = 2.121320341608168e+00,diff = 1.951474537520426e-09,percent error = 9.199339191957163e-08
Richardson = 2.121320341608168e+00,percent error = 9.199339191957163e-08
h0 = 1.000000000000000e-07
Exact      = 2.121320343559643e+00
Naive step = 2.121320341608168e+00,percent error = 9.199339191957163e-08
Richardson = 2.121320340868019e+00,diff = 2.691623368633600e-09,percent error = 1.268843424240664e-07
h0 = 1.000000000000000e-06
Exact      = 2.121320343559643e+00
Naive step = 2.121320343606570e+00,diff = -4.692690680485612e-11,percent error = -2.212155601454860e-09
Richardson = 2.121320343643577e+00,diff = -8.393419292929138e-11,percent error = -3.956695799581460e-09
h0 = 1.000000000000000e-05
Exact      = 2.121320343559643e+00
Naive step = 2.121320343584365e+00,diff = -2.472244631235299e-11,percent error = -1.165427295665677e-09
Richardson = 2.121320343595468e+00,diff = -3.582467655860455e-11,percent error = -1.688791448560268e-09
h0 = 1.000000000000000e-04
Exact      = 2.121320343559643e+00
Naive step = 2.121320343340116e+00,diff = 2.195266191051815e-10,percent error = 1.034858406801534e-08
Richardson = 2.121320343561791e+00,diff = -2.148059508044753e-12,percent error = -1.012604963020456e-10
h0 = 1.000000000000000e-03
Exact      = 2.121320343559643e+00
Naive step = 2.121320321462283e+00,diff = 2.209735949776359e-08,percent error = 1.041679516479040e-06
Richardson = 2.121320343559311e+00,diff = 3.317346397579968e-13,percent error = 1.563812088849040e-11
h0 = 1.000000000000000e-02
Exact      = 2.121320343559643e+00
Naive step = 2.121318133840577e+00,diff = 2.209719065504601e-06,percent error = 1.041671557157002e-04
Richardson = 2.121320343601055e+00,diff = -4.141265108614789e-11,percent error = -1.952211093995174e-09
h0 = 1.000000000000000e-01
Exact      = 2.121320343559643e+00
Naive step = 2.121099269013200e+00,diff = 2.210745464426012e-04,percent error = 1.042155406248691e-02
Richardson = 2.121320759832334e+00,diff = -4.162726914280768e-07,percent error = -1.962328286210455e-05

我相信标准 GSL gsl_deriv_central 采用理查森程序。

现在选择 h0 的常见说法是理论上应该选择尽可能小以提高导数的精度，但在数值上它不应该太小以至于我们达到浮点四舍五入从而破坏精度。人们经常说最佳选择应该在 1e-6 - 1e-8 左右。我的问题是：

泛型导数中 h0 的最佳选择是什么？

我应该逐案检查吗？通常可能无法获得准确的结果进行检查。这种情况下应该怎么办？ 现在在这种特殊情况下，Naive step 的最佳选择似乎是 h0 = 1.000000000000000e-05 而 Richardson h0 = 1.000000000000000e-03。这让我很困惑，因为它们并不小。 关于任何其他有效且精确（double）的好的选择（简单算法/库）的任何建议？

解决方法

这完全符合预期。中心差对精确导数具有二阶误差 O(h^2)。函数评估有一个大小为 mu 的误差，机器常数（对于适度缩放的测试示例）。因此，中心差的评估误差为 mu/h 量级。如果这两个影响大致相等，则总体误差最小，因此 h=mu^(1/3) 给出 h=1e-5，误差约为 1e-10。

理查森外推的相同计算给出了精确导数的误差顺序 O(h^4)，导致 h=mu^(1/5)=1e-3 作为最佳步长，误差约为 1e-12。

两种方法在较大步长样本上的误差的对数图

在实际应用中，您可能希望缩放 h，以便指示的大小与 x 的大小相关。确切的最优值还取决于函数导数的大小，如果它们增长得太厉害。

要获得更精确的导数值，您需要更高或多精度的数据类型，或采用自动/算法微分，其中可以以与函数值相同的精度计算一阶和二阶导数。

,

个人意见警告

根据我的经验，我发现使用与它影响的变量相比较小的步长会更好。

例如，我通常使用这样的东西：

auto dfdx(std::function<double (double)> f,double x,double h0=1e-5){
    h0=std::abs(x)*h0+h0; // this is 1e-5 if x<<1e-5,and |x|*1e-5 otherwise
    
    return ( f(x+h0) - f(x-h0) )/2/h0;
}

这应该有效，因为有限差分是由泰勒展开式激发的。也就是说，只要 x