从Numpy网格生成位置矢量

Question

我将尝试在这里解释我的问题，而不深入实际应用程序的太多细节，这样我们就可以扎根于代码。基本上，我需要对向量场进行运算。我的第一步是将字段生成为

x,y,z = np.meshgrid(np.linspace(-5,5,10),np.linspace(-5,5,10),np.linspace(-5,5,10))

请记住，这是一个广义的情况，在程序中，向量场的边界并不完全相同。总的来说，我希望能说一些类似的话

u,v,w = f(x,y,z)。你知道吗

不幸的是，这个案子需要更多的操作。我需要用一个类似于

其中向量r在程序中定义为np.array([xgrid-x,ygrid-y,zgrid-z])除以它自己的范数。基本上，这是一个从空间中的每个点指向位置（x，y，z）的向量

现在Numpy已经使用np.cross()实现了一个跨积函数，但是我似乎不能像我需要的那样创建一个“向量网格”。 我有一个lambda函数

xgrid,ygrid,zgrid=np.meshgrid(np.linspace(-5,5,10),np.linspace(-5,5,10),np.linspace(-5,5,10)) B(x,y,z) = lambda x,y,z: np.cross(v,np.array([xgrid-x,ygrid-y,zgrid-z]))

现在数组v是从另一个类导入的，看起来工作得很好，但是第二个数组np.array([xgrid-x,ygrid-y,zgrid-z])不是一个合适的形状，因为它是一个“网格向量”，而不是一个“网格向量”。我最大的问题是，我似乎找不到一种方法来格式化meshgrid，使np.cross()函数可以使用位置向量。有办法吗？你知道吗

起初，我以为我可以按照以下思路做一些事情：

x,y,z = np.meshgrid(np.linspace(-2,2,5),np.linspace(-2,2,5),np.linspace(-2,2,5)) A = np.array([x,y,z]) cross_result = np.cross(np.array(v),A)

但是，这返回了以下错误，我似乎无法避免：

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "C:\Python27\lib\site-packages\numpy\core\numeric.py", line 1682, in cross raise ValueError(msg) ValueError: incompatible dimensions for cross product (dimension must be 2 or 3)

Answer 1

我认为你最初的尝试只是缺少横积应该沿着的轴的规范。你知道吗

x, y, z = np.meshgrid(np.linspace(-2, 2, 5),np.linspace(-2, 2, 5), np.linspace(-2, 2, 5))
A = np.array([x, y, z])
cross_result = np.cross(np.array(v), A, axis=0)

我用下面的代码测试了这个。作为np.array([x, y, z])的替代方法，您还可以使用np.stack(x, y, z, axis=0)，它清楚地显示网格沿哪个轴堆叠以形成向量网格，向量与轴0对齐。我还每次打印形状，并使用随机输入进行测试。在测试中，公式的输出在一个随机指标下与输入向量在同一指标下与向量v的叉积进行比较

import numpy as np 

x, y, z = np.meshgrid(np.linspace(-5, 5, 10), np.linspace(-5, 5, 10), np.linspace(-5, 5, 10))
p = np.random.rand(3) # random reference point
A = np.array([x-p[0], y-p[1], z-p[2]]) # vectors from positions to reference
A_bis = np.stack((x-p[0], y-p[1], z-p[2]), axis=0) 
print(f"A equals A_bis? {np.allclose(A, A_bis)}") # the two methods of stacking yield the same

v = -1 + 2*np.random.rand(3) # random vector v

B = np.cross(v, A, axis=0) # cross-product for all points along correct axis
print(f"Shape of v: {v.shape}")
print(f"Shape of A: {A.shape}")
print(f"Shape of B: {B.shape}")


print("\nComparison for random locations: ")
point = np.random.randint(0, 9, 3) # generate random multi-index
a = A[:, point[0], point[1], point[2]] # look up input-vector corresponding to index
b = B[:, point[0], point[1], point[2]] # look up output-vector corresponding to index
print(f"A[:, {point[0]}, {point[1]}, {point[2]}] = {a}")
print(f"v = {v}")
print(f"Cross-product as v x a:         {np.cross(v, a)}")
print(f"Cross-product from B (= v x A): {b}")

结果输出如下所示：

A equals A_bis? True
Shape of v: (3,)
Shape of A: (3, 10, 10, 10)
Shape of B: (3, 10, 10, 10)

Comparison for random locations: 
A[:, 8, 1, 1] = [-4.03607312  3.72661831 -4.87453077]
v = [-0.90817859  0.10110274 -0.17848181]
Cross-product as v x a:         [ 0.17230515 -3.70657882 -2.97637688]
Cross-product from B (= v x A): [ 0.17230515 -3.70657882 -2.97637688]

Answer 2

有一个关于reshape和广播的解决方案：

A = np.array([x_grid, y_grid, z_grid])
# A.shape == (3,5,5,5)

def B(v, p):
    '''
    v.shape = (3,)
    p.shape = (3,) 
    '''
    shape = A.shape

    Ap = A.reshape(3,-1) - p[:,None]

    return np.cross(v[None,:], Ap.reshape(3,-1).T).reshape(shape)

print(B(v,p).shape)
# (3, 5, 5, 5)