# Testing derivatives using numerical differentiation

I've been working recently on the AMPL interface for the GNU Scientific Library. Adding a function to AMPL normally requires providing first and second partial derivatives for it. Computing the derivatives is straightforward and there are plenty of tools, such as Wolfram Alpha, that automate this. However, transferring this to code is a tedious and error-prone process.

For example, consider the function gsl_sf_bessel_I0_scaled(double x) which is defined as $$e^{-|x|} I_0(x)$$, where $$I_0(x)$$ is the regular modified cylindrical Bessel function of zeroth order. It has a simple derivative:

$$e^{-|x|} I_1(x)- \frac{x}{|x|} e^{-|x|} I_0(x)$$

The derivative can be expressed even in a more compact form using $$I_0(x)$$ and $$I_1(x)$$ scaled by $$e^{-|x|}$$, gsl_sf_bessel_I0_scaled and gsl_sf_bessel_I1_scaled.

Transferring this into C results in the following small program:

Let's compile and run the program:

$gcc -pedantic -Wall -Wextra test.c -lgsl -lgslcblas -otest$ ./test
-0.332111


It is easy to check that the printed value is incorrect, -0.148394 should be printed instead. Can you spot the error?

$$f'(x) \approx \frac{f(x + h) - f(x)}{h}$$ for some small $$h$$
Although trivial to implement, this method suffers from several issues such as the difficulty of selecting $$h$$ especially for $$x = 0$$. So I used the Ridders' method described in the book Numerical Recipes in C: The Art of Scientific Computing. This methods starts from large $$h$$ and successfully decreases it estimating error at each iteration.