Last Update: Feb. 9, 2006
/***********************************************/
/* GSL sample program */
/* Numerical Differentiation */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
/* Return relative and absolute errors */
double rel_abs_error(double approx_val, double true_val, double *abserr)
{
static double abs_error, rel_error;
abs_error = fabs(approx_val - true_val);
rel_error = abs_error;
if(true_val >= 1.0e-15)
rel_error /= fabs(true_val);
if(abserr != NULL)
*abserr = abs_error;
return rel_error;
}
/* Original function: f(x) = sin(x) */
double func(double x, void *param)
{
return sin(x);
}
/* Main function */
int main(void)
{
int i;
double x, h;
double results[3], rel_errors[3], abs_errors[3], gsl_abs_errors[3];
/* Definition of gsl_function: */
/* struct gsl_function_struct { */
/* double (* function) (double x, void * params); */
/* void * params; */
/* }; */
gsl_function f;
/* Set the original function */
f.function = &func;
f.params = NULL;
/* Get values of f'(x) by Central diff, absolute and relative errors in [0, 2*PI] */
printf(" x Central Diff True f'(x) gsl_abs Abs_err Rel_err\n");
/* Initial guess for seeking the optimal stepsize */
h = 1;
for(i = 0; i <= 20; i++)
{
x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;
/* Central Difference with 5-points */
gsl_deriv_central(&f, x, h, &results[0], &gsl_abs_errors[0]);
rel_errors[0] = rel_abs_error(results[0], cos(x), &abs_errors[0]);
printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[0], cos(x), gsl_abs_errors[0], abs_errors[0], rel_errors[0]);
}
printf("\n");
/* Get values of f'(x) by Forward diff, absolute and relative errors in [0, 2*PI] */
printf(" x Forward Diff True f'(x) gsl_abs Abs_err Rel_err\n");
/* Initial guess for seeking the optimal stepsize */
h = 1;
for(i = 0; i <= 20; i++)
{
x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;
/* Forward Difference with 5-points */
gsl_deriv_forward(&f, x, h, &results[1], &gsl_abs_errors[1]);
rel_errors[1] = rel_abs_error(results[1], cos(x), &abs_errors[1]);
printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[1], cos(x), gsl_abs_errors[1], abs_errors[1], rel_errors[1]);
}
printf("\n");
/* Get values of f'(x) by Backward diff, absolute and relative errors in [0, 2*PI] */
printf(" x Backward Diff True f'(x) gsl_abs Abs_err Rel_err\n");
/* Initial guess for seeking the optimal stepsize */
h = 1;
for(i = 0; i <= 20; i++)
{
x = 0.0 + (double)i * (2.0 * M_PI) / 20.0;
/* Backward Difference with 5-points */
gsl_deriv_backward(&f, x, h, &results[2], &gsl_abs_errors[2]);
rel_errors[2] = rel_abs_error(results[2], cos(x), &abs_errors[2]);
printf("%9.2e %24.17e %24.17e %7.1e %7.1e %7.1e\n", x, results[2], cos(x), gsl_abs_errors[2], abs_errors[2], rel_errors[2]);
}
return 0;
}
$ ./diff
x Central Diff True f'(x) gsl_abs Abs_err Rel_err
0.00e+00 1.00000000000000000e+00 1.00000000000000000e+00 3.1e-11 0.0e+00 0.0e+00
3.14e-01 9.51056516293700027e-01 9.51056516295153531e-01 5.8e-11 1.5e-12 1.5e-12
6.28e-01 8.09016994378582099e-01 8.09016994374947451e-01 8.6e-11 3.6e-12 4.5e-12
9.42e-01 5.87785252296237792e-01 5.87785252292473137e-01 9.6e-11 3.8e-12 6.4e-12
1.26e+00 3.09016994371405118e-01 3.09016994374947451e-01 7.7e-11 3.5e-12 1.1e-11
1.57e+00 1.29526019539601580e-16 6.12323399573676604e-17 1.1e-15 6.8e-17 6.8e-17
1.88e+00 -3.09016994372883880e-01 -3.09016994374947340e-01 7.9e-11 2.1e-12 2.1e-12
2.20e+00 -5.87785252276747383e-01 -5.87785252292473026e-01 7.1e-11 1.6e-11 1.6e-11
2.51e+00 -8.09016994387561250e-01 -8.09016994374947340e-01 1.1e-10 1.3e-11 1.3e-11
2.83e+00 -9.51056516326770796e-01 -9.51056516295153531e-01 1.0e-10 3.2e-11 3.2e-11
3.14e+00 -9.99999999991971755e-01 -1.00000000000000000e+00 4.7e-11 8.0e-12 8.0e-12
3.46e+00 -9.51056516315975542e-01 -9.51056516295153753e-01 1.0e-10 2.1e-11 2.1e-11
3.77e+00 -8.09016994366233533e-01 -8.09016994374947562e-01 8.5e-11 8.7e-12 8.7e-12
4.08e+00 -5.87785252327138186e-01 -5.87785252292473248e-01 1.3e-10 3.5e-11 3.5e-11
4.40e+00 -3.09016994375665599e-01 -3.09016994374947562e-01 8.5e-11 7.2e-13 7.2e-13
4.71e+00 -9.25185853854297035e-17 -1.83697019872102969e-16 1.1e-15 9.1e-17 9.1e-17
5.03e+00 3.09016994363356501e-01 3.09016994374947229e-01 6.9e-11 1.2e-11 3.8e-11
5.34e+00 5.87785252278960058e-01 5.87785252292472915e-01 8.9e-11 1.4e-11 2.3e-11
5.65e+00 8.09016994360889807e-01 8.09016994374947340e-01 9.5e-11 1.4e-11 1.7e-11
5.97e+00 9.51056516251071793e-01 9.51056516295153531e-01 2.2e-11 4.4e-11 4.6e-11
6.28e+00 1.00000000000613243e+00 1.00000000000000000e+00 6.1e-11 6.1e-12 6.1e-12
x Forward Diff True f'(x) gsl_abs Abs_err Rel_err
0.00e+00 9.99999999999996669e-01 1.00000000000000000e+00 4.0e-14 3.3e-15 3.3e-15
3.14e-01 9.51056510849759618e-01 9.51056516295153531e-01 9.7e-08 5.4e-09 5.7e-09
6.28e-01 8.09017012678312830e-01 8.09016994374947451e-01 2.1e-07 1.8e-08 2.3e-08
9.42e-01 5.87785249964989376e-01 5.87785252292473137e-01 2.7e-07 2.3e-09 4.0e-09
1.26e+00 3.09016983129187905e-01 3.09016994374947451e-01 3.1e-07 1.1e-08 3.6e-08
1.57e+00 8.40261297783238669e-09 6.12323399573676604e-17 3.4e-07 8.4e-09 8.4e-09
1.88e+00 -3.09017001982583472e-01 -3.09016994374947340e-01 3.1e-07 7.6e-09 7.6e-09
2.20e+00 -5.87785245014908586e-01 -5.87785252292473026e-01 2.8e-07 7.3e-09 7.3e-09
2.51e+00 -8.09016952336751749e-01 -8.09016994374947340e-01 2.4e-07 4.2e-08 4.2e-08
2.83e+00 -9.51056528962701320e-01 -9.51056516295153531e-01 1.2e-07 1.3e-08 1.3e-08
3.14e+00 -9.99999962493470784e-01 -1.00000000000000000e+00 3.7e-08 3.8e-08 3.8e-08
3.46e+00 -9.51056548024574067e-01 -9.51056516295153753e-01 1.4e-07 3.2e-08 3.2e-08
3.77e+00 -8.09017011907626538e-01 -8.09016994374947562e-01 2.2e-07 1.8e-08 1.8e-08
4.08e+00 -5.87785197956959493e-01 -5.87785252292473248e-01 2.1e-07 5.4e-08 5.4e-08
4.40e+00 -3.09017005732066974e-01 -3.09016994374947562e-01 3.3e-07 1.1e-08 1.1e-08
4.71e+00 -6.16060791038268051e-09 -1.83697019872102969e-16 3.4e-07 6.2e-09 6.2e-09
5.03e+00 3.09016957995891850e-01 3.09016994374947229e-01 3.6e-07 3.6e-08 1.2e-07
5.34e+00 5.87785191208399227e-01 5.87785252292472915e-01 3.4e-07 6.1e-08 1.0e-07
5.65e+00 8.09016941989682770e-01 8.09016994374947340e-01 2.5e-07 5.2e-08 6.5e-08
5.97e+00 9.51056551687453067e-01 9.51056516295153531e-01 9.7e-08 3.5e-08 3.7e-08
6.28e+00 9.99999948293885943e-01 1.00000000000000000e+00 5.6e-08 5.2e-08 5.2e-08
x Backward Diff True f'(x) gsl_abs Abs_err Rel_err
0.00e+00 9.99999999999996669e-01 1.00000000000000000e+00 4.0e-14 3.3e-15 3.3e-15
3.14e-01 9.51056515947545700e-01 9.51056516295153531e-01 1.3e-07 3.5e-10 3.7e-10
6.28e-01 8.09016975541636763e-01 8.09016994374947451e-01 2.1e-07 1.9e-08 2.3e-08
9.42e-01 5.87785266353461622e-01 5.87785252292473137e-01 2.6e-07 1.4e-08 2.4e-08
1.26e+00 3.09016987824628775e-01 3.09016994374947451e-01 3.3e-07 6.6e-09 2.1e-08
1.57e+00 -8.40261297783236518e-09 6.12323399573676604e-17 3.4e-07 8.4e-09 8.4e-09
1.88e+00 -3.09017003676207824e-01 -3.09016994374947340e-01 3.3e-07 9.3e-09 9.3e-09
2.20e+00 -5.87785279428411078e-01 -5.87785252292473026e-01 2.9e-07 2.7e-08 2.7e-08
2.51e+00 -8.09016949019604303e-01 -8.09016994374947340e-01 1.5e-07 4.5e-08 4.5e-08
2.83e+00 -9.51056555594841480e-01 -9.51056516295153531e-01 1.4e-07 3.9e-08 3.9e-08
3.14e+00 -9.99999962493472783e-01 -1.00000000000000000e+00 3.7e-08 3.8e-08 3.8e-08
3.46e+00 -9.51056513728432695e-01 -9.51056516295153753e-01 1.3e-07 2.6e-09 2.6e-09
3.77e+00 -8.09017008244804159e-01 -8.09016994374947562e-01 1.8e-07 1.4e-08 1.4e-08
4.08e+00 -5.87785171062316869e-01 -5.87785252292473248e-01 3.6e-07 8.1e-08 8.1e-08
4.40e+00 -3.09016979155126492e-01 -3.09016994374947562e-01 3.3e-07 1.5e-08 1.5e-08
4.71e+00 6.16060791038269127e-09 -1.83697019872102969e-16 3.4e-07 6.2e-09 6.2e-09
5.03e+00 3.09016965742043959e-01 3.09016994374947229e-01 2.9e-07 2.9e-08 9.3e-08
5.34e+00 5.87785213112425931e-01 5.87785252292472915e-01 2.3e-07 3.9e-08 6.7e-08
5.65e+00 8.09017033908107752e-01 8.09016994374947340e-01 2.4e-07 4.0e-08 4.9e-08
5.97e+00 9.51056439228550321e-01 9.51056516295153531e-01 6.6e-08 7.7e-08 8.1e-08
6.28e+00 9.99999948293890162e-01 1.00000000000000000e+00 5.6e-08 5.2e-08 5.2e-08