How to Compute Derivative Using Numpy

How do I compute derivative using Numpy?

You have four options

1. Finite Differences
2. Automatic Derivatives
3. Symbolic Differentiation
4. Compute derivatives by hand.

Finite differences require no external tools but are prone to numerical error and, if you're in a multivariate situation, can take a while.

Symbolic differentiation is ideal if your problem is simple enough. Symbolic methods are getting quite robust these days. SymPy is an excellent project for this that integrates well with NumPy. Look at the autowrap or lambdify functions or check out Jensen's blogpost about a similar question.

Automatic derivatives are very cool, aren't prone to numeric errors, but do require some additional libraries (google for this, there are a few good options). This is the most robust but also the most sophisticated/difficult to set up choice. If you're fine restricting yourself to numpy syntax then Theano might be a good choice.

Here is an example using SymPy

In [1]: from sympy import *
In [2]: import numpy as np
In [3]: x = Symbol('x')
In [4]: y = x**2 + 1
In [5]: yprime = y.diff(x)
In [6]: yprime
Out[6]: 2⋅x

In [7]: f = lambdify(x, yprime, 'numpy')
In [8]: f(np.ones(5))
Out[8]: [ 2.  2.  2.  2.  2.]

Calculating the derivative of points in python

The idea behind what you are trying to do is correct, but there are a couple of points to make it work as intended:

1. There is a typo in calc_diffY(X), the derivative of X**2 is 2*X, not 2**X:

  def calc_diffY(X):    
      yval_dash = 3*(X**2) + 2*X

By doing this you don't obtain much better results:

yval_dash = [5, 16, 33, 56, 85]
numpyDiff = [10. 24. 44. 70.]

2. To calculate the numerical derivative you should do a "Difference quotient" which is an approximation of a derivative

numpyDiff = np.diff(yval)/np.diff(xval)

The approximation becomes better and better if the values of the points are more dense.
The difference between your points on the x axis is 1, so you end up in this situation (in blue the analytical derivative, in red the numerical):

If you reduce the difference in your x points to 0.1, you get this, which is much better:

Just to add something to this, have a look at this image showing the effect of reducing the distance of the points on which the derivative is numerically calculated, taken from Wikipedia:

Numpy or SciPy Derivative function for non-uniform spacing?

You can create your own functions using numpy. For the derivatives using forward differences (edit thanks to @EOL, but note that NumPy's diff() is not a differentiate function):

def diff_fwd(x, y): 
    return np.diff(y)/np.diff(x)

"central" differences (it is not necessarily central, depending on your data spacing):

def diff_central(x, y):
    x0 = x[:-2]
    x1 = x[1:-1]
    x2 = x[2:]
    y0 = y[:-2]
    y1 = y[1:-1]
    y2 = y[2:]
    f = (x2 - x1)/(x2 - x0)
    return (1-f)*(y2 - y1)/(x2 - x1) + f*(y1 - y0)/(x1 - x0)

where y contains the function evaluations and x the corresponding "times", such that you can use an arbitrary interval.

Get derivative of data in python

You have a loss of significance problem. It means that when adding a large floating point number to a small one, the precision of the small one is partially lost as a numpy double can only hold 64 bits of information.

To solve this issue you have to make sure that the scale of the numbers you add/multiply/divide is not too different. One simple solution is dividing x by 1e9 or multiplying h by 1e9. If you do this you get essentially the same precision as in your example.

Also if x has high enough resolution a simpler way to numerically differentiate the function would be der = np.diff(y) / np.diff(x). This way you do not have to worry about setting h. However in this case note that dy is one element shorter that y, and dy[i] is actually an approximation of the derivative at `(x[i] + x[i+1]) / 2. So to plot it you would do:

der = np.diff(y) / np.diff(x)
x2 = (x[:-1] + x[1:]) / 2
plt.plot(x2, der, 'r', x, y, 'g', x, -np.sin(x),'b')

Partial derivative in Python

np.diff might be the most idiomatic numpy way to do this:

y = np.empty_like(x)
y[:-1] = np.diff(x, axis=0) / dx
y[-1] = -x[-1] / dx

You may also be interested in np.gradient, although this function takes the gradient over all dimensions of the input array rather than a single one.

Why the Derivative of a function computed using scipy is varying with the value of parameter dx?

As for the documentation, the function scipy.misc.derivative.

Ultimately it is calculating

sum(weights[k]*f(x0+(k-order//2)*dx) for k in range(order))/dx**n

For small values it converges to the derivative, but in practice if you use too small values you the roundoff errors will dominate. So you have to find a reasonable value of dx for your purpose. And, as you already know you can calculate the derivatives symbolically.

---

Related Topics