How to calculate derivatives using Numpy?

Question

How to calculate the derivative of a function, for example

y = x ² +1

using numpy ?

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

Answer 1

You have four options.

• You can use Final Differences
• You can use Auto derivatives .
• You can use Symbolic differentiation
• You can calculate derivatives manually.

The final differences do not require external tools, but are subject to a numerical error and, if you are in a multidimensional situation, may take some time.

Symbolic differentiation is ideal if your problem is simple enough. Symbolic methods are now becoming quite stable. SymPy is a great project for this that integrates well with NumPy. Take a look at the autowrap or lambdify features or check out Jensen's blog post about a similar question .

Automatic derivatives are very cool, not prone to numerical errors, but require some additional libraries (there are some good options for this). This is the most reliable, but also the most difficult / difficult choice. If you're perfectly limiting 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.] 

Answer 2

NumPy does not provide general functionality for computing derivatives. However, it can handle a simple special case of polynomials:

 >>> 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 calculate the derivative numerically, you can avoid using the central difference coefficients for the vast majority of applications. For a derivative at one point, the formula will look 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 x abscissa array with a corresponding array of y function values, you can calculate the approximations of the derivatives using

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

Answer 3

The easiest way I can think of is to use the numpy gradient function :

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

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

Answer 4

I will put another method on the heap ...

scipy.interpolate many interpolation splines are able to provide derivatives. So, using a linear spline ( k=1 ), the derivative of the spline (using the derivative() method) should be equivalent to the forward difference. I'm not quite sure, but I believe that using the cubic spline derivative will be similar to the centered difference derivative, since it uses the before and after values ​​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) 

Answer 5

Assuming you want to use numpy , you can numerically compute the derivative of a function at any point using the Strict 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 to improve results:

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

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

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

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

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

Answer 6

Depending on the level of accuracy you require, you can fix it yourself using a 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 cannot actually take the limit of the gradient, but its a curious pleasure. You have to be careful because

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

Answer 7

The machine learning community uses Autograd to calculate gradients:

" Efficiently computes derivatives of numpy code ."

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 calculate gradients of complex functions, such as multidimensional functions.