154

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...

3

9 Answers 9

205

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.]
6
  • 2
    Sorry, if this seems stupid, What is the differences between 3.Symbolic Differentiation and 4.by hand differentiation?? Apr 12, 2012 at 16:55
  • 16
    When I said "symbolic differentiation" I intended to imply that the process was handled by a computer. In principle 3 and 4 differ only by who does the work, the computer or the programmer. 3 is preferred over 4 due to consistency, scalability, and laziness. 4 is necessary if 3 fails to find a solution.
    – MRocklin
    Apr 13, 2012 at 16:51
  • 5
    In line 7 we made f, a function that computes the derivative of y wrt x. In 8 we apply this derivative function to a vector of all ones and get the vector of all twos. This is because, as stated in line 6, yprime = 2*x.
    – MRocklin
    Apr 14, 2012 at 13:45
  • Just for the sake of completeness, you can also do differentiation by integration (see Cauchy's integral formula), it is implemented e.g. in mpmath (not sure however what they exactly do).
    – DerWeh
    Oct 1, 2019 at 2:28
  • Is there an easy way to do finite differences in numpy without implementing it yourself? e.g. I want to find the gradient of a function at predefined points.
    – Alex
    Aug 25, 2020 at 0:07
79

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.

6
  • 2
    What if dx isn't constant?
    – weberc2
    Jul 1, 2015 at 21:22
  • 4
    @weberc2, in that case you should divide one vector by another, but treat the edges separately with forward and backward derivatives manually.
    – Sparkler
    Jul 2, 2015 at 22:06
  • 3
    Or you could interpolate y with a constant dx, then calculate the gradient.
    – IceArdor
    Nov 16, 2016 at 3:43
  • @Sparkler Thanks for your suggestion. If I may ask 2 small questions, (i) why do we pass dx to numpy.gradient instead of x? (ii) Can we also do the last line of yours as follows: dydx = numpy.gradient(y, numpy.gradient(x))?
    – user929304
    Jan 22, 2018 at 11:33
  • 2
    As of v1.13, non uniform spacing can be specified using an array as the second argument. See the Examples section of this page. Mar 29, 2019 at 2:12
36

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)
9
  • 3
    'Computing numerical derivatives for more general case is easy' -- I beg to differ, computing numerical derivatives for general cases is quite difficult. You just chose nicely behaved functions. Mar 26, 2012 at 17:18
  • what does 2 mean after >>>print p ?? (on 2nd line) Mar 26, 2012 at 17:23
  • @DrStrangeLove: That's the exponent. It's meant to simulate mathematical notation. Mar 26, 2012 at 17:26
  • @SvenMarnach is it maximum exponent?? or what?? Why does it think exponent is 2?? We inputted only coefficients... Mar 26, 2012 at 17:29
  • 2
    @DrStrangeLove: The output is supposed to be read as 1 * x**2 + 1. They put the 2 in the line above because it's an exponent. Look at it from a distance. Mar 26, 2012 at 17:31
22

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
3
  • This is a beautiful application of the fundamental definition of the derivative. I've always wondered why I needed to learn that (other than to understand the idea of secant approaching tangent). Well, now I know.
    – rb3652
    Dec 21, 2021 at 0:30
  • 3
    I would add the function as a parameter: def d_func(func, x):
    – LoMaPh
    Jan 9, 2022 at 22:33
  • 1
    @rb3652 First in foremost it is used to derive all the rules of derivatives. Secondly, it is also often used in mathematical proofs. Just to name a few applications.
    – Tera
    Apr 4, 2022 at 15:45
12

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)
2
  • just tried this, i keep getting errors from this function AxisError: axis -1 is out of bounds for array of dimension 0 and I dont see any answers to this on the community either , any help ?
    – Ayan Mitra
    May 2, 2019 at 10:59
  • Post your problem as a new question and link to it here. Providing an example that causes your error to occur will probably be needed. Errors I have with interp functions are usually because the data isn't well formed going in - like repeated values, wrong number of dimensions, one of the arrays is accidentally empty, data isn't sorted against x or when sorted isn't a valid function, etc. It's possible scipy is calling numpy incorrectly, but very unlikely. Check x.shape and y.shape. See if np.interp() works - it may provide a more helpful error if not. May 3, 2019 at 7:06
9

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
1
7

To calculate gradients, the machine learning community uses Autograd:

"Efficiently computes derivatives of numpy code."

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.

2
  • Hi can this function be used to differentiate between two columns of data numerically by providing the step length ? thanks
    – Ayan Mitra
    May 2, 2019 at 11:16
  • 1
    This was answered more than 3 years ago, but autograd is not being developed anymore (just maintained). autograd's Main developers now point you to github.com/google/jax
    – Jose R
    Nov 23, 2021 at 0:17
3

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
1

To compute the derivative of a numerical function, use this second order finite differences scheme as seen in: https://youtu.be/5QnToSn_oxk?t=1804

dx = 0.01
x = np.arange(-4, 4+dx, dx)
y = np.sin(x)
n = np.size(x)

yp = np.zeros(n)
yp[0] = (-3*y[0] + 4*y[1] - y[2]) / (2*dx)
yp[n-1] = (3 * y[n-1] - 4*y[n-2] + y[n-3]) / (2*dx)
for j in range(1,n-1):
    yp[j] = (y[j+1] - y[j-1]) / (2*dx)

Or if you want to use a higher order, use: https://youtu.be/5QnToSn_oxk?t=1374

All that comes from the Nathan Kutz' lectures of the course "Beginning Scientific Computing".

Not the answer you're looking for? Browse other questions tagged or ask your own question.