# mathematical programming – Precise algorithm for finding higher order derivatives

I’m trying to make an algorithm that finds the first 10 or so terms of a function’s Taylor series, which requires finding the nth derivative of the function for the nth term. It’s easy to implement derivatives by following the definition of the derivative:
$$f'(x) = lim_{hto0}dfrac{f(x+h)-f(x)}{h}$$
implemented here in Python:

dx = 0.001
def derivative(f, x):
return (f(x + dx) - f(x)) / dx

The value seems to be even closer to the actual value of the derivative if we define it like this:

dx = 0.001
def derivative(f, x):
return (f(x + dx) - f(x - dx)) / (2 * dx)

which just returns the average of (f(x + dx) - f(x)) / dx and (f(x) - f(x - dx)) / dx.

For higher order derivatives, I implemented a simple recursive function:

dx = 0.001
def nthDerivative(f, n, x):
if n == 0:
return f(x)
return (derivative(f, n - 1, x + dx) - derivative(f, n - 1, x - dx)) / (2 * dx)

I tested the higher order derivatives of $$f$$ at $$1$$, where $$f(x)=x^9$$, and as can be proved by induction,
$$dfrac{d^n}{dx^n}(x^k)=dfrac{k!}{(k – n)!}x^{k-n}$$

Therefore, the nth derivative of $$f$$ at $$1$$ is $$dfrac{9!}{(9 – n)!}$$.
Here are the values returned by the function for n ranging from 0 to 9:

n             Value  Intended value
-----------------------------------
0             1.000               1
1             9.000               9
2            72.001              72
3           504.008             504
4          3024.040            3024
5         15120.252           15120
6         60437.602           60480
7         82298.612          181440
8      32278187.177          362880
9   95496943657.736          362880

As you can see, the values are waaaay off for $$n$$ greater than $$5$$.

What can I do to get closer to the actual values? And is there an algorithm for this that doesn’t have $$O(2^n)$$ performance like mine?