153

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

nbro
  • 15,395
  • 32
  • 113
  • 196
DrStrangeLove
  • 11,227
  • 16
  • 59
  • 72
  • 6
    You need to use Sympy: http://sympy.org/en/index.html Numpy is a numeric computation library for Python – prrao Mar 26 '12 at 16:55
  • 1
    Alternatively, do you want a method for estimating the numerical value of the derivative? For this you can use a finite difference method, but bear in mind they tend to be horribly noisy. – Henry Gomersall Mar 26 '12 at 17:11
  • [Finite-difference approximation of the gradient of a scalar function](https://het.as.utexas.edu/HET/Software/Scipy/generated/scipy.optimize.approx_fprime.html) – JeeyCi Aug 06 '23 at 19:47

9 Answers9

201

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.]
Nathan Davis
  • 5,636
  • 27
  • 39
MRocklin
  • 55,641
  • 23
  • 163
  • 235
  • 2
    Sorry, if this seems stupid, What is the differences between 3.Symbolic Differentiation and 4.by hand differentiation?? – DrStrangeLove Apr 12 '12 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 '12 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 '12 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 01 '19 at 02: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 '20 at 00:07
  • Theano is kindof dead. Look at pytorch or JAX or basically any "deep learning" toolbox. Autodiff is one of its pillars. – visoft Feb 11 '21 at 10:43
77

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.

Sparkler
  • 2,581
  • 1
  • 22
  • 41
  • 2
    What if dx isn't constant? – weberc2 Jul 01 '15 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 02 '15 at 22:06
  • 3
    Or you could interpolate y with a constant dx, then calculate the gradient. – IceArdor Nov 16 '16 at 03: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 '18 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](https://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html#numpy.gradient). – Nathaniel Jones Mar 29 '19 at 02:12
  • `edge_order=2` is more accurate than the default `edge_order=1`, and hardly slower; try `x = np.linspace( 0, 2, 21 ); y = x**3; g = np.gradient( y, x, edge_order=1 ); print( "np.gradient( x^3 ) - 3x^2: \n", g - 3 * x**2 )` – denis Jul 06 '22 at 15:38
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)
Sven Marnach
  • 574,206
  • 118
  • 941
  • 841
  • 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. – High Performance Mark Mar 26 '12 at 17:18
  • what does 2 mean after >>>print p ?? (on 2nd line) – DrStrangeLove Mar 26 '12 at 17:23
  • @DrStrangeLove: That's the exponent. It's meant to simulate mathematical notation. – Sven Marnach Mar 26 '12 at 17:26
  • @SvenMarnach is it maximum exponent?? or what?? Why does it think exponent is 2?? We inputted only coefficients... – DrStrangeLove Mar 26 '12 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. – Sven Marnach Mar 26 '12 at 17:31
  • @SvenMarnach now, i see..) But what if i had x^3+1 and the same coefficients numpy.poly1d([1, 0, 1]) will it give me 1x+x^3+1?? How does it know the right exponent?? – DrStrangeLove Mar 26 '12 at 17:36
  • @DrStrangeLove: To get that, you'd need to use `numpy.poly1d([1, 0, 0, 1])`. – Sven Marnach Mar 26 '12 at 17:40
  • @SvenMarnach Thanks!! Earlier, you said that if i want to get 2*x result - "i can't" get it. but in your example 2*x is there! How come?? – DrStrangeLove Mar 26 '12 at 18:41
  • @DrStrangeLove: You asked for derivatives in general, and you can't do symbolic derivatives with NumPy in general. You can merely do it for the trivial case of polynomials. – Sven Marnach Mar 26 '12 at 19:04
21

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
fabda01
  • 3,384
  • 2
  • 31
  • 37
  • 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 '21 at 00:30
  • 2
    I would add the function as a parameter: `def d_func(func, x):` – LoMaPh Jan 09 '22 at 22:33
  • @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 04 '22 at 15:45
11

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)
flutefreak7
  • 2,321
  • 5
  • 29
  • 39
  • 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 02 '19 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. – flutefreak7 May 03 '19 at 07: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
johnson
  • 3,729
  • 3
  • 31
  • 32
  • 1
    `scipy.misc.derivative` is deprecated from SciPy 1.10 https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html – divenex Jan 18 '23 at 10:58
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.

Gordon Schücker
  • 463
  • 1
  • 6
  • 7
  • Hi can this function be used to differentiate between two columns of data numerically by providing the step length ? thanks – Ayan Mitra May 02 '19 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 https://github.com/google/jax – Jose R Nov 23 '21 at 00: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
fraxel
  • 34,470
  • 11
  • 98
  • 102
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".

Manuel
  • 478
  • 8
  • 20