为什么我的派生函数产生奇怪的结果？

Question

我正在写一个小的微积分库供个人使用.它是标准的微积分工具 - 取一阶导数,n阶导数,黎曼和等.我遇到的一个问题是,n阶导数函数对某些h值(有限差分)返回明显的虚假结果.

代码在这里和下面:

typedef double(*math_func)(double x);
inline double max ( double a, double b ) { return a > b ? a : b; }

//Uses the five-point-stencil algorithm.
double derivative(math_func f,double x){
    double h=max(sqrt(DBL_EPSILON)*x,1e-8);
    return ((-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h));
}

#define NDEPS (value)
double nthDerivative(math_func f, double x, int N){
    if(N<0) return NAN; //bogus value of N
    if(N==0) return f(x);
    if(N==1) return derivative(f,x);
    double* vals=calloc(2*N+9,sizeof(double)); //buffer region around the real values
    if(vals==NULL){ //oops! no memory
        return NAN;
    }
    int i,j;
    //don't take too small a finite difference
    double h=max(sqrt(DBL_EPSILON)*x,NDEPS);
    for(i=-(N+4);i<=N+4;i++){
        vals[i+N+4]=derivative(f,x+h*i);
    }
    //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    for(j=1;j<N;j++){
        double *vals2=calloc(2*N+9,sizeof(double));
        for(i=2;i<2*N+7;i++){
            vals2[i]=(-vals[i+2]+8*vals[i+1]-8*vals[i-1]+vals[i-2])/(12*h);
        }
        free(vals);
        vals=vals2;
        //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    }
    double result=vals[N+4];
    free(vals);
    return result;
}

我给测试这个函数的一个样本问题是当x = pi时sin(x)的五阶导数.众所周知,答案是-1.当我试图改变NDEPS("Nth derivative epsilon")的值时,问题出现了.

• 当NDEPS = 1.5625e-2(1/64.0)时:sin =的五阶导数x = pi:-1.0003e + 00(看起来很好).
• 当NDEPS = 1e-5(1/100000.0)时:sin =的五阶导数x = pi:2.4302e + 11(我在这里称废话).

为什么会这样？它与sin()函数的性质有关吗？还是因为浮点精度问题？

Answer 1

数值微分是一个难题.关键问题是有限差分近似

f'(x) =(approx) (f(x+h)-f(x-h)) / (2*h)

是一个取消的食谱.

这意味着您必须在两个错误之间选择谨慎的权衡:如果h太大,您将在有限差分中产生大的近似误差.如果h太小,则会产生大的(可能是灾难性的)数值误差.维基百科关于数值微分的文章很好地说明了这一点.随着衍生物的顺序增加,这些问题会加剧.

这就是为什么自动微分是如此重要的原因 - 本质上,当所有给出的算法都是计算基函数的算法时,它允许您计算精确的导数.如果函数足够简单,你可以象征性地确定导数,那么你应该这样做.

如果你确实需要进行数值微分,那么你可以应用很多数字技巧.Numerical Recipes有一个很好的处理 - 看看第5.7节,从第186页开始.