Search code examples
pythonnumpyderivativedifferentiation

non-uniform spacing with numpy.gradient

I'm not sure how to specify non-uniform spacing when using numpy.gradient.

Here's some example code for y = x**2.

import numpy as np
import matplotlib.pyplot as plt

x = [0.0, 2.0, 4.0, 8.0, 16.0]
y = [0.0, 4.0, 16.0, 64.0, 256.0]
dydx = [0.0, 4.0, 8.0, 16.0, 32.0] # analytical solution

spacing = [0.0, 2.0, 2.0, 4.0, 8.0] #added a zero at the start to get length matching up with y

m = np.gradient(y, spacing)

plt.plot(x, y, 'bo',
         x, dydx, 'r-', #analytical solution
         x, m, 'ro')    #calculated solution
plt.show()

The length of the spacing array will always be one less than the array I want to calculated the gradient of. Adding in a zero to get the lengths matching up (like in the example code above) gives incorrect answers, with an infinite gradient for one point.

I can't understand / follow the numpy.gradient documentation for non-uniform spacing (https://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html)

How should I specify the spacing between points? Is there an alternative way of doing this?

Numpy version 1.9.2

Solution

  • The API of the function is quite confusing. For non-uniformly spaced sample points, the gradient function takes the coordinates of the point rather than the spacings:

    varargs : list of scalar or array, optional

    Spacing between f values. Default unitary spacing for all dimensions. Spacing can be specified using:

    1. single scalar to specify a sample distance for all dimensions.
    2. N scalars to specify a constant sample distance for each dimension. i.e. dx, dy, dz, …
    3. N arrays to specify the coordinates of the values along each dimension of F. The length of the array must match the size of the corresponding dimension
    4. Any combination of N scalars/arrays with the meaning of 2. and 3.

    I slightly modified your example:

    import numpy as np
import matplotlib.pyplot as plt

x = np.random.rand(10)
x.sort()
y = x**2
dydx = 2*x

dydx_grad = np.gradient(y, x)

plt.plot(x, dydx, 'k-', label='analytical solution')
plt.plot(x, dydx_grad, 'ro', label='calculated solution')
plt.legend(); plt.xlabel('x'); plt.ylabel('dy / dx'); plt.show();