在 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 循环对其进行初始化：

v = np.empty((nx, ny, 2))
cx, cy = 5, 5
s = 2
for i in range(nx):
    for j in range(ny):
        rx, ry = i*dx - cx, j*dy - cy
        r = np.hypot(rx, ry)
        if r == 0:
            v[i,j] = 0,0
            continue
        # (-ry/r, rx/r): the unit vector tangent to the circle centred at (cx,cy), radius r
        v[i,j] = (s * -ry/r, s * rx/r)

但是当我无法使用 numpy 进行矢量化时。我得到的最接近的是

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：

rx, ry = y - cx, x - cy
r = np.hypot(rx, ry)
v2 = s * np.dstack((-ry, rx)) / r[..., None]
v2[np.isnan(v2)] = 0

如果您真的很喜欢，可以将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 而不是rx, ry = x - cx, y - cy？我同意这非常违反直觉——我决定这样做的唯一原因是为了匹配原始代码的输出。

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