GSL 求导数：gsl_deriv_central

GSL 手册里有函数说明和使用范例，函数说明如下：
int gsl_deriv_central(const gsl_function * f, double x, double h, double * result, double * abserr)
This function computes the numerical derivative of the function f at the point x using an adaptive central
difference algorithm with a step-size of h. The derivative is returned in result and an estimate of its absolute
error is returned in abserr.
The initial value of h is used to estimate an optimal step-size, based on the scaling of the truncation error and
round-off error in the derivative calculation. The derivative is computed using a 5-point rule for equally spaced
abscissae at x ? h, x ? h/2, x, x + h/2, x + h, with an error estimate taken from the difference between the
5-point rule and the corresponding 3-point rule x ? h, x, x + h. Note that the value of the function at x does not
contribute to the derivative calculation, so only 4-points are actually used.

下面是手册里的示范代码：

#include 
#include 
#include 
double f (double x, void * params)
{
        (void)(params); /* avoid unused parameter warning */
        return pow (x, 1.5);
}
        int
main (void)
{
        gsl_function F;
        double result, abserr;
        F.function = &f;
        F.params = 0;
        printf ("f(x) = x^(3/2)\n");
        gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
        printf ("x = 2.0\n");
        printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
        printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0));
        gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
        printf ("x = 0.0\n");
        printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
        printf ("exact = %.10f\n", 0.0);
        return 0;
}
~      

存为 test_gsl_deriv_central.cpp, 然后编译运行：

g++ test_gsl_deriv_central.cpp -lgsl -lgslcblas
./a.out

得到结果：

f(x) = x^(3/2)
x = 2.0
f'(x) = 2.1213203120 +/- 0.0000005006
exact = 2.1213203436

x = 0.0
f'(x) = 0.0000000160 +/- 0.0000000339
exact = 0.0000000000