使用Python进行梯度计算

Question

使用Python进行梯度计算

10

我想知道numpy.gradient是如何工作的。我使用梯度来尝试计算群速度（波包的群速度是频率对波数的导数，而不是一组速度）。我将一个三列数组输入给它，前两列是x和y坐标，第三列是该点（x，y）的频率。我需要计算梯度，并且我期望得到一个二维向量，因为梯度的定义如此。

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

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

df/dx*i+df/dy*j

但是我有两个每个有3列的数组，即2个3D向量；起初我以为这两个向量的和会给我所寻找的向量，但是z分量并没有消失。希望我在解释上已经足够清楚了。我想知道numpy.gradient是如何工作的，以及它是否适合我的问题。否则，我想知道是否有其他可以使用的Python函数。

我的意思是：我想计算一个值数组的梯度。

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

在一个均匀网格上，x1和x2是点的坐标（我在布里渊区域上的点），x3是该点频率的值。我还输入了两个方向的导数步骤。

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只能在一个值的正方形网格上正常工作，而不是我所寻找的。即使我将我的数据制作成一个网格，在原始数据的多边形之外会有很多零值，这将给我的梯度添加非常高的向量，从而对计算的精度产生负面影响。在我看来，这个模块更像是一个玩具而不是一个工具，它有严重的限制。在Python中是否有更强大的东西，或者我最好从头开始编写其他东西，也许转向Fortran？

- Avulso Malloru

1

那么问题是什么？你应该使用哪个模块？有什么出错了吗？ - Stephan

问题是梯度做什么？为什么给我两个3D向量而不是1个2D向量？梯度实际上是否真正计算了梯度？从它的输出来看，我不能确定。在我的看来并不精确。 - Avulso Malloru

我认为很清楚，我的输入的第三个组件是标量场，第三个组件上的每个值都是每个（x，y）点的函数值。 - Avulso Malloru

问题在于您给梯度提供了错误的输入。它不关心 x1，x2 或最后一个示例中的 i，j。它只想要一个 i ** 2 + j ** 2 值的矩阵。您的 i ** 2 + j ** 2 值的矩阵隐含地对应于 xy 平面，并且 gradient 的可选标量参数考虑了步长假设，即如果您的 x 点之间不是相距 1，您的 y 点也是如此。我今天要出城，但我晚上回来后会更新我的答案。 - seth

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- seth · Accepted Answer

你需要给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作为曲面绘制，可以看到：

sinxpy

下面是如何解释梯度的：

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

由于我的数据是从 f(x,y) = sin(x+y) 生成的，因此 gy 看起来也是一样的。

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

f(x,y) enter image description here

和梯度

更新我们来看看 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只能近似表示变化率。