python - 在 Python 中正确计算向量场的散度

Question

我正在尝试计算向量场的散度:

Fx  = np.cos(xx + 2*yy)
Fy  = np.sin(xx - 2*yy)
F = np.array([Fx, Fy])

分析解决方案
基于散度的解析计算，散度 (div(F) = dF/dx + dF/dy ) 应如下所示(参见 Wolfram Alpha here):

dFx/dx = d/dx cos(x+2y) = -sin(x+2y)

dFy/dy = d/dy sin(x-2y) = -2*cos(x-2y)

分歧:

div_analy = -np.sin(xx + 2*yy) - 2*np.cos(xx - 2*yy)

编码:

# Number of points (NxN)
N = 50
# Boundaries
ymin = -2.; ymax = 2.
xmin = -2.; xmax = 2.


# Create Meshgrid
x = np.linspace(xmin,xmax, N)
y = np.linspace(ymin,ymax, N)
xx, yy = np.meshgrid(x, y)

# Analytic computation of the divergence (EXACT) 
div_analy = -np.sin(xx + 2*yy) - 2*np.cos(xx - 2*yy)

# PLOT
plt.imshow(div_analy , extent=[xmin,xmax,ymin,ymax],  origin="lower", cmap="jet")

数值解
现在，我试图获得相同的数字，所以我使用了 this function计算散度

def divergence(f,sp):
    """ Computes divergence of vector field 
    f: array -> vector field components [Fx,Fy,Fz,...]
    sp: array -> spacing between points in respecitve directions [spx, spy,spz,...]
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])

当我使用此函数绘制散度图时:

 # Compute Divergence
 points = [x,y]
 sp = [np.diff(p)[0] for p in points]
 div_num = divergence(F, sp)
    
 # PLOT
 plt.imshow(div_num, extent=[xmin,xmax,ymin,ymax], origin="lower",  cmap="jet")

......我明白了:

问题
数值解不同于解析解!我究竟做错了什么？

Answer 1

让我 1. 解释这种观察背后的原因，以及 2. 如何解决它。
原因:
在计算散度(或一般的梯度)时，需要注意数据的方向，因为沿着正确的轴计算梯度以获得物理上有效的结果很重要。
np.meshgrid 可以通过两种方式输出网格，具体取决于您如何设置 索引参数
索引“xy”:这里，对于每个 y 值，我们扫描 x 值。

for y in ys:
   for x in xs:

索引“ij”:这里，对于每个 x 值，我们扫描 y 值。

for x in xs:
   for y in ys:

最后:
因此，给定特定索引，结果向量场的“x”或“y”维度将与 numpy 数组(向量场)的第一个轴相关联。

对于“xy”索引:第一个轴沿 x 方向

对于“ij”索引:第一个轴沿 y 方向

所以，基本上，你在计算:dFx/dy + dFy/dx而不是 dFx/dx + dFy/dy如何修复:
我重新制定了散度方程以考虑到这一点:

def divergence(f, sp, indexing = "xy"):
    """ 
    Computes divergence of vector field 
    f: array -> vector field components [Fx,Fy,Fz,...]
    sp: array -> spacing between points in respecitve directions [spx, spy,spz,...]
    indexing: "xy" or "ij", see np.meshgrid indexing 

    """
    num_dims = len(f)
    
    if indexing == "xy":
        return np.ufunc.reduce(np.add, [np.gradient(f[num_dims - i - 1], sp[i], axis=i) for i in range(num_dims)])
    if indexing == "ij":
        return np.ufunc.reduce(np.add, [np.gradient(vectors[i], sp[i], axis=i) for i in range(num_dims)])
    

使用这个等式，我们得到: