f(x)的n次导数的数值实现？

Question

我实现了一个 C++ 代码来数值求解函数在一个点中的 n 阶导数 x_0:

double n_derivative( double ( *f )( double ), double x_0, int n )
 {
  if( n == 0 ) return f( x_0 );
  else 
   {
    const double h = pow( __DBL_EPSILON__, 1/3 );
    double x_1 = x_0 - h;
    double x_2 = x_0 + h;
    double first_term = n_derivative( f, x_2, n - 1 );
    double second_term = n_derivative( f, x_1, n - 1);

    return ( first_term - second_term ) / ( 2*h );
   }
 }

我想知道这对您来说是否是一个很好的实现，或者是否可以用 C++ 更好地编写它。问题是我注意到当 n 的值大于 3 时，n 阶导数发散。你知道如何解决这个问题吗？

Answer 1

这不是一个好的实现

至少这些问题。

整数数学

使用 FP 数学，因为 1/3 为零。

1/3 --> 1.0/3

使用最适合n==1

的立方根

但不一定是其他n。

错误的 epsilon

以下代码仅对 |x_0| 约 1.0 有用。当x_0很大时，x_0 - h可能等于x_0。当x_0较小时，x_0 - h可能等于-h.

OP 的 +/- 一些 epsilon 适用于固定点，但 double 是浮动点。

// Bad
const double h = pow( __DBL_EPSILON__, 1.0/3 );
double x_1 = x_0 - h;

需要相对缩放。

#define EPS cbrt(DBL_EPSILON) // TBD code to well select this 
if (fabs(x_0) >= DBL_MIN && isfinite(x_0)) {
  double x_1 = x_0*(1.0 - EP3);
  double x_2 = x_0*(1.0 + EPS);
  double h2 = x_2 - x_1;
  ...
} else {
  TBD_Code for special cases
}

无效代码

f是double ( *f )( int, double )，但是call是f( x_0 )

轻微：混淆名称

为什么 first_term 与 x_2 和 second_term 与 x_1？

f(x)的n次导数的数值实现？

Numerical implementation of n-th derivative of f(x)?

c

c++

math

numerical-integration

derivative