从 Numpy Meshgrid 生成位置向量答案

【问题标题】：Generating Position Vectors from Numpy Meshgrid从 Numpy Meshgrid 生成位置向量
【发布时间】：2020-02-23 18:14:30
【问题描述】：

我将尝试在此处解释我的问题，而不会过多地介绍实际应用程序，以便我们能够以代码为基础。基本上，我需要对向量场进行操作。我的第一步是将字段生成为

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]) 的形状不合适，因为它是“网格网格的向量”而不是“向量的网格网格” ”。我的大问题是我似乎找不到一种方法来格式化网格网格，以使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)

【问题讨论】：

您能否发布一些更具代表性的公式（x、y 和 z 现在没有出现在其中）以及一些示例数据和所需的输出？
使用X, Y, Z = np.meshgrid(x, y, z)，我感觉A = np.array([X, Y, Z ]) 给了你一个位置向量数组。这可以直接馈送到np.cross，您可以使用axisa 或axisb 指定坐标沿哪个轴。
我没有完全听懂你的问题。 v 是与网格 x,y,z 定义的体积相关联的固定矢量场吗？您是否正在尝试计算 B(x,y,z)？
@Liris，这是我最初的直觉，但它给出了我在上面添加的错误。我的猜测是这会产生一个数组数组，numpy 假设它类似于一个无法跨越的 3x50 矩阵。
@QuangHoang，是的，在这种情况下，v 是等式中的向量 dl，我将其作为来自程序不同部分的类的 numpy 数组。真正的问题是构建 r 向量，因为我需要创建一个向量网格网格，从点 (x,y,z) 指向网格网格中的任何给定点。

标签： python arrays numpy matplotlib

【解决方案1】：

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)

【讨论】：

什么是p？ p是这个职位吗？
p 是您的参考点，即您帖子中的p= (x,y,z)。

【解决方案2】：

我认为您最初的尝试只是缺少执行叉积的轴的规范。

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]

【讨论】：