I want to write a subroutine that calculates the nth derivative of a function given by the subroutine of the form:
double func(double x)
// a mathematical expression representing a smooth curve
//returns the value of the function at x
double nthderive(double (*f)(double),double x0,int n)
const double delta=1.0e-8;
double x1=x0 - delta;
double x2=x0 + delta;
Not counting the numerical instability issue, that makes adjustment of
Delta difficult and precludes the use for high order derivatives (say
n > 3 !), the recursive solution is pretty inefficient, as it takes
2^n function evaluations, where
n+1 are enough.
Indeed, you compute
f' = f(+1) - f(-1) f'' = f'(+1) - f'(-1) = f(+2) - f(0) - f(0) + f(-2) f''' = f''(+1) - f''(-1) = f(+3) - f(+1) - f(+1) + f(-1) - f(+1) + f(-1) + f(-1) - f(-3) ...
when the Binomial formula tells you
f' = f(+1) - f(-1) f'' = f(+2) - 2f(0) + f(-2) f''' = f(+3) - 3f(+1) + 3f(-1) - f(-3) ...