用numpy初始化向量场

Question

我想初始化一个numpy数组来表示由以下定义的点组成的100 x 100网格上的二维向量场：

import numpy as np
dx = dy = 0.1
nx = ny = 100
x, y =  np.meshgrid(np.arange(0,nx*dx,dx), np.arange(0,ny*dy,dy))

该字段是关于点cx，cy的恒定速度循环，我可以用常规Python循环初始化它：

^{pr2}$

但是当我和纽比在矢量化上有困难的时候。我最接近的是

v = np.array([s * -(y-cy) / np.hypot(x-cx, y-cy), s * (x-cx) / np.hypot(x-cx, y-cy)])
v = np.rollaxis(v, 1, 0)
v = np.rollaxis(v, 2, 1)
v[np.isinf(v)] = 0

但这并不等同，也不能给出正确的答案。使用numpy初始化向量场的正确方法是什么？在

编辑：好吧-现在我很困惑下面的建议，我试着：

vx = s * -(y-cy) / np.hypot(x-cx, y-cy)
vy = s * (x-cx) / np.hypot(x-cx, y-cy)
v = np.dstack((vx, vy))
v[np.isnan(v)] = 0

但是得到一个完全不同的阵列。。。在

Answer 1

从初始设置开始：

import numpy as np

dx = dy = 0.1
nx = ny = 100
x, y = np.meshgrid(np.arange(0, nx * dx, dx),
                   np.arange(0, ny * dy, dy))
cx = cy = 5
s = 2

您可以这样计算v：

^{pr2}$

如果您真的很喜欢，可以将yx创建为一个3D数组，并在其上广播所有操作：

# we make these [2,] arrays to broadcast over the last output dimension
c = np.array([5, 5])
s = np.array([-2, 2])

# this creates a [100, 100, 2] mesh, where the last dimension corresponds
# to (y, x)
yx = np.mgrid[0:nx * dx:dx, 0:ny * dy:dy].T

yxdiff = yx - c[None, None, :]
r = np.hypot(yxdiff[..., 0], yxdiff[..., 1])[..., None]
v3 = s[None, None, :] * yxdiff / r
v3[np.isnan(v3)] = 0

检查这两个代码给出的答案是否与原始代码相同：

print np.all(v == v2), np.all(v == v3)
# True, True

编辑

为什么rx, ry = y - cx, x - cy而不是{}？我同意这是非常违反直觉的-我决定这样做的唯一原因是为了匹配您原始代码的输出。在

问题是在网格中，连续的x值实际上在x的连续列中找到，而连续的y值出现在y的连续行中，即x[:, j]是第x个x值，y[i, :]是ith y值。但是，在内部循环中，您将dx乘以i，这是您的行索引，dy乘以j，这是您的列索引。因此，您将翻转输出的x和y维。在