用python计算散度 [英] Compute divergence with python

问题描述

来自

代码如下:

将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt# 边界ymin = -2.;ymax = 2.xmin = -2.;xmax = 2.# 点数 (NxN)N = 20# 发散函数定义分歧(f):num_dims = len(f)返回 np.ufunc.reduce(np.add, [np.gradient(f[i],axis=i) for i in range(num_dims)])# 创建网格x = np.linspace(xmin,xmax, N)y = np.linspace(ymin,ymax, N)xx, yy = np.meshgrid(x, y)# 定义二维向量场Fx = np.cos(xx + 2*yy)Fy = np.sin(xx - 2*yy)F = np.array([Fx, Fy])# 计算散度g = 散度(F)打印(最大:"，np.max(g.flatten()))plt.imshow(g)plt.colorbar()

创建情节:

# %%a = []对于范围内的 N(20,100):# 点数 (NxN)# = 20# 边界ymin = -2.;ymax = 2.xmin = -2.;xmax = 2.# 发散函数定义分歧(f):num_dims = len(f)返回 np.ufunc.reduce(np.add, [np.gradient(f[i],axis=i) for i in range(num_dims)])# 创建网格x = np.linspace(xmin,xmax, N)y = np.linspace(ymin,ymax, N)xx, yy = np.meshgrid(x, y)# 定义二维向量场Fx = np.cos(xx + 2*yy)Fy = np.sin(xx - 2*yy)F = np.array([Fx, Fy])# 计算散度g = 散度(F)打印(最大:"，np.max(g.flatten()))a.append(np.max(g.flatten()))plt.plot(a)

解决方案

我在

From this answer, the divergence of a numeric vector field can be computed as such:

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

However, I have noticed that the output seems to depend a lot on the grid resolution, so there seems to be something wrong!

If I look at an example:

We have the following vector field F:

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

If we compute the divergence (using Mathematica):

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

we get:

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

which has a maximum value in the range of y [-2,2] and x [-2,2]:

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

Using the divergence equation given here, we get the following plot, for max value vs. resolution (NxN: number of values in x and y-direction). None of these are even close to 3.

Here is the code:

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

Edit: To create the plot:

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

解决方案

I realized what the issue was with the help of this answer. The default spacing assumed between two consecutive values in numpy.gradient() is 1. It needs to be changed if there is a different grid.

Hence the divergence function needs to be adapted as such:

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

Example

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)

We can see that with increasing resolution, the maximum of the divergence really tends to 3.

这篇关于用python计算散度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！