c - 为什么我的导数函数会产生奇怪的结果？

Question

我正在编写一个供个人使用的小型微积分库。其中包含标准微积分工具 - 采用一阶导数、n 阶导数、黎曼和等。我面临的一个问题是，对于 h 的某些值（有限差分），n 阶导数函数返回明显虚假的结果。

代码在这里和下面：

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) 的 5 阶导数。众所周知，答案是-1。当我尝试改变 NDEPS（“N 次导数 epsilon”）的值时，问题就来了。

• 当 NDEPS=1.5625e-2 (1/64.0) 时：x=pi:-1.0003e+00 处 sin() 的 5 阶导数（虽然看起来不错）。
• 当 NDEPS=1e-5 (1/100000.0) 时：x=pi:2.4302e+11 处 sin() 的 5 阶导数（我在这里说废话）。

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

Answer 1

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

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

是取消的秘诀。

这意味着您必须在两个误差之间做出谨慎的权衡：如果 h 太大，您将在有限差分中产生很大的近似误差。如果 h 太小，您将产生很大的（可能是灾难性的）数值错误。这在关于数值微分的维基百科文章中得到了很好的说明。随着导数的阶数增加，这些问题更加严重。

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

如果您确实需要进行数值微分，则可以应用很多数值技巧。Numerical Recipes有一个很好的处理方法——看看从第 186 页开始的第 5.7 节。

Answer 2

这是您的代码的改编版。除了使用更多的空白之外，主要的修改是使函数NDEPS成为一个参数nthDerivative()，以便可以用不同的值调用它并添加丰富的打印。我还必须在普通derivative()功能上发挥创造力；代码编译正确（但我真的不想用 assertion 做出哲学或神学陈述assert(sin == fun);，但这确实意味着代码在没有警告的情况下编译，并且它认识到这个衍生函数的局限性）。

代码

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define max(a, b)    (((a) > (b)) ? (a) : (b))

#define PRIe_double    "%21.15e"

typedef double(*math_func)(double x);

static double derivative(math_func fun, double x) { assert(sin == fun); return cos(x); }

static double nthDerivative(math_func f, double x, int N, double NDEPS)
{
    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);
    printf("h = " PRIe_double "\n", h);

    for (i = -(N+4); i <= N+4; i++)
    {
        vals[i+N+4] = derivative(f, x+h*i);
        printf("%2d: deriv(" PRIe_double ") = " PRIe_double "\n", i, x+h*i, vals[i+N+4]);
    }

    for (j = 1; j < N; j++)
    {
        printf("Iteration %d\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("%2d: " PRIe_double "\n", i, vals[i]);
    }
    double result = vals[N+4];
    free(vals);
    return result;
}

int main(void)
{
    double val = M_PI;
    double eps;
    double r;

    eps = 1.0 / 64.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    eps = 1.0 / 100000.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    return(0);
}

输出

在 Mac OS X 10.7.5 上使用 GCC 4.7.1，输出为：

h = 1.562500000000000e-02
-9: deriv(3.000967653589793e+00) = -9.901285883701071e-01
-8: deriv(3.016592653589793e+00) = -9.921976672293290e-01
-7: deriv(3.032217653589793e+00) = -9.940245152582091e-01
-6: deriv(3.047842653589793e+00) = -9.956086864580017e-01
-5: deriv(3.063467653589793e+00) = -9.969497940760287e-01
-4: deriv(3.079092653589793e+00) = -9.980475107000991e-01
-3: deriv(3.094717653589793e+00) = -9.989015683384429e-01
-2: deriv(3.110342653589793e+00) = -9.995117584851364e-01
-1: deriv(3.125967653589793e+00) = -9.998779321710066e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.157217653589793e+00) = -9.998779321710066e-01
 2: deriv(3.172842653589793e+00) = -9.995117584851364e-01
 3: deriv(3.188467653589793e+00) = -9.989015683384429e-01
 4: deriv(3.204092653589793e+00) = -9.980475107000991e-01
 5: deriv(3.219717653589793e+00) = -9.969497940760287e-01
 6: deriv(3.235342653589793e+00) = -9.956086864580017e-01
 7: deriv(3.250967653589793e+00) = -9.940245152582091e-01
 8: deriv(3.266592653589793e+00) = -9.921976672293291e-01
 9: deriv(3.282217653589793e+00) = -9.901285883701071e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -1.091570566584531e-01
 3: -9.361273104952932e-02
 4: -7.804555123490846e-02
 5: -6.245931771829009e-02
 6: -4.685783565504131e-02
 7: -3.124491392324735e-02
 8: -1.562436419383969e-02
 9: -1.776356839400250e-15
