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.

In your case:

```
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 calculate gradients, the machine learning community uses Autograd:

To install:

```
pip install autograd
```

Here is an example:

```
import autograd.numpy as np
from autograd import grad
def fct(x):
y = x**2+1
return y
grad_fct = grad(fct)
print(grad_fct(1.0))
```

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
dydx = numpy.gradient(y, dx)
```

This way, dydx will be computed using central differences and will have the same length as y, unlike numpy.diff, which uses forward differences and will return (n-1) size vector.

Source: Stackoverflow.com