python - 用python计算分歧
问题描述
根据这个答案,数值向量场的散度可以这样计算:
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)
解决方案
在这个答案的帮助下,我意识到了问题所在。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。
推荐阅读
- javascript - PJAX 与 tawk.to
- google-cloud-platform - 用于监控谷歌云 pubsub 中未传递消息的 REST API
- vba - 仅使用 VBA 选择可见单元格的范围
- javascript - 用一个字符串替换多个相同的字符
- bash - Bash 空闲会话时间大于 15 分钟
- r - 混合设计方差分析错误:“闭包”类型的对象不是子集
- pintos - Pintos - 系统调用项目 2
- c# - 使用 TextReader 读取文件内容需要文件的完整路径 - 初学者
- java - 如何使用类型擦除来使用继承方法?
- jhipster - zuul可以配置监听80端口吗