10: 1.562436419384087e-02
11: 3.124491392324558e-02
12: 4.685783565504250e-02
13: 6.245931771828653e-02
14: 7.804555123490846e-02
15: 9.361273104953050e-02
16: 1.091570566584471e-01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.577900251527073e+00
 3: 1.660540592568783e+00
 4: 9.969497901146779e-01
 5: 9.980475067341610e-01
 6: 9.989015643694564e-01
 7: 9.995117545137316e-01
 8: 9.998779281977730e-01
 9: 9.999999960264082e-01
10: 9.998779281978486e-01
11: 9.995117545137319e-01
12: 9.989015643693869e-01
13: 9.980475067340947e-01
14: 9.969497901149178e-01
15: 1.660540592568512e+00
16: -3.577900251527123e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 6.553266640232313e+01
 3: 1.898706817407991e+02
 4: -5.267598134705870e+01
 5: 3.608762837830827e+00
 6: 4.685783548517186e-02
 7: 3.124491378285654e-02
 8: 1.562436412277357e-02
 9: 3.226456139297321e-12
10: -1.562436412279785e-02
11: -3.124491378869069e-02
12: -4.685783548888563e-02
13: -3.608762837816180e+00
14: 5.267598134704988e+01
15: -1.898706817408119e+02
16: -6.553266640231030e+01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 8.382087654791743e+03
 3: -5.062815705775390e+03
 4: -7.597917560836845e+03
 5: 3.261984801532626e+03
 6: -4.336626618856812e+02
 7: 1.791410702361821e+01
 8: -9.998779233550453e-01
 9: -9.999999914294619e-01
10: -9.998779238596159e-01
11: 1.791410702341709e+01
12: -4.336626618847606e+02
13: 3.261984801532444e+03
14: -7.597917560838102e+03
15: -5.062815705774387e+03
16: 8.382087654792240e+03
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -9.999999914294619e-01 (eps = 0.015625)
h = 1.000000000000000e-05
-9: deriv(3.141502653589793e+00) = -9.999999959500000e-01
-8: deriv(3.141512653589793e+00) = -9.999999968000000e-01
-7: deriv(3.141522653589793e+00) = -9.999999975500000e-01
-6: deriv(3.141532653589793e+00) = -9.999999982000000e-01
-5: deriv(3.141542653589793e+00) = -9.999999987500000e-01
-4: deriv(3.141552653589793e+00) = -9.999999992000000e-01
-3: deriv(3.141562653589793e+00) = -9.999999995500000e-01
-2: deriv(3.141572653589793e+00) = -9.999999998000000e-01
-1: deriv(3.141582653589793e+00) = -9.999999999500000e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.141602653589793e+00) = -9.999999999500000e-01
 2: deriv(3.141612653589793e+00) = -9.999999998000000e-01
 3: deriv(3.141622653589793e+00) = -9.999999995500000e-01
 4: deriv(3.141632653589793e+00) = -9.999999992000000e-01
 5: deriv(3.141642653589793e+00) = -9.999999987500000e-01
 6: deriv(3.141652653589793e+00) = -9.999999982000000e-01
 7: deriv(3.141662653589793e+00) = -9.999999975500000e-01
 8: deriv(3.141672653589793e+00) = -9.999999968000000e-01
 9: deriv(3.141682653589793e+00) = -9.999999959500000e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -7.000000116589669e-05
 3: -5.999999941330713e-05
 4: -5.000000598739025e-05
 5: -3.999999683331386e-05
 6: -2.999999600591015e-05
 7: -2.000000257999327e-05
 8: -1.000000082740371e-05
 9: 0.000000000000000e+00
