10

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. There is anything more powerful in python or i'd better start writing from scratch something else maybe switchin to fortran?

Improve this question
4
  • 1
    Soooo what's the question? What module you should use? Is something going wrong? Commented Jul 27, 2013 at 19:14
  • Question is what do gradient do? Why gives me 2 3d vectors instead of 1 2d vector? Do gradient actually compute really a gradient? By its output i cannot say. Doesn't look exact to me. Commented Jul 27, 2013 at 19:23
  • I thought it was clear, 3rd component of my input is the scalar field, each value on 3rd component is the value of my function for each (x,y) point. Commented Jul 27, 2013 at 19:39
  • the problem is that you are giving gradient the wrong input. It doesn't care about x1,x2 or in the last example i,j. It just wants a matrix of i**2+j**2 values. It's implicit that your matrix of i**2+j**2 values correspond to the xy plane, and the optional scalar arguments of gradient account for step size assumptions, i.e. if your x points are not 1 away from eachother, and the same for your y points. I'm going out of town today but I'll update my answer when I get back tonight. Commented Jul 28, 2013 at 16:18

1 Answer 1

Reset to default
31

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:

sinxpy

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:

gx

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) enter image description here

and the gradients

g2

g1

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.

Improve this answer
Sign up to request clarification or add additional context in comments.

11 Comments

Grid is already in the file i give to numpy.gradient. I cannot calculate the function as you did because i don't know the equation. I got those values from an ab initio calculation. And again: i don't know why should i have 2 3d vectors when for a 2 variables function gradient is a 2d vector...
My function was part of the contrived example, a means to generate Z data. You say in your question you have x,y,z data already, so you don't need a function...just make Z appropriately with your frequency data. gx and gy give x and y derivatives for every point in Z, which should be an X by Y matrix.
I already tried to do what you did in your example but gradient still give me 2 3-d vectors... Ab initio calculation gives me 8 frequencies for each point (it's a phonon dispersion relation) grid is a sampling of the brillouin zone. Don't know how to make Z appropriately, it should be already appropriate... Maybe gradient simply don't calculate gradient? I tryed with x2+y2, gradiend should be 2d vector: (2x,2y), but i still get 2 3-d vectors. How can i obtain the real gradient from such output?
I think your confusion comes from your expectation that gradient returns a vector of functions, when it cannot. numpy doesn't do algebra like mathematica or maple. You give gradient a matrix of your Z values and it computes step-wise the slope between each X,X+1 and Y,Y+1, giving you for every point the rate of change of x and y at that point. i.e. it approximates the evaluation of f(x,y) = (2x,2y).
Why 3 components then? There are 2 directions so i would expect a 2 component vector...What should i do to get gradient? I didn't expect to get 2x,2y as output, but vectors with such components... in point (2,2) i would expect a vector (4,4), in point 1,3 i would expect (2,6) and so on. All i need is gradient module on each point, but i don't know how to work it out from numpy.gradient output... Have you got any suggestion?
|

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.