Python计算向量场散度时数值解与解析解不一致问题排查

错误原因 #

核心问题是np.meshgrid生成的网格轴顺序和你散度函数中指定的求导维度不匹配：

• np.meshgrid(x,y)默认采用indexing='xy'模式，生成的xx、yy数组形状为(N_y, N_x)，第0轴对应y坐标变化，第1轴对应x坐标变化。
• 你的原divergence函数对第i个向量分量直接沿axis=i求梯度，相当于对x方向分量Fx沿y轴求导、y方向分量Fy沿x轴求导，完全颠倒了求导方向，自然得到的结果和解析解差异极大。

修正方案 #

你可以直接调整散度函数的求导轴适配meshgrid的默认索引规则，修正后的代码如下：

def divergence(f,sp):
    """ Computes divergence of 2D vector field 
    f: array -> vector field components [Fx,Fy]
    sp: array -> spacing between points [spx, spy]
    """
    # Fx沿x方向求导对应axis=1，Fy沿y方向求导对应axis=0
    dFx_dx = np.gradient(f[0], sp[0], axis=1)
    dFy_dy = np.gradient(f[1], sp[1], axis=0)
    return dFx_dx + dFy_dy

如果你需要保留原函数的多维度通用性，也可以在生成网格时指定indexing='ij'：

xx, yy = np.meshgrid(x, y, indexing='ij')

修改后数值解和解析解的误差会非常小，仅在边界处因为梯度计算的边界处理规则存在可忽略的微小差异。

内容的提问来源于stack exchange，提问作者henry