10: 1.000000082740371e-05
11: 2.000000165480742e-05
12: 2.999999508072429e-05
13: 3.999999590812801e-05
14: 5.000000413701854e-05
15: 5.999999663774957e-05
16: 6.999999746515327e-05
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.583333244325556e+00
 3: 1.666666318844711e+00
 4: 1.000000128999663e+00
 5: 1.000000691821058e+00
 6: 9.999995738881511e-01
 7: 9.999997049561470e-01
 8: 1.000000198388602e+00
 9: 1.000000075030488e+00
10: 1.000000144419428e+00
11: 9.999996509869722e-01
12: 9.999995893079155e-01
13: 1.000000645561765e+00
14: 1.000000028771196e+00
15: 1.666666187776715e+00
16: -3.583333074708150e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 1.027777535146502e+05
 3: 2.972222191231725e+05
 4: -8.263881528669108e+04
 5: 5.555518108303881e+03
 6: -6.636923522614542e-02
 7: 4.677328483600658e-02
 8: 1.991719546327412e-02
 9: -3.148201866790915e-03
10: -2.319389535987426e-02
11: -4.176186143567406e-02
12: 6.726872146719150e-02
13: -5.555525175695851e+03
14: 8.263880834779718e+04
15: -2.972222015189417e+05
16: -1.027777456120210e+05
17: 0.000000000000000e+00
18: 0.000000000000000e+00

疯狂的迹象...

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 2.050347140226725e+10
 3: -1.240740057099195e+10
 4: -1.858796490195886e+10
 5: 7.986101364079453e+09
 6: -1.059021715700324e+09
 7: 4.630176289959386e+07
 8: -3.687893612405425e+03
 9: -2.136279835945886e+03
10: -2.968840021291520e+03
11: 4.630204773690500e+07
12: -1.059022490436399e+09
13: 7.986100736580715e+09
14: -1.858796331554353e+10
15: -1.240739964045201e+10
16: 2.050347017082775e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -2.136279835945886e+03 (eps = 0.000010)

请注意，在最后一次迭代中，结果如何以数十亿的值变得疯狂。您至少有一个数值稳定性问题，或者您可能需要查看导数公式。请注意，即使具有较大 epsilon 的运行也倾向于在以后的迭代中使用较大的值。

---

使用“官方”derivative()功能

将问题中现在存在的导数函数插入上面的代码中，在第 4 次迭代中会产生更不稳定的答案：

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -6.248925477935563e+09
 3: 1.729549900845405e+11
 4: -2.600544559219368e+11
 5: -2.755286326619338e+11
 6: 8.100546069731433e+11
 7: -2.961111189495415e+09
 8: -7.936480686806423e+11
 9: 2.430177384467434e+11
10: 2.389084910067162e+11
11: -6.461168564124718e+10
12: -4.574822745530297e+10
13: -8.923883451146609e+10
14: 7.042030613792160e+10
15: 4.988386306820556e+10
16: -1.793262395787471e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = 2.430177384467434e+11 (eps = 0.000010)

我想知道数组索引 0、1、17 和 18 处出现的零在多大程度上加剧了这个问题。

Answer 3

这是一个棘手的话题。数值微分是病态的，会导致取消和舍入误差。在您进行重复区分的情况下，该问题被显着放大。

通过使用 m 个点确定特定微分阶 N 的系数，而不是进行重复微分，您将获得更好的高阶导数估计。

例如，就像您似乎使用一阶导数的系数一样：

{1, -8, 8, -1} / (12*h)

使用 5 个点近似二阶导数的系数是：

{-1, 16, -30, 16, -1} / (12*h²)

您可以通过围绕样本点进行泰勒级数展开来求解一般问题的系数，找到这些展开的线性组合，从而得到所需的导数。

对于 m 个点 {x + n[1]*h, x+n[2]*h, x+n[3]*h, ..., x+n[m]*h}，系数 (k) 估计N次导数可以通过以下方程组计算：

M*k=b

其中 M 是 am*m 矩阵：

M[i,j] = n[j]^(i-1), i,j = 1..m

和

b = [0, 0, 0, …, N!/h^N, …, 0, 0, 0],

其中N！表示 N 的阶乘。