如何进行 f(x) 的 n 阶导数的数值计算？

Question

如何进行 f(x) 的 n 阶导数的数值计算？

4

我实现了一个C++代码，用于数值求解在点x_0处的函数的n阶导数:

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阶导数在n大于3时会发散。你知道如何解决吗？

- Gianluca Bianco

2

也许您可以实现其中一个？https://math.stackexchange.com/q/130192 - konsolas

请查看 https://en.wikipedia.org/wiki/Numerical_differentiation 了解技术。 - Ahmed Masud

7

注意：1/3 == 0 的意思是“1除以3等于0”。 - Ranoiaetep

1

pow(__DBL_EPSILON__, 1/3); 不会按预期工作，而且每次函数调用都计算会非常昂贵，最好将其存储为静态常量。 - phuclv

1

@GianlucaBianco 我的意思是这样。将一个常量存储在静态变量中，或者更好的是，在编译时计算一次，而不是每次函数调用都计算。 - phuclv

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- chux - Reinstate Monica · Accepted Answer

这不是一个好的实现。

至少存在以下问题：

整数运算

使用FP（浮点）计算，因为1/3会变成零。

1/3 --> 1.0/3

在n==1的情况下使用立方根最优

但并不一定适用于其他n。 @Eugene

错误的ε值

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

OP（原始发布者）的+/-某些ε值对于固定点来说是好的，但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)，但调用却为 f(x_0)

小问题：名称混淆

为什么用 first_term 表示 x_2，而用 second_term 表示 x_1？