Gradient calculation with python

Question

I would like to know how does numpy.gradient work. I used gradient to try to calculate group velocity (group velocity of a wave packet is the derivative of frequencies respect to wavenumbers, not a group of velocities). I fed a 3 column array to it, the first 2 colums are x and y coords, the third column is the frequency of that point (x,y). I need to calculate gradient and I did expect a 2d vector, being gradient definition

df/dx*i+df/dy*j+df/dz*k 

and my function only a function of x and y i did expect something like

df/dx*i+df/dy*j 

But i got 2 arrays with 3 colums each, i.e. 2 3d vectors; at first i thought that the sum of the two would give me the vector i were searchin for but the z component doesn't vanish. I hope i've been sufficiently clear in my explanation. I would like to know how numpy.gradient works and if it's the right choice for my problem. Otherwise i would like to know if there's any other python function i can use.

What i mean is: I want to calculate gradient of an array of values:

data=[[x1,x2,x3]...[x1,x2,x3]]

where x1,x2 are point coordinates on an uniform grid (my points on the brillouin zone) and x3 is the value of frequency for that point. I give in input also steps for derivation for the 2 directions:

stepx=abs(max(unique(data[:,0])-min(unique(data[:,0]))/(len(unique(data[:,0]))-1)

the same for y direction. I didn't build my data on a grid, i already have a grid and this is why kind examples given here in answers do not help me. A more fitting example should have a grid of points and values like the one i have:

data=[]
for i in range(10):
  for j in range(10):
    data.append([i,j,i**2+j**2])

data=array(data,dtype=float)

gx,gy=gradient(data)

another thing i can add is that my grid is not a square one but has the shape of a polygon being the brillouin zone of a 2d crystal.

I've understood that numpy.gradient works properly only on a square grid of values, not what i'm searchin for. Even if i make my data as a grid that would have lots of zeroes outside of the polygon of my original data, that would add really high vectors to my gradient affecting (negatively) the precision of calculation. This module seems to me more a toy than a tool, it has severe limitations imho.

Problem solved using dictionaries.

Answer 1

You need to give gradient a matrix that describes your angular frequency values for your (x,y) points. e.g.

def f(x,y):
    return np.sin((x + y))
x = y = np.arange(-5, 5, 0.05)
X, Y = np.meshgrid(x, y)
zs = np.array([f(x,y) for x,y in zip(np.ravel(X), np.ravel(Y))])
Z = zs.reshape(X.shape)

gx,gy = np.gradient(Z,0.05,0.05)

You can see that plotting Z as a surface gives:

Here is how to interpret your gradient:

gx is a matrix that gives the change dz/dx at all points. e.g. gx[0][0] is dz/dx at (x0,y0). Visualizing gx helps in understanding:

Since my data was generated from f(x,y) = sin(x+y) gy looks the same.

Here is a more obvious example using f(x,y) = sin(x)...

f(x,y)

and the gradients

update Let's take a look at the xy pairs.

This is the code I used:

def f(x,y):
    return np.sin(x)
x = y = np.arange(-3,3,.05)
X, Y = np.meshgrid(x, y)
zs = np.array([f(x,y) for x,y in zip(np.ravel(X), np.ravel(Y))])
xy_pairs = np.array([str(x)+','+str(y) for x,y in zip(np.ravel(X), np.ravel(Y))])
Z = zs.reshape(X.shape)
xy_pairs = xy_pairs.reshape(X.shape)

gy,gx = np.gradient(Z,.05,.05)

Now we can look and see exactly what is happening. Say we wanted to know what point was associated with the value atZ[20][30]? Then...

>>> Z[20][30]
-0.99749498660405478

And the point is

>>> xy_pairs[20][30]
'-1.5,-2.0'

Is that right? Let's check.

>>> np.sin(-1.5)
-0.99749498660405445

Yes.

And what are our gradient components at that point?

>>> gy[20][30]
0.0
>>> gx[20][30]
0.070707731517679617

Do those check out?

dz/dy always 0 check. dz/dx = cos(x) and...

>>> np.cos(-1.5)
0.070737201667702906

Looks good.

You'll notice they aren't exactly correct, that is because my Z data isn't continuous, there is a step size of 0.05 and gradient can only approximate the rate of change.