Calculating gradient in 3D

Question

I have the following set of points in 3d-space and D'd like to calculate the gradient everywhere, i.e. have a vector field returned.

points = []
for i in np.linspace(-20,20,100):   
    for j in np.linspace(-20,20,100):
        points.append([i,j,i**2+j**2])
points = np.array(points)

It's an elliptic paraboloid.

Using np.gradient(points), http://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html

I neither get the correct values nor the dimension I would expect. Can anyone give me a hint?

Answer 1

You are mixing together the indices and the values in 'points', so the gradient is giving you wrong results. Here is a better way to construct the points with numpy and calculate the gradient:

x, y = np.mgrid[-20:20:100j, -20:20:100j]
z = x**2 + y**2
grad = np.gradient(z)

The resulting gradient is a tuple with two arrays, one for the gradient on the first direction, another for the gradient on the second direction. Note that this gradient doesn't take into account the separation between points (ie, delta x and delta y), so to get the derivative you need to divide by it:

deriv = grad/(40./100.)

If you want to reconstruct your 'points' as before, you just need to do:

points = np.array([x.ravel(), y.ravel(), z.ravel()]).T

You may also be interested in numpy's diff function, that gives the discrete difference along a given axis.