时间戳

首页 > 解决方案 > C++ 数值导数

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 f(double,void*), double h0){
 return(4.0*derivative1(var,f,h0) - derivative1(var,f,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,f,h0);
  double gtest=gderivative(x,f,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, diff = 1.951474537520426e-09, percent error = 9.199339191957163e-08
h0 = 1.000000000000000e-07
Exact      = 2.121320343559643e+00
Naive step = 2.121320341608168e+00, diff = 1.951474537520426e-09, 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而 for Richardson h0 = 1.000000000000000e-03。这让我很困惑，因为这些都不小。对任何其他既有效又精确（double）的好选择（简单算法/库）有什么建议吗？

标签： c++numericnumerical-methodsderivative

解决方案

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

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

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

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

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

推荐阅读