在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)

代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 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")

数值解

现在，我正在尝试通过数值计算来得到相同的结果，因此我使用这个函数来计算散度。

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")

我得到：

问题

数值解与解析解不同！我做错了什么？

- henry

Answer 1

两张图都有问题，因为您错误地使用了 np.meshgrid。

您的代码其他部分期望 xx[a, b], yy[a, b] == x[a], y[b]， 其中 a、b 是在您的情况下介于 0 和 49 之间的整数。

但是，您写成了

xx, yy = np.meshgrid(x, y)

这会导致xx[a, b], yy[a, b] == x[b], y[a]。此外，div_analy[a, b]的值变为-sin(x[b]+2y[a]) - 2cos(x[b]+2y[a])，div_num[a, b]的值也变成了其他错误的值。

您可以通过编写以下代码来解决问题：

yy, xx = np.meshgrid(x, y)

然后您将看到两种解决方案都得到了正确的结果。

顺便说一下，您设置了N = 50，这样np.linspace(xmin,xmax, N)将从区间[-2, 2]中取50个点，相邻点之间的距离将是1/49。设置N=51可能会得到更好的线性空间。

- Sora