用python计算分歧

Question

根据this 的回答，数值向量场的散度可以这样计算：

def divergence(f):
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

但是，我注意到输出似乎很大程度上取决于网格分辨率，所以似乎有问题！

如果我看一个例子：

我们有以下向量场F：

F(x) = cos(x+2y)
F(y) = sin(x-2y)

如果我们计算散度（使用 Mathematica）：

Div[{Cos[x + 2*y], Sin[x - 2*y]}, {x, y}]

我们得到：

-2 Cos[x - 2 y] - Sin[x + 2 y]

在y [-2,2]和x [-2,2]的范围内有最大值：

N[Max[Table[-2 Cos[x - 2 y] - Sin[x + 2 y], {x, -2, 2 }, {y, -2, 2}]]] = 2.938

使用此处给出的散度方程，我们得到以下图表，最大值与分辨率的关系（NxN：x 和 y 方向的值数）。这些都没有接近 3。

代码如下：

import numpy as np
import matplotlib.pyplot as plt

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

# Divergence function
def divergence(f):
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

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


# Define 2D Vector Field
Fx  = np.cos(xx + 2*yy)
Fy  = np.sin(xx - 2*yy)

F = np.array([Fx, Fy])
# Compute Divergence
g = divergence(F)

print("Max: ", np.max(g.flatten()))

plt.imshow(g)
plt.colorbar()

编辑： 创建情节：

# %%
a = []
for N in range(20,100):
    # Number of points (NxN)
    # = 20
    # Boundaries
    ymin = -2.; ymax = 2.
    xmin = -2.; xmax = 2.
    
    
    # Deivergence function
    def divergence(f):
        num_dims = len(f)
        return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])
    
  
    
    # Create Meshgrid
    x = np.linspace(xmin,xmax, N)
    y = np.linspace(ymin,ymax, N)
    xx, yy = np.meshgrid(x, y)
    
    
    # Define 2D Vector Field
    Fx  = np.cos(xx + 2*yy)
    Fy  = np.sin(xx - 2*yy)
    
    F = np.array([Fx, Fy])
    # Compute Divergence
    g = divergence(F)
    
    print("Max: ", np.max(g.flatten()))
    a.append(np.max(g.flatten()))
plt.plot(a)

Answer 1

在this answer 的帮助下，我意识到了问题所在。 numpy.gradient() 中假定的两个连续值之间的默认间距为 1。如果有不同的网格，则需要更改。

因此散度函数需要这样调整：

散度函数

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

示例

a = []
for N in range(20,100):
    # Number of points (NxN)
    # = 20
    # Boundaries
    ymin = -2.; ymax = 2.
    xmin = -2.; xmax = 2.
    
    
    # Divergence function
    def divergence(f,sp):
        num_dims = len(f)
        return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])

    
    # Create Meshgrid
    x = np.linspace(xmin,xmax, N)
    y = np.linspace(ymin,ymax, N)
    xx, yy = np.meshgrid(x, y)
    
    
    # Define 2D Vector Field
    Fx  = np.cos(xx + 2*yy)
    Fy  = np.sin(xx - 2*yy)
    
    F = np.array([Fx, Fy])
    # Compute Divergence
    sp_x = np.diff(x)[0]
    sp_y = np.diff(y)[0]
    sp = [sp_x, sp_y]
    g = divergence(F, sp)
    
    print("Max: ", np.max(g.flatten()))
    a.append(np.max(g.flatten()))
plt.plot(a)

我们可以看到，随着分辨率的增加，最大的散度确实趋于 3。