[python] How do I compute derivative using Numpy?

How do I calculate the derivative of a function, for example

y = x2+1

using `numpy`?

Let's say, I want the value of derivative at x = 5...

This question is related to `python` `math` `numpy`

NumPy does not provide general functionality to compute derivatives. It can handles the simple special case of polynomials however:

``````>>> p = numpy.poly1d([1, 0, 1])
>>> print p
2
1 x + 1
>>> q = p.deriv()
>>> print q
2 x
>>> q(5)
10
``````

If you want to compute the derivative numerically, you can get away with using central difference quotients for the vast majority of applications. For the derivative in a single point, the formula would be something like

``````x = 5.0
eps = numpy.sqrt(numpy.finfo(float).eps) * (1.0 + x)
print (p(x + eps) - p(x - eps)) / (2.0 * eps * x)
``````

if you have an array `x` of abscissae with a corresponding array `y` of function values, you can comput approximations of derivatives with

``````numpy.diff(y) / numpy.diff(x)
``````

Depending on the level of precision you require you can work it out yourself, using the simple proof of differentiation:

``````>>> (((5 + 0.1) ** 2 + 1) - ((5) ** 2 + 1)) / 0.1
10.09999999999998
>>> (((5 + 0.01) ** 2 + 1) - ((5) ** 2 + 1)) / 0.01
10.009999999999764
>>> (((5 + 0.0000000001) ** 2 + 1) - ((5) ** 2 + 1)) / 0.0000000001
10.00000082740371
``````

we can't actually take the limit of the gradient, but its kinda fun. You gotta watch out though because

``````>>> (((5+0.0000000000000001)**2+1)-((5)**2+1))/0.0000000000000001
0.0
``````

You can use `scipy`, which is pretty straight forward:

`scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)`

Find the nth derivative of a function at a point.

``````from scipy.misc import derivative

def f(x):
return x**2 + 1

derivative(f, 5, dx=1e-6)
# 10.00000000139778
``````

Assuming you want to use `numpy`, you can numerically compute the derivative of a function at any point using the Rigorous definition:

``````def d_fun(x):
h = 1e-5 #in theory h is an infinitesimal
return (fun(x+h)-fun(x))/h
``````

You can also use the Symmetric derivative for better results:

``````def d_fun(x):
h = 1e-5
return (fun(x+h)-fun(x-h))/(2*h)
``````

Using your example, the full code should look something like:

``````def fun(x):
return x**2 + 1

def d_fun(x):
h = 1e-5
return (fun(x+h)-fun(x-h))/(2*h)
``````

Now, you can numerically find the derivative at `x=5`:

``````In [1]: d_fun(5)
Out[1]: 9.999999999621423
``````

I'll throw another method on the pile...

`scipy.interpolate`'s many interpolating splines are capable of providing derivatives. So, using a linear spline (`k=1`), the derivative of the spline (using the `derivative()` method) should be equivalent to a forward difference. I'm not entirely sure, but I believe using a cubic spline derivative would be similar to a centered difference derivative since it uses values from before and after to construct the cubic spline.

``````from scipy.interpolate import InterpolatedUnivariateSpline

# Get a function that evaluates the linear spline at any x
f = InterpolatedUnivariateSpline(x, y, k=1)

# Get a function that evaluates the derivative of the linear spline at any x
dfdx = f.derivative()

# Evaluate the derivative dydx at each x location...
dydx = dfdx(x)
``````

To install:

``````pip install autograd
``````

Here is an example:

``````import autograd.numpy as np

def fct(x):
y = x**2+1
return y

``````

It can also compute gradients of complex functions, e.g. multivariate functions.

The most straight-forward way I can think of is using numpy's gradient function:

``````x = numpy.linspace(0,10,1000)
dx = x[1]-x[0]
y = x**2 + 1