如何根据其值计算矢量场的梯度？

Question

我想要一些关于 numpy 和数组的帮助。我想计算一个向量场的梯度。

假设我有一个函数 foo，它接受一个坐标元组 (x,y,z) 并返回一个向量 (u,v,w)。

如果我有一个坐标数组 POS = [[x1,y1,z1],[x2,y2,z2],[x3,y3,z3],etc] 我可以生成一个向量数组，其原点位于DIR 中的位置和方向 = [[u1,v1,w1],[u2,v2,w2],[u3,v3,w3],etc]。

现在我如何计算这个向量场在 POS 的每个点的梯度？我最终需要的是另一个数组 GRAD = [grad1, grad2, grad3, etc]，其中每个 grad 将是 POS 中相应点中向量场的偏导数的 3x3 数组。

PS：我知道我可以手动推导函数 foo，然后在 python 中实现导数，但在我的情况下，函数 foo 真的很复杂，我想以这种方式进行导数:)

EDIT1：现在我正在以这种方式提出 POS：

parts = 100
limit = 10
xs = linspace(-limit, limit, parts)
ys = linspace(-limit, limit, parts)
zs = linspace(-limit, limit, parts)

POS = array([(x, y, z) for z in zs for y in ys for x in xs])
DIR = array([foo(pos) for pos in POS])

如果需要，我可以这样做：

POS = POS.reshape(parts,parts,parts,3)
DIR = DIR.reshape(parts,parts,parts,3)

Answer 1

您可以通过以下方式使用numpy.gradient：

import numpy as np

N = 100
limit = .1

def vec(x,y,z): # Example vector field
    return np.array([x,x,z])

x = np.arange(-limit, limit, 2*limit/N) # np.arange takes the spacing as 3. arg
y = np.arange(-limit, limit, 2*limit/N)
z = np.arange(-limit, limit, 2*limit/N)

# Create 3D grid from 1D arrays, indexing is important!
X,Y,Z = np.meshgrid(x,y,z,indexing='ij') 
V = vec(X,Y,Z) # Get vector field, shape: (3,N,N,N)

D = np.gradient(V, x, y, z, axis=(1,2,3)) # Get gradient, this is a list!
D = np.array(D).transpose((1,0,2,3,4)) # Make it an array and flip first axes.

结果数组的形状为(3,3,N,N,N)，其中第一个轴索引向量场的分量，第二个轴索引坐标方向w.r.t。导数是计算出来的。