python - 用python计算梯度

Question

我想知道 numpy.gradient 是如何工作的。我使用梯度来尝试计算群速度(波包的群速度是频率相对于波数的导数，而不是一组速度)。我给它提供了一个 3 列数组，前两列是 x 和 y 坐标，第三列是该点 (x,y) 的频率。我需要计算梯度，我确实期望一个 2d 向量，梯度定义

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

我的函数只是 x 和 y 的函数，我确实期待类似的东西

df/dx*i+df/dy*j 

但是我得到了 2 个数组，每个数组有 3 个列，即 2 个 3d 向量；起初我以为两者的总和会给我我正在搜索的向量，但 z 分量不会消失。我希望我的解释已经足够清楚了。我想知道 numpy.gradient 是如何工作的，以及它是否适合我的问题。否则我想知道是否有任何其他 python 函数我可以使用。

我的意思是:我想计算值数组的梯度:

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

其中 x1、x2 是均匀网格上的点坐标(我在布里渊区上的点)，x3 是该点的频率值。我还输入了 2 个方向的推导步骤:

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

y 方向相同。我没有在网格上构建我的数据，我已经有了一个网格，这就是为什么这里给出的答案示例对我没有帮助。一个更合适的例子应该有一个点和值的网格，就像我有的那样:

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)

我可以补充的另一件事是，我的网格不是正方形，而是多边形形状，是二维晶体的布里渊区。

我知道 numpy.gradient 只能在值的方形网格上正常工作，而不是我要搜索的内容。即使我将我的数据制作成一个网格，在我的原始数据的多边形之外会有很多零，这也会给我的梯度添加非常高的向量，从而影响(负面)计算的精度。这个模块在我看来更像是玩具而不是工具，恕我直言，它有严重的局限性。

使用字典解决问题。

Answer 1

您需要为 gradient 提供一个矩阵，该矩阵描述您的 (x,y) 点的角频率值。例如

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)

您可以看到将 Z 绘制为曲面会给出:

这里是解释渐变的方法:

gx 是一个矩阵，它给出了所有点的变化 dz/dx。例如gx[0][0] 是 dz/dx 在 (x0,y0)。可视化 gx 有助于理解:

因为我的数据是从 f(x,y) = sin(x+y) 生成的，所以 gy 看起来是一样的。

这是一个更明显的例子，使用 f(x,y) = sin(x)...

f(x,y)

和渐变

更新让我们看一下 xy 对。

这是我使用的代码:

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)

现在我们可以看看到底发生了什么。假设我们想知道什么点与 Z[20][30] 处的值相关联？然后……

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

重点是

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

是吗？让我们检查一下。

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

是的。

那时候我们的梯度分量是什么？

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

那些检查了吗？

dz/dy 始终为 0 检查。dz/dx = cos(x) 和...

>>> np.cos(-1.5)
0.070737201667702906

看起来不错。

你会注意到它们并不完全正确，那是因为我的 Z 数据不是连续的，步长为 0.05 并且 gradient 只能近似变化率。